Deepening the role of excitation/inhibition balance in human iPSCs-derived neuronal networks coupled to MEAs during long-term development

We used two previously characterized hiPSCs lines genetically modified to generate homogeneous populations of excitatory and inhibitory neurons thanks to the forced expression of the transcription factors neurogenin-2 (Ngn2) and Achaete-scute homolog 1 (Ascl1) (Mossink et al 2022). We received the hiPSCs lines in frozen vials from Dr Frega (University of Twente). Both lines were generated from fibroblasts. Control line 1 (C1, healthy 30 years-old female) was reprogrammed via episomal reprogramming (Coriell Institute for medical research, GM25256). Control line 2 (C2, healthy 51 years-old male) was reprogrammed via a non-integrating Sendai virus (KULSTEM iPSC core facility Leuven, Belgium, KSF-16-025). The thawed hiPSCs cultures were resuspended into E8Flex medium (Thermo Fisher Scientific) supplemented with E8 supplements (2%, Thermo Fisher Scientific), penicillin/streptomycin (1%, Sigma-Aldrich), G418 (50 μg ml⁻¹, Sigma-Aldrich) and puromycin (0.5 μg ml⁻¹, Sigma-Aldrich). Glutamatergic cortical Layer 2/3 neurons were derived from C1 by overexpressing mouse neuronal determinant Ngn2 upon doxycycline (4 μg ml⁻¹, Sigma-Aldrich) treatment (Frega et al 2017). GABAergic neurons were derived from C2 by overexpressing mouse neuronal determinant Ascl1 (Addgene 97 329) upon doxycycline treatment and supplementation with forskolin (10 μM, Sigma Aldrich) (Mossink et al 2022). The hiPSCs can be considered differentiated into neurons after 3 weeks of doxycycline and forskolin treatment (Shi et al 2016, Frega et al 2017). The cultures were maintained into the incubator at stable condition (37 °C, 5.5% CO₂, 95% humidity atmosphere). The medium was refreshed at 50% every two days.

The neuronal cultures were treated with the adapted pipeline presented in (Wang et al 2023). Briefly, PDMS rings were used to favor the attachment and to both confine and protect the cells without creating patterned networks and without affecting the cells' viability (Liu et al 2019, Ming et al 2021). Each ring had a diameter of 6 mm and height of about 1 mm. The PDMS rings were obtained using a mixture of PDMS base (Syligard 184) and curing agent at a 10:1 (w/w) ratio, which was polymerized in oven (80 °C for 15 min). The day before the neuronal cultures plating, the rings were positioned and let adhere to the MEA surface. The devices were sterilized in an oven (120 °C for 3 h) and pre-coated over night with poly-L-Ornithine (50 μg ml⁻¹, Sigma-Aldrich) and human laminin (20 μg ml⁻¹, BioLamina) (Hartmann et al 2023). At DIV 0, neurons (iNs) were co-plated with rat astrocytes (A) to favor neuronal growth and maturation (Tang et al 2013) in a ratio equal to 70:30 (iNs:A) (figure 1(c)). The glutamatergic (E) and GABAergic (I) neurons were mixed in different proportions to obtain 5 configurations: E:I 100:0, 75:25, 50:50, 25:75, 0:100 (figure 1(b)). The neuronal cultures density was equal to 1200 cells mm⁻². During the first week, the neuronal cultures were maintained in Neurobasal medium (Thermo Fisher Scientific) supplemented with B27 supplements (2%, Thermo Fisher Scientific), penicillin/streptomycin (1%, Sigma-Aldrich), stable L-Glutamine (1% GlutaMax 100x, Thermo Fisher Scientific), human Brain-Derived Neurotrophic Factor (BDNF, 10 ng ml⁻¹, Sigma-Aldrich), human Neurotrophin-3 (NT-3, 10 ng ml⁻¹, Sigma-Aldrich), doxycycline (4 μg ml⁻¹, Sigma-Aldrich) and forskolin (4 μg ml⁻¹, Sigma-Aldrich). After 7 DIV, fetal bovine serum (FBS, 2%, Thermo Fisher Scientific) was added to the aforementioned supplemented medium to support astrocytes. After 14 DIV, doxycycline and forskolin were removed from the medium. The neuronal cultures on MEAs were stably maintained in the incubator (same conditions specified in the previous paragraph) up to 98 DIV. The medium was refreshed at 50% each two days.

