Simulation of developing human neuronal cell networks

Background Microelectrode array (MEA) is a widely used technique to study for example the functional properties of neuronal networks derived from human embryonic stem cells (hESC-NN). With hESC-NN, we can investigate the earliest developmental stages of neuronal network formation in the human brain. Methods In this paper, we propose an in silico model of maturating hESC-NNs based on a phenomenological model called INEX. We focus on simulations of the development of bursts in hESC-NNs, which are the main feature of neuronal activation patterns. The model was developed with data from developing hESC-NN recordings on MEAs which showed increase in the neuronal activity during the investigated six measurement time points in the experimental and simulated data. Results Our simulations suggest that the maturation process of hESC-NN, resulting in the formation of bursts, can be explained by the development of synapses. Moreover, spike and burst rate both decreased at the last measurement time point suggesting a pruning of synapses as the weak ones are removed. Conclusions To conclude, our model reflects the assumption that the interaction between excitatory and inhibitory neurons during the maturation of a neuronal network and the spontaneous emergence of bursts are due to increased connectivity caused by the forming of new synapses.

external input or stimulus, as observed in MEA experiments. Each neuron has either an inhibitory (negative synaptic strength) or an excitatory (positive synaptic strength) effect on its neighbors. These models of synaptic communications can be considered to include all types of interactions between the neurons. The activity of a neuron depends on its previous spiking history.
Here, the INEX model is used to simulate the developing hESC-NNs on MEAs. The model and its parameters are tuned to mimic the activity measured from in vitro hESC MEA data from six measurement time points during neuronal network activity development and maturation. The activity level is defined as various spike and burst parameters. Thus, the modeled neuronal networks will produce statistically similar spike and burst activity as the in vitro actual neuronal system. Therefore, the main question we aim to answer with the simulations is: Which aspects of the maturation process contribute to the development of stable burst patterns?

Cell cultures
Human embryonic stem cells (hESCs) [cell lines Regea 08/023 and 11/013] were differentiated into neuronal cells using the previously published method [9] and plated on MEAs as described in Heikkilä et al. [6]. Briefly, cells were differentiated for 8 weeks in differentiation medium containing D-MEM/F-12 and Neurobasal (1:1, both from Gibco Invitrogen, Carlsbad, CA, USA), N2 supplement (Gibco Invitrogen, Carlsbad, CA, USA), B27 supplement (Gibco Invitrogen, Carlsbad, CA, USA), 2 mM GlutaMax (Gibco Invitrogen, Carlsbad, CA, USA), and 25 U/ml penicillin/streptomycin (Cambrex, Verviers, Belgium) in the presence of basic fibroblast growth factor (4 ng/ml, FGF, Sigma-Aldrich, St. Louis, MO, USA) in neurosphere culture. Next, 10-15 small aggregates dissected from neurospheres (50,000-150,000 cells in total) and plated to MEA or dissociated into single cell suspension using TrypLe Select (Sigma-Aldrich, St. Louis, MO, USA) and thereafter plated on MEA dishes. The dishes were coated with polyethyleneimine (0.05 % solution, Sigma-Aldrich, St. Louis, MO, USA) and subsequently with mouse laminin (20 μg/ml, Sigma-Aldrich, St. Louis, MO, USA). Differentiation medium supplemented with FGF (4 ng/ml) and brain-derived growth factor (5 ng/ml, BDNF, Gibco Invitrogen, Carlsbad, CA, USA) was replaced three times a week for the MEA cultures. All the MEAs with cells were kept in an incubator (+37 °C, 5 % CO 2 , 95 % air) prior to and between recordings. All recordings were made using MEAs and equipment purchased from Multi Channel Systems (MCS GmbH, Reutlingen, Germany). Figure 1 shows the neuron distribution on 7, 12, and 19 days in vitro (DIV) in MEAs. In addition, cultures grown on the cell culture well plates, were stained with Gamma-aminobutyric acid (GABA) antibody (Rabbit anti-GABA IgG, 1:1000, Sigma Aldrich, St. Louis, MO, USA). Cells were calculated from at least two wells, at least five images and repeated at least twice for each different measurement time point. Additionally, a portion of the cultures were stained either with neuronal marker Mouse anti-β-tubulin III IgG (1:1200, Sigma Aldrich, St. Louis, MO, USA), with GABA synthesizing enzyme glutamate decarboxylase Mouse anti-GAD67 IgG (1:100, Chemicon International Inc., Temecula, CA, USA) or with calcium-binding protein calretinin Rabbit anti-calretinin IgG (1:800, Swant, Marly, Switzerland). The immunocytochemical protocol has been published previously [9]. hESC experiments were performed at the Institute of Biomedical Technology (University of Tampere, Tampere, Finland). Approval was given to culture the hESC lines (Skottman, R05116) by the Ethics Committee of the Pirkanmaa Hospital District.

