Human neural organoid microphysiological systems show the building blocks necessary for basic learning and memory

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.

Fig. 1: Schematic overview of the experimental approach.

a Experimental timeline. Created in BioRender. Alam El Din, D. (2025) https://BioRender.com/v4k2lpzb Overview of avalanche and network connectivity analysis for time series electrophysiology data from organoids plated on HD-MEAs. Created in BioRender. Alam El Din, D. (2025) https://BioRender.com/trj7ehfc Schematic representation of synaptic transmission modulation by pharmacological and electrical stimuli to induce synaptic plasticity. Adapted from Kim, S. (2025). Long-Term Potentiation. https://app.biorender.com/biorender-templates/details/t-61006a6924e0d000a40de3a1-long-term-potentiation. Created in BioRender. Alam El Din, D. (2025) https://BioRender.com/eqpwhdw.

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.

Fig. 2: Expression of glutamatergic and GABAergic receptor and synaptic plasticity-related genes in neural organoids over course of differentiation.
figure 2

a Representative immunocytochemistry images of organoids showing postsynaptic marker (HOMER1) and presynaptic marker (SYP) in 8- and 12-week cultures. In composite images, HOMER1 is shown in blue, and SYP is shown in yellow. Scale bars are 100 µm, 50 µm, and 10 µm, respectively. b Presence of inhibitory post-synaptic marker (Gephyrin), pre-synaptic marker (SYN1) and dendrites (MAP2) in 8- and 12-week organoids. In composite images, Gephyrin is shown in blue, SYN1 in yellow, and MAP2 in grey. Scale bars are 100 µm and 50 µm, respectively. For a, and b, all images were taken at 20x, 100x, and 100x + 4x zoom and processed with ImageJ for visualization. c Gene expression of Gamma-Aminobutyric Acid Type A Receptor Subunit Alpha1 (GABRA1), Glutamate Ionotropic Receptor NMDA Type Subunit 1 (GRIN1), Glutamate [NMDA] Receptor Subunit Epsilon-1 (GRIN2A), and Glutamate [NMDA] Receptor Subunit Epsilon-2 (GRIN2B), Glutamate Ionotropic Receptor AMPA Type Subunit 1 (GRIA1), homer scaffold protein 1 (HOMER1) in organoids over the course of differentiation. d Representative immunocytochemistry images of weeks 8 and 12 organoids stained for Neuronal Pentraxin 2 (NPTX2), Activity-Regulated Cytoskeleton-associated protein (ARC), cAMP response element-binding protein (CREB), and Brain-Derived Neurotrophic Factor (BDNF). Scale bar is 100 µm. e Gene expression over the course of differentiation of immediate early genes (IEGs) ARC, BDNF, Neuronal PAS Domain Protein 4 (NPAS4), NPTX2, Fos proto-oncogene AP-1 transcription factor subunit (FOS), and Early Growth Response 1 (EGR1). f Gene expression of synaptic slasticity- related genes: CREB, calcium/calmodulin-dependent protein kinase II A (CAMK2A), Synaptic Ras GTPase-activating protein 1 (SYNGAP1). g Gene expression of synaptic plasticity -related miRNAs. For all gene expression plots, data is shown as a box and whisker plot (with the box extending from the 25th to 75th percentiles) and represented as log2(Fold Change) normalized to NPCs from 2-3 independent experiments with 3 technical replicates each. In all qPCR experiments, ACTB was used as a reference gene.

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.

Fig. 3: Neural organoid calcium oscillatory dynamics across different time points to show maturation of spontaneous network bursting.
figure 3