Zoom In Zoom Out Reset image size

The neuronal cultures were plated on commercial MEAs, i.e. 60 TiN electrodes arranged in an 8 × 8 grid layout on a glass surface (Multi Channel Systems—MCS, Reutlingen, Germany). The electrodes were 30 μm in diameter and equally spaced with a 200 μm of distance. The electrophysiological data were acquired with the MEA2100 recording system (MCS) with the MC_rack software (MCS). After a 10 min period of acclimatation out of the incubator, the electrophysiological spontaneous activity was recorder for 15 min in stable condition (37 °C, 5% CO₂) and sampled at 10 kHz. The neuronal networks activity was recorded from DIV 35 until DIV 98. The number of devices recorded for each configuration in each temporal step is reported in table 1. All collected data were obtained from neuronal networks derived from three different experimental batches. It is worth noticing that we observed the entire development of the 56%, 86%, 63%, 63% and 83% of the cultures for the 0:100, 25:75, 50:50, 75:25, and 100:0, respectively.

In the following subsections, we explain in detail the analyses carried out on the data and, in particular, the parameters we extrapolated to characterize the networks. The parameters are summarized in the glossary (supplementary).

2.4.1. Spiking and bursting activity

2.4.2. Network bursting activity

From the network burst detection, the following parameters were derived: mean network bursting rate (NBR) computed for each culture individually as the number of network bursts detected in a minute, and the network burst duration (NBD). Then, we computed the cumulative spike time histograms (STHs): the cumulative instantaneous firing rate (bin = 10 ms) related to all detected network bursts were aligned with respect to the longest network burst peak. From the cumulative STH, we inferred the peak amplitude by considering the maximum value of the shape (STH_peak). We obtained the average cumulative STH for each configuration by averaging the cumulative STHs of each culture belonging to the considered configuration. We computed the 90% confidence interval by considering all the shapes of the same configuration to evaluate the robustness of the dataset. We evaluated the evolution of the cumulative STH by analyzing the trend of the rising and decaying phases. The rising phase of the cumulative STH was considered as the part of the curve ranging from the time required to pass from the 10% to the 90% of the peak amplitude. The rising phase has been fitted with a double exponential function:

$$\begin{align}{\boldsymbol{r}}\left( {\boldsymbol{t}} \right) = {{\boldsymbol{a}}_1}{{\mathbf{e}}^{\left( {{{\boldsymbol{b}}_1} - 1} \right){\boldsymbol{t}}}} + {{\boldsymbol{a}}_2}{{\mathbf{e}}^{\left( {{{\boldsymbol{b}}_2} - 1} \right){\boldsymbol{t}}}} + {\boldsymbol{c}}\end{align} \tag{ 1 }$$

where ${a_1}$, ${a_2}$ and c are the coefficients of the fitting equation. The smaller exponent between ${b_1}$ and ${b_2}$ represents the neuron recruitment rate (τ_rise) (Van De Vijver et al 2019). The decay phase of the cumulative STH was considered as the part of the curve ranging from the time required to pass from 90% to 10% of the peak amplitude of the curve after the peak itself. The decay phase has been fitted with a first-degree exponential function:

$$\begin{align}{\boldsymbol{{d}}}\left( {\boldsymbol{{t}}} \right) = {{\boldsymbol{{a}}}_1}{{\text{e}}^{{\boldsymbol{{\delta t}}}}} + {{\boldsymbol{{a}}}_2}\end{align} \tag{ 2 }$$

where ${a_1}$ and ${a_2}$ are the coefficients of the fitting equation and δ represents the network modulation (τ_decay) during the network bursts events (Marom and Shahaf 2002). All the aforementioned network bursting statistics (i.e. NBR, NBD, STH_peak, τ_rise, and τ_decay) related to the different configurations (i.e. E:I 100:0, 75:25, 50:50, 25:75, 0:100) at each DIV were obtained by averaging the parameters of each culture belonging to the considered configuration.

2.4.3. Functional connectivity