Electrophysiological recordings
Electrical activities were recorded using two 1-well (60MEA200/30 Ti, datasets #8 and #9) and eight 6-well MEAs (60-6wellMEA200/30iR-Ti-w/o; all from MCS GmbH, Reutlingen, Germany). All MEAs had internal reference electrodes. Signals were sampled at 20 or 50 kHz, and stored on a standard PC using MC Rack software (MCS GmbH, Reutlingen, Germany). During the measurements, the culture temperature was maintained at +37 °C using a TC02 temperature controller (MCS GmbH, Reutlingen, Germany). Recordings were visually inspected for artifacts and the measurements or channels likely to contain artifacts were excluded from further analysis.
MEA recordings from ten hESC-NNs were used with an approximated spike train (sequence of spikes) length of 300 s. The hESC-NNs were measured as follows: the first measurement time point was at 7 DIV when the neurons in at least 10 % of the channels of the MEA were active, and when at least 100 spikes within the 300 s were found in the active channels during the recording period. To make the hESC-NN datasets #1-#10 comparable, they were grouped according to the measurement time points (MTP) 1-6, which correspond to 7-26 DIV (see Table 1). The spontaneous activity developed by the hESC-NNs is important in neural development and includes differentiation, maturation, and generation of neuronal processes and connections [6,9]. Channels were considered as inactive when less than 20 spikes/min [10] were recorded at the last measurement time point (measurement time point 5 or 6). Additionally, if less than two channels per well were active, the well data were excluded from further analysis.
To get a reference for the simulation, we calculated the medians and lower and upper quartiles of spike rate, burst rate, burst duration, and average number of spikes per burst separately for all electrodes and all measurement time points, as shown in Fig. 3. Briefly, the burst analysis algorithm, which was used to examine the intrinsic bursting, relies on the cumulative moving average (CMA) and the skewness (α) of the interspike interval (ISI) histogram. For bursting, the ISI threshold was found at the ISI closest to the value of α · CMA m , where CMA m is the average of CMA. Additionally, three or more spikes  Table 1) on the MEA for three points in time (a 7 days in vitro (DIV), b 12 DIV, and c 19 DIV). It is clearly visible that the number of neuronal connections increases and the neurons move over time. The black dots indicate the MEA electrodes. The scale is 100 μm had to be in a row. The CMA algorithm does not use a fixed ISI but adapts to the dynamics of the studied spike trains. Burst duration means the time between the first spike's peak and the last spike's peak. Kapucu et al. [10] have demonstrated the functionality of the tool for highly variable network structures and time-varying dynamics such as in hESC-NNs. In 78 % of all electrodes, the spike rate increased from measurement time point 1 to measurement time point 5. In 16 % of electrodes it decreased, and in 6 % it remained stable or zero. In 70 % of all electrodes, the burst rate increased from measurement time point 1-6. In 20 % of electrodes it decreased, and in 10 % it remained stable or zero. The datasets showed a large variability. For model validation, the means of spike rate and burst rate per well were calculated. The wells were grouped according to the spike rate at measurement time point 5 in low (<50 spikes/min), medium (between 50 and 250 spikes/min), and high (>250 spikes/min) activity (Table 2). This is a kind of normalization to be able to compare the measurements. To get some similarity of the varying cultures, we used only the medium activity datasets for the analysis and simulations. Figure 3 displays the development of neuronal network activity in the medium range.