a Representative changes in fluorescence over resting fluorescence (∆F/F) graphs across 360 seconds for each time point from week 2 (W2) to week 14 (W14) of differentiation. b Average rise time, peak amplitude, firing rate, decay time, burst duration, number of peaks, and percentage of active organoids shown across different time points. At least 8 individual organoids across at least 3 independent experiments were imaged and quantified for each time point. Data is shown as box and whisker plots (with the box extending from the 25th to 75th percentiles). Statistics were performed using one-way ANOVA and a Tukey post-hoc test. Changes over time were significant for rise time (p < 0.05), burst firing rate (p < 0.0001), peak amplitude (p < 0.0001), decay time (p < 0.01), burst duration (p < 0.001), and total number of peaks per organoid (p < 0.0001). Pairwise comparisons are shown on the figure: # = Significant difference from week 4, Ŧ = Significant difference from week 6, $ = Significant difference from week 8, ¥ = Significant difference from week 10, ȼ = Significant difference from week 12, • = Significant difference from week 14, * = Significant difference from all weeks. For exact p values see Supplementary Tables 4–9. See also Supplementary Fig. 7 for single neuron calcium imaging analysis.

Fig. 4: Changes in spontaneous electrical activity in neural organoids throughout development.
figure 4

Representative raster plots and active area plots from HD-MEA recordings showing spontaneous electrical activity over time during a weeks 6-to-9 and b weeks 10-to-13 of differentiation. DOM: Days on MEAs. c Network dynamic metrics from both organoid age groups over time (blue line represents 6-to-9 week organoids, red line – 10-to-13 week organoids. The line shown represents mean and the shaded region represents the standard deviation plotted from 2 independent experiments with 5 to 6 HD-MEA wells per group per experiment with 2–5 organoids per well (n = 11–12 wells per age group). Statistics were performed using a mixed-effects model with matching and a Tukey post-hoc test. p < 0.05 was considered significant. For exact p values from pairwise comparisons, see the Supplementary Data 3 file. ISI: Interspike Interval. IBI: Interburst Interval.

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).

Fig. 5: Neural organoids show highly interconnected neuronal networks and criticality throughout development.
figure 5

a Representative plots of functional connectivity at day on MEA (DOM) 3, 9, 15, and 21 for the week 6-to-9 and week 10-to-13 old organoids. For clarity of visualization, only the 200 connections (edges) with the highest mutual information are shown. Each red dot represents an electrode, and the lines indicate the connections between electrodes. The thickness of the line indicates the weight of connectivity. b Average number of nodes; c Average fraction of total possible edges; d Average modularity over time in week 6-to-9 and week 10-to-13 organoids. e Deviation from Criticality Coefficient (DCC). f Branching Ratio (BR) g Shape collapse error (SCe) over time in 6-to-9 week and 10-to-13 week old organoids. In bd the line shown represents the mean and the shaded region represents the standard deviation plotted from 2 independent experiments with 5 to 6 HD-MEA wells per group per experiment with 2–5 organoids per well (n = 11–12 wells per age group). Panels e–g show regression lines (blue line- 6–9 week old organoids, red line – 10-13 week old organoids) with a 95% confidence interval. Data plotted is from 2 independent experiments with 5-6 HD-MEA wells per group per experiment (n = 11–12 wells per age group). Statistics were performed using a two-way ANOVA and a Tukey post-hoc test. p < 0.05 was considered significant. For exact p values from pairwise comparisons, see the Supplementary Data 3 file.

Fig. 6: Pharmacological characterization of synaptic transmission changes of neuronal spiking and bursting activity and Immediate Early Gene expression.
figure 6

