Human pluripotent stem cell-derived neural constructs for predicting neural toxicity

Results

Cell types representing distinct components of the developing brain were derived from the H1 human ES cell line (SI Appendix, Fig. S1) (21), including neural progenitor cells (24), endothelial cells, mesenchymal stem cells (25), and microglia/macrophage precursors (SI Appendix, SI Materials and Methods). Vascular and microglial components were incorporated into the neural constructs to mimic the repertoire of cell types and cell–cell interactions of the developing human brain that might be susceptible to toxic exposure (16⇓–18). The vascular component included both endothelial cells and mesenchymal stem cells, because mesenchymal support cells play an important role in the maintenance and long-term stability of vascular networks in vitro and in vivo (26). Neural tissue constructs were produced on peptide functionalized polyethylene glycol (PEG) hydrogels that were formed in 24-well Transwell inserts using “thiol-ene” photopolymerization (SI Appendix, Fig. S1) (27, 28). For these experiments, the time point when neural progenitor cells were seeded on the PEG hydrogels was defined as day 0. To mimic the in vivo recruitment of blood vessels and microglia after the initial formation of the neural tube (10⇓⇓⇓⇓–15), vascular cells and microglia/macrophage precursors were added at day 9 and day 13, respectively (SI Appendix, Fig. S1).

RNA-Seq was used to quantitatively assess sample reproducibility and to identify global gene expression patterns within the neural constructs. Replicate samples (n = 4) were characterized by Spearman’s correlation coefficients (ρ) ≥ 0.97 at days 16 and 21 of growth on PEG hydrogels (Dataset S1) (29). Spearman’s rank comparisons to RNA-Seq data for human samples from the Allen Brain Atlas (30, 31) demonstrated that the neural constructs were most correlated to early developmental time points and least correlated to later adult time points (Dataset S1). For example, Spearman’s coefficients for day 16 and day 21 neural constructs were higher for all comparisons to eight postconception week (PCW) samples ([ρ] ≥ 0.82, all brain regions) than 30-y-old adult samples ([ρ] ≤ 0.76, all brain regions) (30, 31). However, given the timing of our differentiation protocol, the neural constructs likely represent developmental time points before the earliest available RNA-Seq data (8 PCW) from the Allen Brain Atlas (30, 31).

RNA-Seq data were then analyzed by EBSeq (32) to identify genes up-regulated within the neural constructs compared with undifferentiated human ES cells (Dataset S2). Characteristic gene ontology (GO) clusters were identified from the resulting gene sets using the DAVID Bioinformatics Database Functional Annotation Tool (Dataset S2) (33, 34). Genes threefold up-regulated with an EBSeq false discovery rate (FDR) ≤ 0.005 for day 21 neural constructs relative to H1 ES cells were enriched within GO categories that included neurogenesis (GO:0022008, 206 genes), forebrain development (GO:0030900, 40 genes), hindbrain development (GO:0030902, 26 genes), synaptic transmission (GO:0007268, 112 genes), and vasculature development (GO:0001944, 61 genes) (Dataset S2). RNA-Seq also identified expressed genes for phenotypes important to neurogenesis (Datasets S2 and S3), such as GABAergic neurons (e.g., GABA receptors), glutamatergic neurons (e.g., VGLUT2 and VGAT), cortical neurons (reelin/RELN, BRN2/POU3F2, CTIP2/BCL11B, etc.), synaptic markers (e.g., synapsins and synaptic vesicle components), and glial cells (GFAP, PDGFRA, GLAST/SLC1A3, etc.) (9, 35⇓⇓⇓–39).

