Complex Oscillatory Waves Emerging from Cortical Organoids Model Early Human Brain Network Development

Abstract

Structural and transcriptional changes during early brain maturation follow fixed developmental programs defined by genetics. However, whether this is true for functional network activity remains unknown, primarily due to experimental inaccessibility of the initial stages of the living human brain. Here, we developed human cortical organoids that dynamically change cellular populations during maturation and exhibited consistent increases in electrical activity over the span of several months. The spontaneous network formation displayed periodic and regular oscillatory events that were dependent on glutamatergic and GABAergic signaling. The oscillatory activity transitioned to more spatiotemporally irregular patterns, and synchronous network events resembled features similar to those observed in preterm human electroencephalography. These results show that the development of structured network activity in a human neocortex model may follow stable genetic programming. Our approach provides opportunities for investigating and manipulating the role of network activity in the developing human cortex.

Keywords: cortical organoids; network oscillations; phase-amplitude coupling; preterm electroencephalography; single-cell transcriptomics; stem cells.

Figure 1.. Cellular and molecular development of human cortical organoids.

(A) Overview of human neural network formation and dynamics evaluation using cortical organoids. (B) Schematic of the protocol used to generate cortical organoids. Scale bar, 200 μm. (C) Organoid growth during different developmental stages. (D) Representative immunostainings showing proliferating NPCs (Ki67 and Nestin), lower (TBR1 and CTIP2) and upper (SATB2) cortical layer neurons, glial cells (GFAP) and GABAergic (CR) neurons overtime. Scale bar, 50 μm. (E) Population analysis of specific markers indicating stages of maturation and multiple neuronal subtypes. The data are shown as mean ± s.e.m. (n = 8). (F) Representative image of a pyramidal neuron; dendritic structures are observed in cells transduced with the SYN:EGFP reporter (scale bar, 5 μm). (G) Electron microscopy of synaptic ultrastructure in 4-month cortical organoids (blue). Scale bar, 200 nm. (H) Uniform manifold approximation and projection (UMAP) plot of 15,990 cells from integrating datasets on 1-, 3-, 6-, and 10-month cortical organoids. Colors denote cells sampled from four different time points. (I) UMAP plot of the integrated datasets colored by seven main cell clusters. Red as GABAergic neurons, orange as glutamatergic neurons, blue as glia cells, green as intermediate progenitors, purple as progenitors, green blue as mitotic cells, and grey as others. (J) Separate UMAP plots of integrated data by different time points. Same color scheme used for main cell clusters. (K) Dot plots showing cluster specific gene expression across main cell clusters. (L) UMAP plots showing expression levels of cell-type specific markers (see Figure S1 for additional markers). (M) Barplots of proportion of cell types at individual time points.

Figure 2.. Oscillatory network dynamics in long-term cortical organoids.

(A) Schematic of the organoid signal processing pipeline. Raw MEA data is analyzed as population spiking and LFP separately. Synchronous network events are highlighted in yellow. (B) Raster plot of network spiking activity after 1.5 and 6 months of maturation. A 3-s interval of activity over 5 channels is shown in the upper right corners. (C) Cortical organoids show elevated and continuously increasing mean firing rate compared to 2D monolayer neurons (n = 8 for organoid cultures, and n = 12 for 2D neurons). Inset, correlation of the firing rate vector over 12 weeks of differentiation (from 8 to 20) between pairs of cultures showing reduced variability among organoid replicates. (D) Temporal evolution of cortical organoid network activity. Detailed definitions and further parameters are presented in Figure S2. (E) Time series of population spiking and LFP during network events in cortical organoid development. Each overlaid trace represents a single event during the same recording session. The number of subpeaks during an event (F) is maximized and stereotypical at 6-months, developing nonlinearly and following an inverted-U trajectory. (G) Network variability, measured as the coefficient of variation of the inter-event interval, increases monotonically throughout development. (H) 1-4 Hz oscillatory power in the LFP increases up to the 25th week in culture and plateaus at 30 weeks. Inset, oscillatory power is calculated by fitting a straight line (dashed) over the aperiodic portion of the PSD and taken as the height of narrow peaks rising above the linear fit. (I) Pairwise correlation of LFP across all electrodes (coherence) within a well during network events initially increase, then decreases after 30 weeks. (J) An example of sequential frames during a network event show the spatial propagation of wave spreading, then disappearing again after 100ms. The data shown in C, D, F, G, H, and I are presented as mean ± s.e.m. *P < 0.05, **P < 0.01, ***P < 0.001, unpaired Student’s t-test (C), quadratic (F, H and I) and linear (G) regression.