Table 1 Sorted measurement time points (MTP) of the cultured hESC-NNs
The sign × means no measurement was done on this measurement time point. The first MTP was on the 7th day in vitro (DIV). MTP 2 was between 9 and 12 DIV, MTP 3 was between 14 Even if the spike rate and the burst rate showed high variability, the general tendency in both features is an increase.

INEX model
To simulate the maturing hESC-NN, we used our spiking neuronal model called INEX [24]. Briefly, the phenomenological model is a cellular automaton whose cells are neurons with two possible states: ON or OFF. Each neuron obtains several inputs and produces exactly one output (spike or no spike). In order to simulate spontaneous activity, we assume that the spikes obey an inhomogeneous Poisson distribution [25]. The momentary firing rate i of neuron i in time slice t k is calculated as follows: where c i denotes the basic activity (which include all kind of noise sources such as thermal noise), y ji the synaptic strength of all neurons j connected to neuron i and s j the particular spike of the previous time slice of neuron j (1 for a spike and 0 for no spike). To find appropriated values for the parameters types c i , y + ji and y − ji , a brute force approach was used. The parameter values were randomly chosen from a triangular distribution. The values lie between zero and an upper boundary that is at most 1. For c i , the upper boundary varies from 0.01, 0.02, …, 0.09, for the excitatory synaptic strength y + ji from 0.1, 0.2, …, 0.9 and for the inhibitory synaptic strength y − ji from −0.1, −0.2, …, −0.9. For the evaluation of the parameter space search, the mean values of the basic activities and synapses strengths of all neurons were calculated. The objective functions of the parameter space search are the spike and burst rate obtained from the experimental data. This means that they are approximately in the range of the MEA data (see Table 3). The brute force method was applied to the simulated data of each virtual measurement time points (vMTP). The vMTPs are considered to resemble the actual measurement time points.
The probability P i for the occurrence of a spike in time slice t is defined as follows: The time slice t is chosen with a length of 5 ms to cover the temporal length of the action potential and the subsequent refractory period. For each time slice, the algorithm tests if x i < P i , where x are uniformly distributed random values. A spike time history was added so that x decreases by factor f = 0.1 when a spike in the last time slice occurred. This history method corresponds biologically to an enforced state of excitation of a neuron after spiking and ensures synchronous bursting of the neurons. To summarize, our model has four parameters: c i , y + ji , y − ji , and f.

Simulation of maturing neuronal networks
In our in vitro MEA experiments with hESC-NN, circa 50,000 to 150,000 cells were plated on each well. Based on calcium imaging assessment (data not shown) an estimated 1000-4000 neurons were active and could be recorded. Based on these findings, we chose to simulate 1000 neurons. In the MEA data, one electrode signal is the sum of activity of possible one or several neurons detected by the electrode. In the INEX model, we can consider that one computational neuron corresponds to the activity shown by one electrode. Thus, the model depicts the activity seen by the measurement system as in many other neuronal network models [22,23]. In the brain, the common proportion of excitatory pyramidal cells and inhibitory interneurons is considered to be 80 and 20 %, respectively [26]. The inhibitory interneurons are mainly GABAergic neurons (reviewed by Chattopadhyaya et al. [27]). The proportion of GABAergic cells in hPSC-derived neuronal cultures has not been studied to any great extent but, based on the immunocytochemical analysis, the portion of GABA positive cells varies between 35 and 90 %, depending on the differentiation method used [28][29][30]. Here, we performed GABA analysis of cultures paralleling the measurement time points. The portion of GABA positive cells varied between 13 and 19 % of the total neuronal cells (Fig. 4). Thus, for the simulation model, we used the common proportion of 80 % of excitatory neurons and 20 % of inhibitory neurons.
We assumed that there are no connections between neurons on the day of plating and no autapses [31,32], which are self-connections of a neuron. The INEX model only allowed the addition of connections. Therefore, no reduction of connections [11] was