Immunofluorescence imaging was used to investigate cellular organization within the neural constructs. Neural progenitor cells differentiated and self-assembled into layered βIII-tubulin⁺ and GFAP⁺ cells that extended around the circumference of the neural constructs by day 9 of culture on PEG hydrogels (SI Appendix, Fig. S2). By day 21, sample thickness varied from ∼50 μm in the center to >350 μm in outer regions of the neural constructs (Fig. 1), where complex cellular organization was most pronounced. Ki67⁺ proliferative layers (SI Appendix, Fig. S3) and radially organized cells that expressed GFAP, phospho-vimentin (p-VIM), and/or SOX2 were each identified within the neural constructs (SI Appendix, Fig. S4) (9, 38, 39). Immunofluorescence imaging revealed markers for GABAergic neurons (GABA; Fig. 2A) and glutamatergic neurons (VGLUT2; Fig. 2B), and reelin⁺ (RELN, Fig. 2C), calretinin⁺ (CALB2, Fig. 2C), Brn2⁺ (POU3F2, Fig. 2D), Ctip2⁺ (BCL11B, Fig. 2E), and Tbr1⁺ layers were also evident (Fig. 2E) (9, 35⇓–37). Thus, neural phenotypes identified by immunofluorescence imaging were in agreement with the RNA-Seq results (Datasets S2 and S3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Morphological characteristics of neural constructs. Human embryonic stem cell-derived precursor cells were cocultured on polyethylene glycol (PEG) hydrogels in 24-well Transwell inserts. Neural progenitor cells (NPCs) were seeded on synthetic PEG hydrogels (day 0), followed by endothelial cells (ECs) and mesenchymal stem cells (MSCs) at day 9 and microglia/macrophage precursors (MGs) at day 13 (SI Appendix, Fig. S1). (A and B) Maximum projection Z stack (525-μm thickness) and slice views (NIS Elements) illustrating βIII-tubulin (green), GFAP (red), and DAPI (blue) for a day 21 neural construct. XZ and YZ cross-sections are illustrated in the regions indicated by dashed lines. The boxed region in A is illustrated in B. (C and D) Volume view images (NIS Elements) corresponding to (C) the full neural construct shown in A (6,300 μm × 6,300 μm × 550 μm) and (D) the region shown in B (1,570 μm × 2,290 μm × 300 μm). (Scale bar in A, 1,000 μm and B, 500 μm.)

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Neural differentiation within the neural constructs. Cryosectioned samples illustrating neuronal phenotypes after 21 d of growth on PEG hydrogels: (A) βIII-tubulin (red) and GABA (green), (B) βIII-tubulin (red) and VGLUT2 (green), (C) reelin (red) and calretinin (green), (D) MAP2 (red) and Brn2 (green), and (E) Ctip2 (green) and Tbr1 (red). (Scale bars in A–C, 50 μm and D and E, 100 μm.)

Human ES cell-derived endothelial cells and mesenchymal stem cells were added to the neural constructs after neural progenitor cells had self-assembled into multilayers with radially organized neural and glial populations (SI Appendix, Fig. S2) (9, 37, 38) to mimic the recruitment of blood vessels by the neuroepithelium (15). The neural constructs were characterized by extensive capillary networks by day 21 (Fig. 3A; SI Appendix, Fig. S5), which associated with radially oriented glial cells (Fig. 3 B–D) and penetrated into neuronal layers (SI Appendix, Fig. S5) (13⇓–15). RNA-Seq demonstrated that the neural constructs expressed many vascular genes [day 21: 122 genes, GO:0001944 vasculature development, transcripts per million (TPM) > 16; Dataset S4], including blood-vessel–promoting growth factors (e.g., VEGFA) that were also detected in control samples lacking vascular cells (Dataset S3). Because factors such as VEGF were not added exogenously, cellular signaling within the neural constructs provided the necessary cues for the formation of interconnected vascular networks, which is consistent with the recruitment of capillaries by cells within the neuroepithelium (15).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Vascular network formation within the neural constructs. (A and B) Immunofluorescence for endothelial cells (CD31, green), glial cells (GFAP, red), and nuclei (DAPI, blue) for a day 21 neural construct. (B) Zoom of the boxed region shown in A to illustrate association and alignment for a capillary tubule and radially oriented glial cells (arrows). The cells in B are shown as single channel grayscale images for (C) CD31 and (D) GFAP. [Scale bars in A, 250 μm and B–D, 100 μm (shown in B).]

Microglia/macrophage precursors were derived by differentiating H1 ES cells through mesendoderm and hemogenic endothelium lineages (40), which resemble early precursors in the yolk sac that contribute to microglia in vivo (41). The microglia/macrophage precursors were CD11b⁺CD14⁺ by FACS analysis (SI Appendix, Fig. S1) and phagocytosed yeast particles (Movie S1). Iba1⁺ microglia begin to populate the human brain as early as 4.5 gestational weeks (GW), but minimally interact with early blood vessels (11). Therefore, the microglia/macrophage precursors were added to the neural constructs after initial vascular network assembly (11). Characteristic microglia genes were up-regulated relative to control samples (Fig. 4A) and were expressed at days 16 and 21 (SI Appendix, Fig. S6 and Dataset S3) when microglia/macrophage precursors were incorporated into the neural constructs (10⇓–12). TREM2, IBA1/AIF1, and other microglial genes were detectable by RNA-Seq only when microglia/macrophage precursors were incorporated into the neural constructs (Fig. 4A). By day 21, Iba1⁺ cells were dispersed throughout the neural constructs (SI Appendix, Fig. S6), associated with endothelial tubules (Fig. 4B and SI Appendix, Fig. S6), and adopted both amoeboid and ramified morphologies (Fig. 4B and SI Appendix, Fig. S6) (10⇓–12).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Incorporation of microglia into neural constructs. (A) Gene expression for neural constructs with or without microglia (Dataset S3, Quality Control Experiments; N.D., not detected). Statistical analysis was conducted using a Student’s t test (TPM ± SD; ***P < 0.001; n = 4 replicate samples each). (B) Iba1 (microglia, red) and CD31 (endothelial cells, green) expression for a day 21 neural construct. Microglia adopt ramified morphologies (e.g., closed arrow) and associate with capillary tubules (e.g., open arrows). (Inset) Iba1 (red) and DAPI (blue) expression for the cell pointed out by the closed arrow (Bottom, Right corner) and surrounding nuclei. Image is brightened for clarity. (Scale bar, 100 μm.)

RNA-Seq and linear support vector machines were then used to build a predictive model for neurotoxicity based on changes in global gene expression by neural constructs exposed to known toxins and nontoxic controls (Fig. 5 and SI Appendix, Fig. S7 and Datasets S5 and S6) (42⇓⇓⇓–46). Neurotoxicity was evaluated using a set of 31 control compounds and 39 toxins with previous literature support for toxicity (Dataset S5). Control chemicals included common food additives or pharmaceuticals with no known neurotoxicity. Dosage was based on previous in vitro or in vivo data for toxic and nontoxic compounds, which included blood serum concentrations (such as for pharmaceuticals) or previously published screening conditions when available or set at 10 μM when data were not available. Replicate day 14 neural constructs were continuously exposed to test compounds through day 16 (2 d of exposure) or day 21 (7 d of exposure), and then harvested for RNA-Seq (Datasets S7 and S8) and machine learning analysis (Dataset S6).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Machine learning predictions. (A) A linear support vector machine (SVM) for a 2D problem, where an (n − 1)-dimensional hyperplane reduces to a line that separates the classes (filled vs. open circles) and maximizes the closest points between classes (the support vectors, which fix the position and orientation of the hyperplane). The x_is are the examples (points in A), the y_is are their labels (filled or open in A), and w is the weight vector, or vector of coefficients on the features (the dimensions). The linear SVM’s output is the weight vector w and the other coefficient b. To make a prediction, the SVM computes the number w′x_i − b, and outputs the label 0 (nontoxic, for our application) if this number is less than 0, and 1 otherwise. The extensions required for the soft margin version of the SVM are highlighted in pink in the equation, which minimizes the sum of the distances between incorrectly classified training points (ξ_i) in addition to the margin, and is used when the data are not linearly separable (43). (B) The ROC curves are shown for the data averaged from days 16 (2-day dosing) and 21 (7-day dosing) (Dataset S6). The ROC curve plots true positive rate on the y axis against the false positive rate (1 − specificity) on the x axis as the threshold is varied (SI Appendix, SI Materials and Methods shows additional details).

We used two standard hold-out testing methods for evaluation to avoid overly optimistic prediction of accuracy (45⇓–47): (i) A nearly unbiased (slightly pessimistic) leave-one-out cross-validation and (ii) an unbiased blinded trial with a single hold-out set. For leave-one-out cross-validation, there were 60 compounds in the training set and the method proceeded in 60 steps. In each step a different data point was held out of the training set, the support vector machine was trained on the remaining data points, and a prediction was made for the held-aside data point. Hence every data point was a test case exactly once, for a model trained without that data point. Results were aggregated over all of the folds, or test cases, to estimate how well the support vector machine trained on all data will perform on a new data point (compound). Performance estimates are shown in the form of receiver operating characteristic (ROC) curves (Fig. 5 and SI Appendix, Fig. S7). For the second method, the unbiased blinded trial used the predictive model generated from the training set to make predictions on a separate hold-out set, including estimates of accuracy and area under the ROC curve (AUC). The leave-one-out cross-validation method has lower variance than a single train/test split because it tests on all of the compounds of the training set, but it is a slightly pessimistic estimate of future performance because each training set is slightly smaller (one less) than the actual training set.

Leave-one-out cross-validation was used to evaluate neural constructs exposed to a training set of 34 toxins and 26 nontoxic controls (Datasets S5 and S6). The area under the ROC curves for the training compounds were 0.86 on day 16 (SI Appendix, Fig. S7), 0.88 on day 21 (SI Appendix, Fig. S7), and 0.91 for data averaged from both days (Fig. 5B). Thus, the support vector machine produced an estimate of future data ≥0.86 for each day individually and 0.91 using data from two developmental time points. Additionally, a wrapper-based approach of recursive feature elimination (48) demonstrated that the number of genes could be reduced to 1,000 without harming accuracy (Dataset S6). Unbiased hold-out testing was then used to predict toxicity for a set of 10 blinded compounds that were not in the initial training set (five toxins, five nontoxic controls, Dataset S5) and were unknown to researchers generating the support vector machine model until after the predictions were made. The support vector machine model was generated using the full set of ∼19,000 genes and data were averaged from days 16 and 21 training sets to make predictions, which produced probabilities for ranking the blinded compounds from most likely toxic to least likely toxic (Dataset S6). In addition, we used a threshold of 0.5 to make definitive predictions, labeling every chemical with probability ≤0.5 as nontoxic and all others toxic. The area under the ROC curve generated for the ranking of the blinded set was 0.92, which is in agreement with the training experiment. Importantly, all compounds except oleic acid (a false positive) were properly assigned as toxic or nontoxic (Dataset S6). Therefore, the predictive model correctly classified 9 of 10 compounds in the blinded trial.