The statistical interdependence between neuron pairs was evaluated by extracting the functional connectivity. In order to distinguish both excitatory and inhibitory functional connections, we inferred connectivity by adopting the Total Spiking Probability Edges (TSPE) algorithm (De Blasi et al 2019). Briefly, once the cross-correlograms were obtained, an edge filter was applied to them. Such procedure allowed identifying the local maxima and minima, corresponding to excitatory and inhibitory connections, respectively. After the connectivity matrix has been computed, we used a threshold to discard those values that did not reasonably identify actual functional connections or that were not compatible with a physiological propagation speed. First, we applied a spatial filter to keep only the functional connections that guarantee a propagation speed between 30 and 300 mm s⁻¹ (Jacobi and Moses 2007). Then, we used a hard threshold independently evaluated for the excitatory and inhibitory connections (Pastore et al 2018b), defined as follows:

$$\begin{align}{\boldsymbol{{T}}}{{\boldsymbol{{H}}}_{{\boldsymbol{{e}}},{\boldsymbol{{i}}}}} = \,{{\boldsymbol{{\mu }}}_{{\boldsymbol{{e}}},{\boldsymbol{{i}}}}} \pm {{\boldsymbol{{n}}}_{{\boldsymbol{{e}}},{\boldsymbol{{i}}}}}\,{{\boldsymbol{{\sigma }}}_{{\boldsymbol{{e}}},{\boldsymbol{{i}}}}}\end{align} \tag{ 3 }$$

where μ and σ are the mean and the standard deviation of all non-zero cross-correlation matrix elements, and n is an integer. The subscripts e and i are relative to the excitatory and inhibitory links, respectively. We set for all experiments ${n_e}$ = 1 and ${n_i}$ = 2 (Brofiga et al 2022).

2.4.4. Synchronization

The effect of the different percentages of inhibitory neurons was evaluated in terms of synchronization level of both spiking and bursting activity. For each pair of electrodes k and j, the normalized cross-correlogram $({C_{kj}}\left( \tau \right))$ was computed within a window of 50 ms (T) using a bin of 0.1 ms (${{\Delta }}\tau$). The normalized cross-correlation was obtained using the following equation:

$$\begin{align}{{\boldsymbol{{C}}}_{{\boldsymbol{{kj}}}}}\left( {\boldsymbol{{\tau }}} \right) = \frac{1}{{\sqrt {{{\boldsymbol{{N}}}_{\boldsymbol{{k}}}}{{\boldsymbol{{N}}}_{\boldsymbol{{j}}}}} {\boldsymbol{{\,}}}}}\mathop \sum \limits_{{\boldsymbol{{s}}} = 1}^{{{\boldsymbol{{N}}}_{\boldsymbol{{k}}}}} \mathop \sum \limits_{{{\boldsymbol{{t}}}_{\boldsymbol{{i}}}} = {\boldsymbol{{\tau }}} - \frac{{{\boldsymbol{\Delta }}{\boldsymbol{{\tau }}}}}{2}}^{{\boldsymbol{{\tau }}} + \frac{{{\boldsymbol{\Delta }}{\boldsymbol{{\tau }}}}}{2}} {\boldsymbol{{k}}}\left( {{{\boldsymbol{{t}}}_{\boldsymbol{{s}}}}} \right){\boldsymbol{{j}}}\left( {{{\boldsymbol{{t}}}_{\boldsymbol{{s}}}} - {{\boldsymbol{{t}}}_{\boldsymbol{{i}}}}} \right)\end{align} \tag{ 4 }$$

where ${t_s}$ indicates the timing of an event in the k train.

From the cross-correlation function, we computed the coincidence index ($C{I_0}$), defined as follows:

$$\begin{align}{\boldsymbol{{C}}}{{\boldsymbol{{I}}}_{\boldsymbol{{o}}}} = \frac{{\mathop \sum \nolimits_{{\boldsymbol{{\tau }}} = {\boldsymbol{{\,}}} - \frac{{{\boldsymbol{{i}}}{\boldsymbol{\Delta }}{\boldsymbol{{\tau }}}}}{2}}^{\frac{{{\boldsymbol{{i}}}{\boldsymbol{\Delta }}{\boldsymbol{{\tau }}}}}{2}} {{\boldsymbol{{C}}}_{{\boldsymbol{{kj}}}}}\left( {\boldsymbol{{\tau }}} \right)}}{{\mathop \sum \nolimits_{{\boldsymbol{{\tau }}} = - {\boldsymbol{{T}}}}^{\boldsymbol{{T}}} {{\boldsymbol{{C}}}_{{\boldsymbol{{kj}}}}}\left( {\boldsymbol{{\tau }}} \right)}}\end{align} \tag{ 5 }$$