Figure 3.. Cortical organoid serves as a model of functional oscillations and their synaptic mechanisms.

(A-C) Phase-amplitude coupling is observed in organoid LFP during network events. (A) Example of raw LFP during a network event decomposed into its low-frequency component (1-4 Hz delta) and the amplitude envelope of the high-frequency, broadband gamma component (200-400 Hz). Analysis was repeated for 100-200 Hz with near identical effect size and significance. (B) Normalized gamma amplitude binned by delta phase during network events (black) shows greater modulation depth by low frequency delta than during non-event periods (red). (C) Phase-amplitude coupling during network events is significantly greater than non-event periods in all batches. (D) Effect of selective drug treatments on neuronal electrical activity in 6-month organoids. Representative raster plots and burst measurements of untreated and treated organoids. The pharmacological manipulation was performed using cortical organoid plated on 4 MEA wells (n = 4, cortical organoid culture for each treatment). Scale bar, 20 s. Exposure to AP5 + CNQX, baclofen and muscimol reversibly extinguish the network bursts (synchrony), while no changes were promoted by bicuculline. (E-F) Pharmacological perturbation of oscillatory activity during network events in 6-month organoids. Pre and post refer to before treatment administration and after administration, respectively. Application of bicuculline and picrotoxin increases the number of network events, while CNQX + AP5 and baclofen completely abolish synchronized network events. Bicuculline blocks oscillatory network activity but not the network event itself. Data are shown as mean ± s.e.m.; unpaired Student’s t-test.

Figure 4.. Cortical organoid network dynamics mimic premature neonates after 28 weeks of maturation.

(A) Representative LFP trace from cortical organoid, highlighting instances of network events (yellow). Comparable events between periods of quiescence (discontinuous network dynamics) are shown in human preterm neonate EEG at 35 weeks gestational age, while a different pattern of continuous activity is observed in adult EEG. SAT: spontaneous activity transient. (B) Examples of analogous features in preterm neonate EEG and organoid LFP show various levels of similarity throughout development. RMS: root mean square; 50%/95% refer to 50th and 95th percentile of feature distribution within a recording. (C) Pearson’s correlation coefficients between age and electrophysiological features (mean ± std of bootstrapped distribution) for both organoid and premature neonates show different degrees of developmental similarity for individual features (12 total selected). For example, SATs (events) per hour shows remarkable similarity over time between organoid and neonates. (D) Schematic of machine learning procedure for age prediction: EEG features from 39 premature neonates (n = 567 recordings) between 25- and 38-weeks PMA (post-menstrual age) were used to train and cross-validate a regularized regression model (ElasticNet), optimizing for preterm neonate age prediction based on their EEG features only (top). The model was clamped after training, and applied directly on organoid LFP features and control datasets, including held-out preterm neonate data, mouse primary culture, 2D iPSC culture, and human fetal brain culture. (E) Model-predicted developmental time (y-axis, age in weeks) follows actual weeks-in-culture (x-axis) for organoids (orange and blue), as well as true age of held-out preterm neonate data points (black). Dashed line represents unity, signifying perfect prediction. Large circles on solid lines and shaded regions denote mean ± std of prediction, respectively, while dots indicate per-sample prediction (n = 8 for organoids at all time points). (F) Pearson’s correlation coefficient between predicted and actual developmental time for organoid, and control datasets. Significant positive correlations indicate the model’s ability to capture developmental trajectory in a particular dataset.