Table 3 Lower quartile (Q1), median (M) and upper quartile (Q3) of the calculated features for simulated (INEX) and experimental (MEA) data on measurement time point (MTP) 1-6
The table shows as well the selected upper boundary as result of the parameter space search which resembled best the spike and burst rate of the experimental data (see "Methods" section/ INEX model). The table is visualized in Fig. 3 The spike rate SR is given in spikes per minute, the burst rate BR in bursts per minute, the burst duration BD in seconds and the average number of spikes per burst SB as count. simulated. Connections appeared simultaneously between two sequential vMTPs. The model did not take into account apoptosis or proliferation, and we did not include transmission delays or cell movement in the model. In order to model the maturing process and the developing connectivity of the neural network, we started with a few randomly chosen connections with a probability of 1 % of all possible connections and weak synaptic strength for vMTP 1, respectively. Thus, the neuronal network was not inactive at the first simulation step (vMTP 1). Then, the connection probability was increased to 2, 4, 6, 8 %, and up to 10 % of all possible synaptic connections (corresponding to vMTP 2 to vMTP 6) [22]. The 10 % connection probability corresponded to the connection probability in matured neuronal networks. The arrangement of connections between the neurons was selected randomly. For every vMTP, the connections in the simulated neuronal network were redefined. The values of the synaptic strengths were automatically varied with a brute force approach, as presented above. Additionally, we simulated according to the following scenario: (1) an increase of the activity between vMTP 1 and vMTP 6; (2) an increase of the activity between vMTP 1 and vMTP 5, and a decrease at vMTP 6, as seen in Fig. 3. All resulting spike trains had a length of 300 s. The simulation tool was then run ten times with these constraints to get statistically significant data.

Validation of the simulated spike trains
For validation, we calculated four features [spike rate (spikes/minute), burst rate (bursts/ minute), burst duration (in seconds), and the average number of spikes per burst] for each of the simulated spike trains using the burst analysis tool described by Kapucu et al. [10]. The results were then compared with the same features obtained from the ten previously mentioned MEA experiments with hESC-NNs. The spike rate and burst rate were selected as goal functions for the parameter search. Too many features would lead to an over-fitting, and thus produce unstable points. The other two parameters, burst duration and average number of spikes per burst, described the burst structure and seemed to undergo typical changes during the maturation of the network.

Results
As a basis for our simulations, we conducted 10 MEA experiments (two 1-well MEAs each with 60 electrodes and eight 6-well MEAs each with nine electrodes) with hESC-NNs. The datasets were grouped according to six measurement time points that correspond to 7-26 days in vitro in MEAs ( Table 1). The INEX model generated a large-scale network of 1000 neurons that corresponds to the number of active cells in the experiments with hESC-NNs. For the vMTP 1-6 used in the simulations, we created a neuronal network with increasing connection probability over time. We applied a brute force method to each obtained dataset to find one parameter set (comprising the basic activity, the excitatory and inhibitory synaptic strengths, and a factor for the spike time history) that produced neuronal activity that best resembled the experimental data.
We kept the basic activity, which was modeled as the random noise of each neuron in the system, as constant as possible for vMTP 1-6 with the hypothesis that during maturation only the network properties will change. Thus, only the inhibitory and excitatory synaptic strengths were more variable (in comparison to the basic activity which remains stable over the measurement time points). The simulated network showed an increase in the excitatory synaptic strengths over time (Table 3). This increase continued until the final vMTP where a decrease in excitatory synaptic strengths was observed. The inhibitory strengths remained stable over the simulated time duration. For each vMTP, we simulated ten datasets, each with 1000 neurons. For the first nine neurons (corresponds to the number of electrodes on a 6-well MEA), we calculated the lower and upper quartile as well as the median of four features, in particular spike rate, burst rate, average number of spikes per burst, and burst duration. Table 3 and Fig. 3 show both the development of the four features from measurement time point 1-6 for both experimental and simulated data. The validation showed that all calculated median values of the spike rate in the INEX data are within the lower and upper quartile of the MEA data. This was also the case for the burst rate with the exception of vMTP 6. Nevertheless, the upper quartile of the simulated data was within the quartile range of the experimental data. In three out of six measurement time points, the median and quartiles of the burst duration in the simulated data were higher than in the MEA data. The median of the average number of spikes per burst was mostly within the quartile range of the experimental data. For the spike and burst rate as well as for the average number of spikes per burst, we saw an increase in the features over time in the experimental data and correspondingly in the simulated data. The spike rate and burst rate dropped at the last measurement time point in the experimental, and thus also in the simulated data. The alternating burst duration over maturation can be seen in both the experimental and simulated data.
The spike trains of five sample electrodes and five example neurons are displayed in Fig. 2a. The experimental and simulated spike trains of the first measurement time point showed just a few spikes. The overall number of spikes increased with the number of connections and with the number of measurement time points (Figs. 2a, 3). The simulated activity of the last measurement time point exhibited typical spike and burst patterns as recorded from hESC-NNs (see Table 3) [6]. Partly synchronous spiking and intrinsic bursting was recorded for maturated hESC-NNs and could also be seen in the corresponding simulated spike trains. Figure 2a also displays the raw voltage traces of channel 63 of the same hESC-NN. Figure 2b shows the ISI histograms of one experimental and one simulated neuron at measurement time point 5. Both histograms show a similar ISI distribution. By varying the inhibitory and excitatory parameters, the model produced similar spiking characteristics to those measured. Figure 2b displays also the population ISI histograms of dataset #9 and one simulated neuronal network at (v)MTP 5.