where the numerator represents the area of the correlogram around zero and i is the number of bins around the zero-bin ($\pm \,$1 ms).

2.4.5. Statistics

At DIV 70, cells were fixed with 4% paraformaldehyde (PFA) for 10 min at room temperature. Subsequently, the cultures were permeabilized with 0.2% Triton-X100 diluted in PBS for 10 min. Non-specific binding sites were blocked with a 1-hour treatment of blocking buffer solution (BBS) composed by PBS, 0.2% Triton-X100, and 5% of fetal bovine serum (FBS, Sigma-Aldrich). The primary antibodies were diluted in BBS and incubated overnight at 4 °C. We used anti-GABA (1:500, Poly Rabbit, A2052, Sigma-Aldrich) and anti-NeuN (1:500, Poly Guinea Pig, Synaptic Systems) as primary antibodies. The day after, we rinsed twice the samples with PBS to remove the excess primary antibodies solution. The secondary antibodies were diluted in BBS and incubated for one hour at room temperature. The following secondary antibodies were used: Alexa Fluor 546 (1:1000, Goat anti-Guinea Pig, ThermoFisher Scientific), and Alexa Fluor 488 (1:700, Goat anti-Rabbit, Invitrogen). We imaged with fluorescence equipped microscope (Olympus BX-51). Images were acquired at 50X magnification. Fluorescent signals were analyzed using Fiji software (Schindelin et al 2012). All the samples were acquired with the same settings, to compare the signals between the different experimental configurations. The number of glutamatergic neurons was defined as the difference between the number of total neurons (NeuN) and the number of GABAergic neurons (GABA A2052). We collected 240 images (80 for each heterogeneous configuration) from different samples.

We morphologically characterized the heterogeneous configurations (i.e. 25:75, 50:50, and 75:25) by performing immunocytochemistry to prove that we can effectively control the different percentages of the plated glutamatergic and GABAergic neurons. Figures 2(a)–(c) shows a representative image of each heterogeneous configuration. Since NeuN primary antibody has a specificity for neuronal nuclear proteins, it is expressed in the nuclei of any differentiated neuron. The presence of glutamatergic neural nuclei marked with NeuN can be noticed in red (the only marker for that type of cell), while the presence of NeuN-marked GABAergic neuronal nuclei can be observed in yellow, since it overlaps with anti-GABA primary antibody (labeled in green, GABA neurotransmitter specific). Thus, the reported images indicate the right differentiation in both GABAergic and glutamatergic neurons. Then, to ensure that the nominal proportions of glutamatergic and GABAergic neurons were those actually obtained, we counted the differently labeled excitatory and inhibitory neurons. We found ratios close to the nominal values (figure 2(f)). In particular, the 75:25 configuration showed an average value of inhibition equal to 29% ± 2.5% (mean ± standard error of mean—SEM). The 50:50 exhibited a percentage of GABAergic neurons equal to 51% ± 2.6% (mean ± SEM), while the 25:75 showed 68% ± 2.3% (mean ± SEM) as average fraction of inhibitory neurons. Significant differences were found among all the configurations, thus showing a clear distinction between the E/I ratios (25:75 vs 50:50, p = 6.6⋅10⁻⁶; 25:75 vs 75:25, p = 1.1⋅10⁻¹⁷; 50:50 vs 75:25, p = 9.1⋅10⁻⁹).

Zoom In Zoom Out Reset image size