Expression of ARC, NPAS4, FOS, and EGR1 after 2 hours of exposure to 20 µM AP5 + 20 µM NBQX (pink box), 10 µM bicuculline (green box) and 100 µM 4-AP (purple box) in a 8-week and b 13-week-old organoids, represented as box and whisker plots (25th to 75th percentiles) and as log2(Fold Change) normalized to negative control (organoids with no chemical treatment = 2 h control). ACTB was used as a reference gene. The data represents 3 independent experiments with 2 technical replicates each for 8 weeks and 4-5 independent experiments with 2 technical replicates each for the 13-week time point. Statistics were calculated based on the replicate average from each independent experiment, with one-way ANOVA and post-hoc Dunnett’s tests *p < 0.05, ***p < 0.001, ****p < 0.0001 c Representative raster plots from MEA recordings in 13-week-old organoids (from 6 wells per condition) before and after treatment with bicuculline, 4-AP, and NBQX + AP5. d Burst frequency, Interburst interval coefficient of variation, burst duration, and percentage of spikes within bursts plotted as box and whisker plot (with the box extending from the 25th to 75th percentiles) for bicuculline, 4-AP, and NBXQ + AP5 treated wells prior to (baseline), 0 mins, 2 hours, and 4 hours after exposure. The data represents 3 independent experiments with 2 HD-MEA wells per experiment per chemical (n = 6). Statistical significance was calculated with repeated measures ANOVA with post-hoc Dunnett tests. p < 0.05 was considered significant. Pairwise comparisons can be seen in the Supplementary Tables 10-21 and significant groups are shown in the figure. ARC – Activity-Regulated Cytoskeleton-Associated Protein; NPAS4 – Neuronal PAS Domain Protein 4; FOS – Fos proto-oncogene AP-1 transcription factor subunit; EGR1 – Early Growth Response Protein 1; AP5 – 2-Amino-5-phosphonopentanoic acid (an NMDA receptor antagonist); NBQX – 2,3-dihydroxy-6-nitro-7-sulfamoyl-benzo[f]quinoxaline (an AMPA receptor antagonist).

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.

Fig. 7: Theta-burst stimulation modulated short-term plasticity for Experiment A.
figure 7

a Graphical summary of TBS protocol. i-The TBS was performed four times spaced by 13 minutes. ii-Within each TBS there are 10 trials with four spikes per trial. iii-The schematic of each trial. b Percent active area before and after stimulation across all 6 wells. Wells 4A–6A show consistent increase or decrease in active area in response to stimulation while wells 1A–3A show little change. c Representative heat map evoked activity response for wells 4A–6A. Bin size is equal to 10 ms. The stimulation pulses are the light grey vertical lines, and the dashed orange lines indicate the start/stop time of the analysis window for calculating evoked activity. d percentage of active electrodes, total spikes, and evoked activity for wells 1A-3A and then 4A-6A. Purple circle represents well 1A, pink square – well 2 A, turquoise triangle – well 3A, green circle – well 4A, blue square – well 5A, and yellow triangle – well 6A. The data represents the mean with 10th to 90th percentile for each well. The 90th percentile response of a well treated with NBQX/AP5 before and during stimulation is shown with a blue dashed line overlayed on all graphs. The mean response of a well-treated with NBQX/AP5 before and during stimulation is shown in a black dashed line overlayed on all graphs. The 10th percentile response of a well treated with NBQX/AP5 before and during stimulation is shown in a red dashed line overlayed on all graphs. Responses above this NBQX/AP5 region indicate responses generated by glutamatergic receptors. e Histograms of total evoked activity per bin (bin size of 10 ms), total spikes, and total active area. The top three graphs show data aggregated across all electrodes for all 4 TBS for wells 1A–3A, and the bottom three graphs show data aggregated across all electrodes for all 4 TBS for wells 4A-6A. Wells 1A-3A show little to no response while wells 4A-6A indicate evoked responses on the millisecond timescale.

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.

Fig. 8: Theta-burst stimulation drives short-term changes in connectivity and criticality and long-term potentiation and depression of neuronal units.
figure 8

a Connectivity metrics for all wells that demonstrated STP. b Representative connectivity graph before and immediately following TBS #1 for well 1B. c Criticality metrics for all wells that demonstrated STP. The data represents the mean of seven biological replicates from two independent experiments. A Wilcoxon matched-pairs signed rank test was performed to determine statistical significance for a and c, **p < 0.01, ***p < 0.001, ****p < 0.0001. Exact p-values are listed in Supplementary Tables 22–27. For a and c green circle represents well 4A, blue square – well 5 A, yellow triangle -well 6A, red line – well 1B, purple circle – 2B, open blue circle – 3B, open pink square – 4B. d) Quantification of input-specific long-term potentiation (LTP) and depression (LTD) by measuring firing rate over time in neuronal units. Two example units demonstrating either LTP (blue line with 95% confidence intervals depicted with dashed lines) or LTD (red line, with 95% confidence intervals depicted with dashed lines) are shown on the left. The proportion of neuronal units that demonstrated LTP (red) or LTD (blue) across wells is demonstrated and quantified on the right.

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.

Continue Reading