Stem cell data
The potential of human pluripotent stem cells and their neural derivatives in the fields of neurotoxicity, drug screening, developmental biology, and tissue engineering is well known [1,2,33]. In these applications, stem cells need to be differentiated into pure neuronal populations and show neuronality in genotype and phenotype as well as at the functional level [33]. Thus, it is also important to study these cells in vitro at the functional level [34]. MEAs are used for the characterization of the network activity of these cells as well as for studying drug and neurotoxic effects on the cells [6,8]. However, little is known about the development of the network processes that generate the signaling patterns in hESC-NN. Earlier, Heikkilä et al. [6] observed single spike activity in hESC-NN cultured on MEA in the first week followed by the development of spike trains during the next two weeks. From the fourth week onwards, they observed synchronous bursts. Our study had similar results (see spike trains and voltage traces in Fig. 2 and the statistics in Fig. 3) with the exception that the data points used were up to 26 DIV, and thus later points of network maturation were not studied. Here, as a larger data set was analyzed, we found quite high variability in the spike and burst behavior across network maturation. The observed variability can be explained by the different number of cells in the networks and the various fractions of neuronal and glial cells on these spontaneously formed neuronal networks. Furthermore, there is evidence that the neuronal networks are not fully matured even at measurement time point 5 or 6, which corresponds to 21-26 DIV, respectively, and that the networks we used are still in different developmental stages [6,35], since the signaling of these measurement time points differ from others in terms of both spike and burst behavior.
In addition to synaptic activity, several other activity pathways exists especially during development [36]. Especially, gap junction mediated activity is important during development [37] and has been also studied in dissociated neuronal networks cultured on MEAs [20]. In this work, we focus only on synaptic mediated activity, which exists in these human neuronal cultures as proved with pharmacological modification of neurotransmitter receptors [6].
For the burst analysis, we did not use the traditional burst analysis approach with fixed ISI that had been used earlier with similar cultures (e.g., Heikkilä et al. [6]). As Kapucu et al. [10] demonstrated, the traditional approach quite often fails when examining hESCs. Thus, the authors developed the cumulative moving average approach that adapts the ISI threshold for bursts to the network behavior [10]. The method also  Table 3 Lenk et al. BioMed Eng OnLine (2016) 15:105 finds statistically burst-like behavior in spike data from spike trains with quite low firing activity. Here, we use the CMA tool for analysis on both simulated and measured data, resulting in comparable statistical data. The synchronous population bursts behavior described earlier by Heikkilä et al. [6] was not taken into consideration because the used data sets did not cover later time points (1 month onwards).
The field of in vitro experiments with hESC-NNs is quite new and not all of the previously conducted experiments were suitable as a basis for our simulations because we modeled the maturation over a relatively long time period. Even with a limited number of data sets, we can see the tendency of first an increase and later a decrease in neuronal activity, especially in the spike and burst rate (see Fig. 3). Johnson et al. [38] also report that the neuronal activity is reduced during the course of the maturing process.
The in vitro cultures are meant to mimic the neuronal network in vivo. Even the in vitro developed neuronal network might lack certain network structural functions as seen in brain and possible effect, like the electrical field effect, between the neurons might not be observed in the cultured neurons [39]. However, the hESC-NNs provide us a way to model in vitro the human neuronal system that has bot been available earlier.