We first analyzed the spontaneous electrophysiological activity of the homogeneous networks composed by GABAergic neurons (i.e. 0:100 configuration). Over the entire development, the cultures showed a tonic activity, and an increasing signal amplitude (A_{0:100_DIV35} = 47.30 ± 0.08 μV (mean ± SEM); A_{0:100_DIV98} = 86.19 ± 0.10 μV (mean ± SEM)) in terms of single channel activity was observed (figures 3(a) and (c)), revealing a percentage of variation (Var%_35vs98) equal to +82.20%. At first glance, the electrophysiological activity of the same representative network at DIV 35 (figure 3(b)) and DIV 98 (figure 3(d)), displayed an increase in the firing at the network level as emerges by looking at the average firing rate profile. Then, we computed the MFR to quantitatively evaluate the evolution of the network firing rate. The MFR showed statistically significant differences with increasing values during the weeks (figure 3(e), MFR_{0:100_DIV35} = 1.13 ± 0.73 sp s⁻¹, n = 9; MFR_{0:100_DIV98} = 8.15 ± 5.04 sp s⁻¹, n = 5; p = 0.008) with Var%_35vs98 = +621%. Moreover, the neuronal activity was not characterized by organized bursting or network bursting patterns which were completely missing. For this reason, in figures 3(f) and (g) we only show the extreme time points, i.e. DIV 35 and DIV 98, to avoid redundant information. The MBR and NBR showed very low values, even close to zero at both DIV 35 and DIV 98 (figures 3(f) and (g)), demonstrating that few bursts or network bursts were detected during the development, even in more advanced temporal steps (MBR_{0:100_DIV35} = 0.11 ± 0.27 bursts min⁻¹, MBR_{0:100_DIV98} = 0.75 ± 0.75 bursts min⁻¹, NBR_{0:100_DIV35} = 0.11 ± 0.27 NBs min⁻¹, NBR_{0:100_DIV98} = 0.11 ± 0.16 NBs min⁻¹). Consequently, most spikes were classified as RS, i.e. not belonging to bursts (figure 3(h), RS_{0:100_DIV35} = 98.83 ± 2.61%, RS_{0:100_DIV98} = 97.25 ± 1.79% with Var%_35vs98 = −1.6%). Given the lack of bursting and network bursting activity in the 0:100 configuration, we did not consider it for further analysis.

Zoom In Zoom Out Reset image size

We monitored the activity of our cultures from DIV 35 up to DIV 98, focusing on parameters related to the spiking and bursting activity (figure 4). In the heterogeneous conditions (i.e. 25:75, 50:50, 75:25), we observed an increase in the MFR during the cultures' growth regardless of the presence of the inhibitory component. Across the network development, the MFR displayed a monotonous increase (figure 4(a)). Considering DIV 35 and DIV 98 (first and last temporal checkpoints), the networks showed a percentage of variation (Var%_35vs98) equal to +1840% in 25:75 (n₃₅ = 7, n₉₈ = 6), +1560% in 50:50 (n₃₅ = 8, n₉₈ = 5) and +1230% in 75:25 (n₃₅ = 8, n₉₈ = 5). Simultaneously, an increase in the MBR was observed as well (figure 4(b), Var%_35vs98 +85% in 25:75, +100% in 50:50, +250% in 75:25). Moreover, a temporal redistribution of the spiking activity was highlighted by a reduction of both BD and RS values (figures 4(c) and (d). Focusing on the BD, in the heterogeneous configurations a trend with ascending and then descending phase was observed. In the 25:75 and 75:25 a maximum value was reached at DIV 49 (BD_{25:75_DIV49} = 1390 ± 480 ms and BD_{75:25_DIV49} = 2124 ± 624 ms). On the other hand, the 50:50 showed a slower increase in BD profile, reaching the highest value at DIV 63 (BD_{50:50_DIV63} = 1567 ± 262 ms). Then, the BD profile decreased in all configurations until DIV 98, reaching the final values of: BD_{25:75_DIV98} = 725 ± 65 ms, BD_{50:50_DIV98} = 853 ± 104 ms, and BD_{75:25_DIV98} = 776 ± 116 ms. On the other hand, the fully excitatory networks (100:0, n₃₅ = 6, n₉₈ = 5) showed a statistical increase in MFR only during the two first weeks, reaching an average value of MFR_{100:0_DIV49} = 6.25 ± 5.6 sp s⁻¹. After DIV 49, the MFR_100:0 values did not show any significant variation. Moreover, in 100:0 configuration there was a reorganization of the electrophysiological activity induced by a temporal redistribution of the spikes: regardless the constant values of MFR (figure 4(a)), it can be qualitatively observed an increase in the MBR (Var%_35vs98 = +210%), indicating a reorganization of spikes into burst. This spikes' redistribution was also demonstrated by a constant trend of the firing in front of the reduction in BD (Var%_35vs98 = −44%) and RS (Var%_35vs98 = −60%). All the statistical differences are reported in the supplementary tables 2 and 3.

