Neural organoids were differentiated from iPSC-derived Neural Progenitor Cells (NPC) for up to 14 weeks and characterized throughout development (Fig. 1a). Gene expression of synaptic plasticity markers was quantified from week 0 to week 12. Calcium signaling development was analyzed from week 2 to week 14. Finally, electrical activity was characterized by High-Density Microelectrode Arrays (HD-MEAs) over two independent time periods, from weeks 6-to-9 and 10-to-13. In addition, pharmacological modulation of neurotransmission was performed at weeks 8 and 13. Lastly, input-specific Theta Burst Stimulation (TBS) was implemented at week 14 to induce synaptic plasticity. To analyze both spontaneous and evoked electrical activity from the HD-MEA data, functional connectivity and criticality analysis were performed. A schematic overview of the neurocomputational investigations is shown in Fig. 1b. In addition, an example of how evoked activity from pharmacological or electrical stimuli can modulate synaptic transmission to induce synaptic plasticity is shown in Fig. 1c. In long-term potentiation, IEGs play a role in trafficking of glutamatergic receptors into the postsynaptic terminal (Fig. 1c), therefore expression of these IEGs can serve as a molecular marker for long-term memory27,37.
Neural organoids develop proper synapse formation and express receptors necessary for synaptic transmission
Neural organoids were differentiated following our in-house established protocol11. The expression of markers for astrocytes (GFAP), oligodendrocytes (MBP, OLIG2) and mature neurons (MAP2) increased in the first 8 weeks of maturation and then plateaued in the following weeks, indicating that the differentiation predominantly occurs rapidly until week 8 and then reaches a more stable, mature cell composition from week 8 to 12 (Supplementary Fig. 1). Hence, two time points were selected for the experiments in this paper.
RNA-sequencing (RNAseq) revealed some interesting trends in gene expression between week 8 and week 12/13 organoids. Firstly, while many GABAergic markers showed no difference in expression (e.g., PVALB, GABRA1, LHX6), we observed a slight trend in downregulation of several GABAergic markers, such as SST, SLC32A1, and GAD 1 and 2 in the week 12/13 organoids compared to the week 8 group, as shown in Supplementary Fig. 2a and b.
By comparing expression levels across main cell lineages and brain regions (Supplementary Figs. 3 and 4), we observe an increase in the expression of astrocyte markers (e.g, AQP4, GFAP, HOPX, S100B, and SLC1A3) in the week 12/13 group compared to week 8. These data suggest that the observed trends in downregulation in expression of GABA and glutamatergic markers (GRIN3A, GRIN3B, and SLC17A7) may be linked to the shift in cellular populations with more astrocytes (based on the higher expression of astrocyte markers in the week 12/13 organoids), (Supplementary Figs. 2 and 3). We also profiled the expression of brain region markers, which showed the presence of genes across the forebrain, hindbrain, and midbrain and no consistent differences in expression between the two age groups (Supplementary Fig. 4).
These findings are further supported by the functional enrichment analysis (Supplementary Fig. 5), which revealed that genes associated with synaptic signaling, neuronal differentiation, and axonogenesis were downregulated and cell cycle genes were upregulated in the week 12/13 group (Supplementary Fig. 5).
We then evaluated the presence of pre- and postsynaptic markers as well as IEGs as the key proteins involved in synaptic plasticity and learning formation. Presence of the presynaptic marker Synaptophysin (SYP) and postsynaptic marker HOMER1 was detected in both week 8 and 12 organoids, displaying typical punctual staining (Fig. 2a). Gephyrin-positive signal was close to background with few positive cells at week 8 and increased at week 12 (Fig. 2b). This indicates that there are more developed inhibitory synapses at the later stage of differentiation. Gene expression of GABRA1, which encodes the inhibitory GABAA receptor, followed the same pattern (Fig. 2c). Gene expression of postsynaptic marker HOMER1 was steady over time (Fig. 2c). Presence of Synapsin1 and MAP2 seemed to be reduced at week 12 vs. week 8, supporting RNAseq data of the lineage shift towards more astrocytes (Fig. 2b), although cell number quantification is needed to drive the final conclusions.
Both AMPA and NMDA receptors play an important role in synaptic plasticity, including STP/LTP24,38,39, therefore showing expression of these receptors was imperative for this study to give insight into the mechanisms of learning and memory in neural organoids. The increase in gene expression over time was the greatest for GRIN1, which plateaued around week 8 to week 12 (Fig. 2c). GRIN2A and GRIN2B both increased over time with a higher increase of GRIN2A expression than GRIN2B, suggesting increasing maturity40(Fig. 2c). GRIA1 expression also increased over time and plateaued after week 8 (Fig. 2c). Thus, plateau in expression of these subunits suggests the organoids reached a mature state between week 8 to 12.
Dynamic expression of immediate early genes associated with synaptic plasticity and cognitive functions over time
IEGs are crucial for cognitive functions as they act directly at the synapse and mediate the cellular processes that are essential for learning and memory consolidation27. To demonstrate that the neural organoids have the necessary cellular components for cognitive processes, we quantified IEG expression during the course of differentiation (Fig. 2d, e, Supplementary Fig. 6). Gene expression of ARC, BDNF, NPAS4, NPTX2, and FOS was significantly increased over time, while EGR1 was already expressed in NPCs and remained at levels close to those in NPCs. Expression of upstream regulators of IEGs (CREB and CAMK2A) also increased over time with the largest increase in expression of CAMK2A (Fig. 2f). In addition, SYNGAP1, which plays a key role in regulating synaptic function and plasticity41, was stably expressed throughout the course of differentiation, starting from NPCs. The levels of IEG proteins (NPTX2, ARC, and BDNF) and upstream IEG transcription factor CREB were comparable between week 8 and 12, confirming the plateau observed in RT-qPCR data (Fig. 2d). These results were consistent with RNAseq data which showed no changes in IEG expression between 8 and 12/13 week old organoids (Supplementary Fig. 6a and Supplementary Fig. 6b).
Finally, we assessed the expression of microRNAs known to be involved in synaptic plasticity (Fig. 2g)42 and observed expected dynamics. A strong increase in expression of mir-124-3p over the course of differentiation was observed. mir-132-3p and mir-134-3p had opposite expression patterns: mir-132-3p was increased over time while expression of mir-134-3p was first strongly induced from NPC to 2 weeks of differentiation and was downregulated thereafter (Fig. 2g).
Evidence of spontaneous electrical activity and highly interconnected neuronal networks in neural organoids
Electrophysiology over the course of organoid development was characterized using calcium imaging and HD-MEAs. In addition to the expression of molecular machinery involved in synaptic plasticity, neural organoids showed highly patterned spontaneous electrical activity (Figs. 3 and 4). Calcium transients were measured using Fluo-4 biweekly from week 2 to 14. Whole organoid change in fluorescence over resting fluorescence intensity (∆F/F) was quantified and compared across age groups (Fig. 3a). From these ∆F/F plots, the average rise time, peak amplitude, firing rate, decay time, burst duration, number of peaks, and percentage of active organoids per time point were calculated (Fig. 3b). Bursts were identified as peaks in calcium transients. Burst firing rate was calculated as the number of burst peaks per second.
Calcium imaging showed that 2-week-old organoids did not exhibit spontaneous oscillatory calcium dynamics. The first signs of oscillatory calcium activity were detected at week 4, with high-frequency oscillations at weeks 4 and 6, as shown by higher burst firing rates and number of peaks (Fig. 3, Supplementary Video 1, and Supplementary Fig. 7). At week 8, oscillation patterns shifted to lower frequency with fewer calcium peaks, lower burst firing rates, higher amplitudes, longer burst durations, and larger decay times (Fig. 3, Supplementary Video 2, and Supplementary Fig. 7). The plateau shape of the oscillations at week 8 indicated multiple neuronal action potentials contributing to the calcium oscillation (Fig. 3a). The decrease in the number of peaks from weeks 6 to 8 suggested more synchronous calcium transients, indicating a more densely connected mature network. From weeks 10 to 14, burst duration, decay time, and number of peaks did not change significantly, but amplitude and percentage of active organoids decreased, suggesting a plateau in differentiation around week 8.
In addition to whole organoid analysis, ∆F/F was quantified in single neurons for at weeks 4–10 (Supplementary Fig. 7). Maximum intensity z-projections of time course videos showed that neuronal networks at weeks 4 and 6 were less connected compared to weeks 8 and 10 (Supplementary Fig. 7). At weeks 4 and 6, neurons spiked at higher frequencies and with less synchronization (Supplementary Fig. 7a and b). By weeks 8 and 10, larger burst amplitudes and longer burst durations are likely contributed to multiple action potentials across different neurons, which were spiking simultaneously (Supplementary Fig. 7c and d). At week 10, the propagation of an action potential across connected neurons was observed by the slightly delayed peak burst amplitude of region of interest (ROI) 1 compared to ROIs 2 and 3 (Supplementary Fig. 7d).
To measure network activity over time, HD-MEAs were used to obtain high spatial and temporal resolution of organoids’ electrical activity across two different time periods (weeks 6-to-9 and 10-to-13) (Fig. 4). Representative raster plots indicated differences in spontaneous electrical activity in organoids depending on their age (Fig. 4a, b). The 6-to-9-week organoids have a significantly higher burst frequency, number of spikes within burst, and percent active area than those in the later time point group (Fig. 4c). They also had significantly shorter interburst intervals compared to the more mature group, consistent with the calcium imaging data in Fig. 3.
To further assess the organoids’ functionality, neuronal connectivity and criticality were quantified from the same HD-MEA time course data (Figs. 5 and 6). In both age groups, changes in functional connections between electrodes were observed over time on the HD-MEA (Fig. 5a). Denser connections and more active electrodes were observed in the 10-to-13-week group compared to the 6-to-9-week group, as denoted by the thickness of the black lines and red electrodes, respectively in the connectivity graphs (Fig. 5a). However, while both groups showed significant increases in the number of nodes over time, the 10-to-13-week group had a significantly lower number of nodes overall in their functional connectivity matrices compared to the 6-to-9-week group (Fig. 5b). To quantify the shift in the strength of the edges over time, an edge weight distribution was calculated by measuring the fraction of total possible edges that are realized (Fig. 5c). Interestingly, most edges were activated across all samples over time (Fig. 5c). The 10-to-13-week group showed no significant changes over time, while the 6-to-9-week group showed a temporary significant decrease in strength of edges at day on MEA (DOM) 7 (Fig. 5c). Finally, the organoid’s modularity was significantly different across age groups and significantly decreased in both age groups over time, indicating that the networks started with multiple communities but then became more of a single community over time (Fig. 5d). The decrease in modularity may also be due to an increased number of nodes. Despite the similarity in modularity, the 10-to-13-week group maintained a significantly higher modularity over time, indicating that it maintained more communities or network connections (Fig. 5d).
Criticality is a state in which complex systems such as a brain operates at the critical point between organization and randomness, demonstrating how neuronal networks may navigate between the two stages of chaos and order43. The critical point state is key for brain functionality, as it operates at its optimal and most efficient computational capacity and is highly sensitive to external stimuli during this stage. Organoids exhibited properties of criticality over the course of differentiation (Fig. 5e–g). The more mature 10-to-13-week group showed a consistently lower and more tightly regulated Deviation from Criticality Coefficient (DCC) value and higher Branching Ratio (BR), approaching 1, compared to the 6-to-9-week group (Fig. 5e). While the BR in the 10-to-13-week group decreased non significantly over the period of 3 weeks on the HD-MEAs, the 6-to-9-week group gradually increased significantly, demonstrating maturation and pursuit of criticality and thus a stable state (Fig. 5f). Additionally, the Shape Collapse error (SCe) for the 10-to-13-week group was significantly lower than that of the 6-to-9-week group, indicating a more accurate scaling of avalanches of varying durations to an universal shape in the 10-to-13-week group (Fig. 5g). This analysis suggests that the 10-to-13-week group was in a more critical state compared to the 6-to-9-week group. However, over time, both the BR and SCe appeared to converge for both groups, suggesting that the 6-to-9-week group exhibited increasingly critical dynamics, while the 10-to-13-week group showed diminishing critical dynamics on the MEA over time.
Pharmacological characterization of synaptic transmission changes neuronal bursting activity and immediate early gene expression
To validate reactiveness to network modulations, pharmacological agents were used to cause neuronal depolarization and disrupt excitatory glutamatergic synaptic transmission. Expression of IEGs and synaptic plasticity-related genes was measured 2 hours after exposure to pharmacological agents and compared to the corresponding untreated control in two age groups (8 weeks and 13 weeks) (Fig. 6). To disrupt excitatory glutamatergic synaptic transmission, organoids were treated with 2,3-dioxo-6-nitro-7-sulfamoyl-benzo[f]quinoxaline (NBQX), an AMPA receptor antagonist, D-2-amino-5-phosphonovalerate (AP5), a NMDA receptor antagonist, 4-Aminopyridine (4-AP), a voltage-gated potassium (Kv) channel antagonist, and bicuculline, a GABA receptor antagonist, were used to enhance neuronal depolarization and synaptic transmission (Fig. 6).
Bicuculline induced a slight increasing trend in gene expression across both age groups (Fig. 6a), while exposure to 4-AP led to significant changes in NPAS4 and FOS expression at both age groups. Expression of ERG1 was significantly induced only at week 13. Lastly, ARC expression showed an increased trend in expression after 4-AP exposure (Fig. 6a). No significant changes in gene expression were seen after exposure to NBQX and AP5 individually or combined (Fig. 6 and Supplementary Fig. 8a–d). Since IEGs were more strongly perturbed at week 13, the effects of these chemicals on electrophysiological activity were assessed in this age group. Organoids were exposed to the pharmacological agents directly on the HD-MEA at DOM 29. Network recordings were taken before the addition of the chemicals as a baseline. Network activity was then recorded immediately after exposure, followed by 2 and 4 hours thereafter, after which the recorded parameters were compared to baseline activity (Fig. 6c and d). 4-AP and bicuculline increased network activity while NBQX + AP5 decreased network activity over time (Fig. 6c). More specifically, bicuculline caused an insignificant increasing trend in mean burst frequency and interburst interval coefficient of variation (CV) over time, a significant increase in percent of spikes within bursts 0 minutes after and an increasing trend in percent of spikes within bursts 2 and 4 hours after exposure. In addition, bicuculline caused no significant changes or trends in burst duration over time. 4-AP exposure caused a significant increase in mean burst frequency and an increasing trend in mean percent of spikes within bursts 0 minutes after. In addition, the percentage of spikes within bursts maintained an increasing trend within 2 and 4 hours after exposure. 4-AP also caused a decreasing trend in burst duration that was maintained over time. Finally, 4-AP caused no significant changes or trends in interburst interval CV over time. Additionally, NBQX + AP5 exposure completely abolished network bursting activity (Figs. 6c and d, Supplementary Fig. 8). Overall, NBQX + AP5 significantly decreased mean burst frequency, interburst interval CV, burst duration, and percentage of spikes within bursts from 0 minutes to 4 hours. Interestingly, we found that NMDA receptors are largely responsible for neuronal network bursting, as exposure to only AP5 was enough to abolish the bursting, while blocking only AMPA receptors with NBQX only partially reduced the bursting (Supplementary Fig. 8). These results agreed with previous reports showing that ketamine and xenon, which act on NMDA receptors, lead to burst silencing and reduction in vitro44,45. No changes in firing rate, spikes per burst, and burst duration were seen after NBQX application alone, but when AP5 or NBQX + AP5 was applied, no bursts were observed; therefore, firing rate, spikes per burst, and burst duration were not quantifiable (Supplementary Fig. 8).
While NBQX + AP5 reduced the network bursting, some spiking activity was still seen (Fig. 6c and Supplementary Fig. 8e). To elucidate which type of neurons contributed to the remaining activity, different modulators of the glutamatergic, GABAergic, dopaminergic, noradrenergic and cholinergic neurons were tested in different combinations (Supplementary Fig. 9). The first set of experiments further blocked NMDA and AMPA receptors, by increasing concentrations of NBQX + AP5 to 40 µM. As for the previous 20 µM exposure (Fig. 6c and Supplementary Fig. 8e), all bursting was abolished and spiking was diminished, but increased concentration did not remove the activity completely. Therefore, next, we blocked different types of neurons. First, 10 µM bicuculline and 10 µM CPG 55845 hydrochloride (CPG 55845) were added to block GABA(A) and GABA(B) receptors in addition to NMDA and AMPA receptors. Upon addition of bicuculline and CPG 55845, a further reduction in spiking activity was observed (Supplementary Fig. 9a, last panel, 10–20 min). The quantification of the effect after blocking excitatory receptors and the subsequent blocking of inhibitory neurons is shown in Supplementary Fig. 9d.
Next, inhibitory and excitatory receptors were blocked simultaneously by adding 40 µM NBQX, 40 µM AP5, 10 µM bicuculline, and 10 µM CPG 55845 at the same time with similar effects as before (Supplementary Fig. 9b). Subsequently, 5 µM 6-hydroxydopamine (6-OHDA) and 10 µM haloperidol were added, which induce lesions in dopaminergic (and noradrenergic) neurons and block dopamine D2 receptors, respectively. This resulted in an even more pronounced reduction in spiking activity (Supplementary Fig. 9b, the last panel). Supplementary Fig. 9e quantifies the corresponding changes in bursting and spiking metrics. These results suggest that most of the remaining activity originated from dopaminergic and noradrenergic neurons.
The last population of neurons examined for chemical modulation effects was cholinergic neurons. To achieve this, 40 µM NBQX, 40 µM AP5, 10 µM bicuculline, 10 µM CPG 55845, 5 µM 6-OHDA, and 10 µM haloperidol were added simultaneously to block excitatory, inhibitory, dopaminergic and noradrenergic neurons. The same effect as shown in Supplementary Fig. 9b was also observed in Supplementary Fig. 9c, demonstrating a drastic reduction in bursting and spiking, as quantified in Supplementary Fig. 9f. Finally, 1 nM chlorpyrifos (CPF)-Oxone was added to that same well to attempt to increase activity in the network by activating cholinergic neurons. However, no changes in bursting or spiking were observed, as shown in Supplementary Fig. 9c and quantified in Supplementary Fig. 9f. Therefore, findings suggest that the remaining activity likely originated from partially blocked AMPA and NMDA receptors, as well dopaminergic and noradrenergic neurons.
Theta-burst stimulation modulated synaptic plasticity
To generate input-specific evoked activity from electrical stimulation, Theta Burst Stimulation (TBS) was delivered to 14-week-old organoids 4 times with 13-minute intervals between TBS (Fig. 7a) on the HD-MEA. Two independent experiments (referred to as A and B in the figures and below) were conducted. In both experiments, four to five organoids were seeded on each well at week 9.5 of differentiation and grown on the MEA until 33 ± 1 DOM before stimulation (Supplementary Fig. 10a and Supplementary Fig. 11a). First, we recorded the baseline network activity across all wells and found that basal activity was lower in wells 1A-3A than wells 4A-6A and 1B-4B (Supplementary Figs. 10b and 11b). The MaxWell HD-MEA has an electrode size of 8.75 × 12.50 µm², and the electrode center-to-center distance is 17.5 µm, allowing one neuron to be recorded by multiple electrodes. For input-specific synaptic plasticity, one neuron from each well was identified based on its footprint (Supplementary Figs. 12a and 13a) and spike-sorted neuron traces using the Axon Tracking assay in the MaxLab Live Software (Supplementary Figs. 12b and 13b). Then, 32 electrodes focusing on a single neuron in each well were stimulated using a modified version of previously described LTP induction protocols46,47,48 (Fig. 7a). To optimize the stimulation of each neuron, electrodes along the entire neuron including the soma and axon were targeted for stimulation.
To investigate short-term changes in evoked activity, total evoked activity per bin (10 ms), total spikes, and total active area were measured. Active area before and after each stimulation are shown for all wells in both experiments (Fig. 7b and Supplementary Fig. 14a). Wells 4A-6A, 1B-4B showed substantial changes in active area in response to the stimulus while wells 1A-3A showed little changes (Fig. 7b, Supplementary Fig. 14a).
Representative evoked activity heatmaps from wells 4A-6A and 1B-4B demonstrated strong short-term responses within milliseconds following stimulation (Fig. 7c, Supplementary Figs. 15, and 16). Wells with lower baseline activity (e.g., 1A-3A) did not show any response.
To determine an activity threshold, we treated one well with NBQX/AP5 to block glutamatergic receptors-dependent synaptic plasticity. The 90th, mean, and 10th percentile responses from the NBQX/AP5-treated well is shown overlayed on the plots as the dotted blue, black, and red lines, respectively (Fig. 7d, Supplementary Fig. 14c). Wells 1A-3A did not exceed this threshold, while wells 4A-6A and 1B-4B consistently did across all four TBS sets (Fig. 7d, Supplementary Fig. 14c).
Aggregated data for active area, total spikes, and evoked activity showed that wells 4A-6A and 1B-4B had a distribution skewed to the right of 0, while wells 1A-3A only exhibited a mode around 0 (Fig. 7e, Supplementary Fig. 14d). Wells 1A-3A, with lower baseline activity and connectivity compared to wells 4A-6A and 1B-4B, did not respond above threshold, whereas the shift to the right of 0 in wells 4A-6A and 1B-4B suggests short-term potentiation, as stimulation led to short-term increases in activity.
To further confirm STP, connectivity and criticality was quantified across each TBS for well that exhibited STP (Well 4A, 5A, 6A, 1B, 2B, 3B, 4B) across both experiments (Fig. 8a–c). We observed a significant increase in the number of nodes and edges while modularity significantly decreased after every stimulation indicating a more connected network of neurons after stimulation (Fig. 8a and b). These results suggest that wells that exhibit evoked electrical activity also became more connected, more of a single community. While metrics of connectivity were extremely consistent across stimulations, metrics of criticality varied across the stimulations including the deviation of criticality coefficient and shape collapse error (Fig. 8c). There were no significant differences in either metric (Fig. 8c). The branching ratio, however, showed a significant increase following stimulation (Fig. 8c), suggesting that the neural organoids are becoming more critical after stimulation.
Long-term effects of TBS on organoids were assessed by quantifying criticality, connectivity, network-level dynamics (including interspike interval, interspike interval CV, and firing rate), and spike sorted unit (aka neuronal) level changes in firing rate in wells that demonstrated input-specific STP (Supplementary Figs. 17 and 18). These metrics were quantified before and at 60 – 180 minutes post-TBS. There were no significant differences before and after stimulation (after stimulation includes grouped data from 60-180 minutes post stimulation) for all connectivity and criticality metrics (Supplementary Figs. 17 and 18). In addition, there were no consistent trends over time for network level interspike interval, interspike interval CV, or firing rate for wells exhibiting STP (Supplementary Figs. 17 and 18). Interspike Interval (ISI) was calculated with a 4 Hz threshold (up to 250 ms) to account for changes in theta entrainment/phase locking. Well 4A, 5A, 1B, 2B, 3B, and 4B showed no long-term changes in ISI after stimulation (Supplementary Figs. 17c and 18c). While well 6A showed a significant increase in ISI 60–180 min post stimulation, apart from at 90 minutes (Supplementary Fig. 17c). The CV was used to measure ISI variability across timepoints7. A CV of 2.5 indicates a perfect Poisson process7,49, while a CV near zero indicates a perfectly periodic spike train. All wells analyzed showed no significant change in ISI CV after stimulation (Supplementary Figs. 17c and 18c). In addition, none of the wells showed any changes in firing rate over time compared to the baseline (Supplementary Figs. 17d and 18d). Since network level ISI, ISI CV, and firing rate graphs were analyzed based on overall trends across grouped units for each well and not on a unit level, the results indicate that there are no network level changes in interspike interval or firing rate. Together, these results suggest that there are no long-term changes to the network dynamics after TBS to an individual neuron. This is expected as the TBS was input-specific; therefore, we focused on unit-level changes rather than network-level changes. Upon further investigation of unit-level firing rates over time, linear regression analysis showed that all wells had units that were either potentiated or depressed and were maintained for 180 minutes post-stimulation (Fig. 8d). To account for random changes in firing rate, the same analysis was performed on a time series recording without stimulation. In this recording, only one out of 40 units showed potentiated firing rate over time, therefore, this can be considered the noise level of changes due to basal synaptic plasticity within the organoids. Wells 4A-6A and 2B-4B show potentiated and depressed units above this level (Fig. 8d). Representative LTP and LTD units are shown in Fig. 8d. Together, these data indicate there are input-specific TBS-induced changes in connected neurons (aka units) over hours but not the overall network, supporting the use of this model to modulate input-specific short- and long-term synaptic plasticity and detect changes in synaptic plasticity in connected neurons. In addition, due to the nature of input-specific synaptic plasticity on a HD-MEA, it is difficult to control what kind of neuronal circuits are being modulated, whether they are increased or decreased, which could explain why the population of LTP compared to LTD is different in each well.