Simulation
The INEX model is a very simple, general and flexible model. Despite its primary application for cortical culture modeling [24], it is not bound only to the simulation of cortical networks in vitro. In this study, we use large-scale networks with 1000 neurons to study spike and burst behavior in hESC-NNs. Here, the neurons are considered as points with neither spatial extension nor bio-physical structure (no axons, soma, or dendrites) and the connections representing synapses are formed randomly between these virtual neurons. We made a number of simplified assumptions which are described in the "Methods" section. There is evidence that neurons are interacting with synapse communication, with gap junction mediated exchange of ions and small molecules, like ATP (adenosine triphosphate), and with electric field effects [40]. Computational models of neuronal networks simulate the synaptic transmission per se; however, we can consider that the model of interaction includes all communication as the parameters of the synaptic interaction models are tuned to provide similar responses as in the actual biological networks. Further, as the communication through the other pathways is not directly mediated by spiking activity modeling such weak and less known pathways, it is not considered the core of the present study. It has been previously shown that these networks and cultures have a minority of astrocytes [9]. The INEX model does not take glial cell effect directly into account. However, the effect is inbuilt with the effect of the spike history.
It is assumed that synapses develop during the maturation process, and that mature in vitro networks have a connectivity of about 10 %. This means that each neuron is connected to 10 % of the other neurons. For the simulation, the starting point is almost no connectivity (1 %), and the end point has 10 % connectivity [22]. The steps in between correspond to virtual measuring points and are determined linearly (1, 2, 4, 6, 8 and 10 %). In line with this, the experimental measuring points are also almost linear. Another approach would be to increase the connectivity exponentially with restricted resources as described by Lai et al. [41]. However, a detailed connectivity analysis of hESC-NN has not yet been carried out. Therefore, we did not follow this approach in this paper.
Present technology such as MEA or patch clamp cannot provide the connectivity analysis reliably, the results from the INEX model strengthen the concept that the maturing hESC-NN and its spiking activity can indeed be explained by the development of connectivity between the neuronal cells. In biological networks, the development of connectivity can be generally explained either as increased synaptic strengths, increased number of synapses between the processes, or an increased number of processes between cells [6,14]. Nevertheless, an overproduction of synaptic connections is followed by the elimination of some synapses and the stabilization of activity [14]. The results indicate that the model can simulate the reduction of synapses [42,43], which is an important feature of the maturating process, by changing the synaptic strengths. Thus, the number of neurons remains the same over all virtual measurement time points. Without thorough biological characterization of the time course of this development in vitro, the separation of these processes using model concepts is in practice very difficult or even impossible. Therefore, these difficulties must be taken into account when evaluating the results presented in this paper.
A steadiness or increase in the excitatory synaptic strengths is seen in the simulations from vMTP 1 to vMTP 5. At vMTP 6, the excitatory strengths are slightly reduced and result in reduced spike and burst activity, as seen in the experimental data. The inhibitory synaptic strengths, however, remain stable over time. From the simulations, we can draw the conclusion that the share of inhibitory neurons is relatively low, as the inhibitory strengths remain low. This can also be observed in the experimental data (see Fig. 4). Moreover, if the proportion of the inhibitory and excitatory neurons is incorrect, the strengths and ratio of excitatory and inhibitory neurons in the simulation can compensate for this situation. As both strengths and the number of inhibitory neurons remain low, we consider the conclusions to be correct. The calculated features adapted from spikes and bursts show that the maturing process of hESC-NNs can be modeled