Zoom In Zoom Out Reset image size

To define the best temporal step to compare the 4 configurations, we fitted the MFR profile for all the configurations with the Hill function (figure S1). We evaluated the accuracy of the fitting model by computing the coefficient of determination (R²). The fitting indicated a relationship between the model and the data for the heterogeneous configurations with R² = 0.93, 0.95, 0.96 for 25:75, 50:50, and 75:25, respectively. On the other hand, the homogeneous configurations showed R² = 0.65 for 0:100 and R² = 0.75 for 100:0. We considered the temporal step in which the fitted curve reached 80% of the plateau as the day of maturation (i.e. corresponding to a stable level of activity). In particular, we obtained: DIV 63 for the 0:100, DIV 70 for the 25:75, DIV 59 for the 50:50, DIV 67 for the 75:25, and DIV 67 for the 100:0. To compare the configurations at the same time point, we chose DIV 70, which was the highest obtained value. Starting with the evaluation of the MFR, the heterogeneous cultures (i.e. 25:75, 50:50, 75:25) did not show statistically significant differences due to the presence of different inhibitory percentage (figure 4(e), MFR_25:75 = 20.61 ± 4.15 sp s⁻¹, MFR_50:50 = 21.60 ± 4.80 sp s⁻¹, MFR_75:25 = 25.21 ± 1.67 sp s⁻¹). The pure excitatory neuronal networks showed a lower firing activity with respect to all the heterogenous configurations (figure 4(e) MFR_100:0 = 5.81 ± 4.01 sp s⁻¹, 25:75 vs. 100:0 p = 0.004, 50:50 vs. 100:0 p = 0.006, 75:25 vs. 100:0 p = 0.009). On the other hand, MBR showed decreasing values with the increasing percentage of the excitation (figure 4(f), MBR_25:75 = 8.17 ± 0.72 burst min⁻¹, MBR_50:50 = 6.18 ± 0.27 burst min⁻¹, MBR_75:25 = 5.35 ± 0.84 burst min⁻¹). Statistical differences arose between 25:75 and 50:50 (p = 0.006), and between 25:75 and 75:25 (p = 0.003). Furthermore, the 100:0 firing activity reduction was reflected in a reduced bursting activity with a value of MBR_100:0 = 3.20 ± 2.74 burst min⁻¹. Subsequently, we investigated how the inhibition percentage affected the BD. Considering the heterogeneous configurations, although the 50:50 did not show statistical differences with the other two, a trend could be identified: the value of BD decreased with the increase of the percentage of the inhibition into the neuronal cultures (figure 4(g), BD_25:75 = 1151 ± 149 ms, BD_50:50 = 1319 ± 154 ms, BD_75:25 = 1578 ± 290 ms, 25:75 vs. 75:25, p = 0.016). The BD_100:0 deviated from the heterogeneous configuration trend and showed values with a higher spread (figure 4(g), BD_100:0 = 931.54 ± 416.50 ms) and statistically significant difference between 100:0 and 75:25 (p = 0.01). By analyzing the distribution of the spikes in the bursting pattern, the percentage of RS among heterogeneous configurations decreased with the increase of the percentage of the excitation (figure 4(h), 25:75 vs. 50:50, p = 0.005, 25:75 vs. 75:25, p = 0.006). On the other hand, the pure excitatory networks differed again from the heterogeneous configurations with statistically significant differences, showing a higher presence of RS (figure 4(h), RS_25:75 = 16.86% ± 2.85%, RS_50:50 = 11.33% ± 2.34%, RS_75:25 = 9% ± 0%, RS_100:0 = 65.50% ± 23.42%. 100:0 vs. 25:75, p = 0.003, 100:0 vs. 50:50, p = 0.004, 100:0 vs. 75:25, p = 0.008). To summarize, we can state that the completely excitatory cultures represent an extreme configuration with a particularly different behavior with respect to the heterogeneous configurations for what concern all the parameters analyzed. We can also assert that the percentage of inhibition does not affect the firing activity, since the MFR is comparable across the heterogeneous configurations. However, the inhibition component affects the activity at bursting level, decreasing the duration of bursts and increasing the frequency of burst events and the percentage of RS.

In this section, we focused on the network bursting activity. In figure 5, we reported the raster plots (top panel) and the cumulative STHs (bottom panel) of four representative neuronal cultures at DIV 70, one for each considered configuration (i.e. 25:75, 50:50, 75:25 and 100:0). At first glance, differences can be observed. First, in the heterogeneous configurations (figures 5(a)–(c)), the frequency of the occurrence of network events decreases as the level of excitation increases. Second, the NBD increases with the reduction of inhibition into the networks. Third, the 100:0 configuration (figure 5(d)) showed a lower activity amplitude and a completely different shape of the cumulative STH than the heterogeneous configurations. To quantify these claims, we extrapolated the NBR (figure 5(e)) and the NBD (figure 5(f)). The heterogeneous configurations show a decreased NBR with the reduction of the percentage of the excitatory neurons, reflecting the MBR trend (figure 5(e), NBR_25:75 = 8.16 ± 0.74 NBs min⁻¹, NBR_50:50 = 5.79 ± 0.81 NBs min⁻¹, NBR_75:25 = 5.12 ± 1.16 NBs min⁻¹, 25:75 vs. 50:50, p = 0.006, 25:75 vs. 75:25, p = 0.003). Despite there are not statistically significant differences between the NBD across the heterogeneous configurations (figure 5(f)), also for this parameter we can observe the same increasing trend of the BD when the amount of excitation increases in the neuronal networks. The 100:0 configuration shows NBR and NBD values (NBR_100:0 = 6.40 ± 2.06 NBs min⁻¹, NBD_100:0 = 2.31 ± 1.20 s) that deviate from the trend of heterogeneous configurations, as seen above for other parameters. In figure 6(a), we reported all the network bursts aligned with respect to the peak of the longest network burst of a representative culture for each configuration at DIV 70 with the respective cumulative STH (overlapped). Despite we chose representative neuronal networks, the n = 32 cultures are characterized by repeatable neuronal network activity, since the 90% confidence interval of the average cumulative STH is restricted for all the considered configurations (figure S2). In figure 6(b), we reported the average cumulative STH of each configuration obtained by averaging the cumulative STHs of each culture belonging to the same configuration. From a qualitative point of view, the heterogeneous cultures show a different amplitude (STH_peak) of the cumulative STH and, again, a diverse NBD, in particular the decay phase duration (figure 6). These qualitative considerations are supported by statistically significant differences: the STH peak amplitude is inversely proportional to the percentage of GABAergic neurons present in the neuronal cultures (figure 6(b), inset. STHpeak_25:75 = 10 367 ± 1497 sp s⁻¹, STHpeak_50:50 = 13 445 ± 786 sp s⁻¹, STHpeak_75:25 = 16 049 ± 1161 sp s⁻¹; 25:75 vs. 50:50, p = 0.005; 25:75 vs. 75:25, p = 0.003; 50:50 vs. 75:25, p = 0.011). The time constant of the rise dynamics (τ_rise, figure 6(b), inset) shows comparable values for the heterogeneous configurations (τ_{rise 25:75} = 2.01 ± 0.69 ms, τ_{rise 50:50} = 1.35 ± 0.64 ms, τ_{rise 75:25} = 1.72 ± 0.30 ms), while the τ_decay shows decreasing values with the increased percentage of the inhibition into the heterogeneous neuronal networks (τ_{decay 25:75} = 66.21 ± 6.45 ms, τ_{decay 50:50} = 83.15 ± 26.21 ms, τ_{decay 75:25} = 102.48 ± 29.40 ms, 25:75 vs. 75:25, p = 0.017). On the other hand, the 100:0 shows different shape of the curve, a reduced NBD and a drastic reduction of the amplitude of the cumulative STH (figure 6). These qualitative considerations are supported by the statistically significant differences of the STH peak amplitude between the 100:0 and the heterogeneous configurations (figure 6(b), inset. STHpeak_100:0 = 5148 ± 4227 sp s⁻¹; 100:0 vs. 25:75, p = 0.045, 100:0 vs. 50:50, p = 0.018 100:0 vs. 75:25, p = 0.004). Furthermore, the 100:0 reveals higher τ_rise values (τ_{rise 100:0} = 54.13 ± 28.42 ms, 100:0 vs. 25:75, p = 0.009, 100:0 vs. 50:50, p = 0.006 100:0 vs. 75:25, p = 0.005) and a minimum value of τ_decay (τ_{decay 100:0} = 37.09 ± 11.91 ms, 100:0 vs. 25:75, p = 0.003, 100:0 vs. 50:50, p = 0.004 100:0 vs. 75:25, p = 0.003). To summarize, even at the network level activity, the 100:0 configuration shows a different behavior from heterogeneous networks. On the other hand, the amount of inhibition in the heterogeneous configurations affects the frequency, the duration (in particular the decay dynamic) and the amplitude of the network-wide events.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the previous sections, we found that the different percentages of inhibitory neurons plated at DIV 0, which modulate differently the network activity, are stable over development from a morphological point of view. In order to verify whether these percentages were found and maintained also from a functional point of view, we investigated the functional connectivity using the TSPE algorithm (see section 2.4.3). Figure 7(a) shows a representative example of the functional connectivity map of a physiological neuronal network (75:25). Each node represents an electrode of the MEA, while links indicate the detected functional connections. The topological characterization of the inhibitory networks (0:100) was not performed as the low spiking and almost absent bursting activity did not allow the extraction of robust cross-correlograms. In fact, considering DIV 70, the connections detected in inhibitory networks represent only 5% of the total functional connections detected in all configurations (figure 7(b)). Moreover, due to the limit imposed by the algorithm for the detection of causal connections from inhibitory to inhibitory and from inhibitory to excitatory neurons in a non-random network, we decided to also neglect the networks of the 25:75 configuration. In fact, due to its unbalanced composition in favor of the inhibitory component, this configuration presents many of these wrongly detected connections. Figure 7(c) quantifies the percentage of excitatory connections that modulate the neuronal activity in the different configurations. We detected 69.6 $\pm$ 10.4% excitatory connections in the 50:50 configuration, 80.32.3% in the 75:25 configuration, and approximately 83.8 $\pm$ 2.3% in the fully excitatory (100:0) configuration. The reported values deviate from the putative percentage of excitatory neurons (figures 2(d) and 7(c)); nonetheless, it is worth noticing how the percentage of functional excitatory connections increased as a function of the increase of excitatory neurons. The 50:50 networks displayed a percentage of excitatory connections lower than the 75:25 condition of about 17.0% (p = 0.004), and lower than the 100:0 configuration of about 13.3% (p = 0.002). Moreover, the 75:25 configuration showed 4.2% fewer excitatory connections than the 100:0 configuration (p = 0.03).

Zoom In Zoom Out Reset image size

Since the inhibitory component of neuronal networks is considered responsible for the synchronization of network activity (Bonifazi et al 2009), it was interesting to study how both the observed functional and structural variations at the core of this study modulated the global temporal organization of the events. We evaluated the synchronization level (CI₀) of the network in terms of both spiking (figure 7(d)) and bursting (figure 7(e)) activity. From the spiking activity point of view, a descending trend of synchronization level was observed as a function of the decrease of inhibitory cells that composed the networks. In fact, inhibitory networks (0:100), which are characterized by tonic activity, showed an average synchronization value of 0.03 $\pm$ 0.0003, which was higher than all the other configurations. Focusing on the heterogeneous configurations, there were no appreciable statistical variations between 25:75 and 50:50 (p = 0.05). Instead, 25:75 showed a sync value higher than 75:25 by 4.1% (p = 0.001). A reduction of 1.3% in the level of synchronization was observed between the 50:50 configuration and 75:25 one (p = 0.01). At this point, a sort of plateau in the synchronization of the spiking activity was reached: purely excitatory networks (100:0) had comparable values with the 75:25 ones (p = 0.5) and lower synchronization values than all other configurations. A different behavior was observed at the bursting activity level (figure 7(e)): although not always statistically relevant, an increase in the synchronization value was observed as the excitatory contribution increases, indicating a greater network synchronization (network burst) in cultures with a higher percentage of excitation.