article: Recently evidence has accumulated that many neural networks exhibit self-organized criticality. In this state, activity is similar across temporal scales and this is beneficial with respect to information flow. If subcritical, activity can die out, if supercritical epileptiform patterns may occur. Little is known about how developing networks will reach and stabilize criticality. Here we monitor the development between 13 and 95 days in vitro (DIV) of cortical cell cultures (n = 20) and find four different phases, related to their morphological maturation: An initial low-activity state (≈19 DIV) is followed by a supercritical (≈20 DIV) and then a subcritical one (≈36 DIV) until the network finally reaches stable criticality (≈58 DIV). Using network modeling and mathematical analysis we describe the dynamics of the emergent connectivity in such developing systems. Based on physiological observations, the synaptic development in the model is determined by the drive of the neurons to adjust their connectivity for reaching on average firing rate homeostasis. We predict a specific time course for the maturation of inhibition, with strong onset and delayed pruning, and that total synaptic connectivity should be strongly linked to the relative levels of excitation and inhibition. These results demonstrate that the interplay between activity and connectivity guides developing networks into criticality suggesting that this may be a generic and stable state of many networks in vivo and in vitro.
  During the last years increasing evidence has accumulated that networks in the brain can exhibit “self-organized criticality” [1]–[3]. Self-organized criticality is one of the key concepts to describe the emergence of complexity in nature and has been found in many systems – ranging from the development of earthquakes [4] to nuclear chain reactions [5]. All these systems exhibit spatial and temporal distributions of cascades of events called avalanches which can be well described by power laws [6]–[8]. This indicates that the system is in a critical state [6], [9] and that similar dynamic behavior exists across many different scales. Several neural network models have predicted that neural activity might also been organized this way [10]–[14] and recently this had been confirmed experimentally [1], [15]–[17]. A recent study by Levina and colleagues [18] addresses the question how self-organized criticality can emerge in such networks in a robust way by using dynamical synapses, which alter their synaptic connection strength on a fast time scale. This contribution, which is able to analytically predict the network behavior, is a theoretical milestone in our understanding of criticality in neural systems. In general, however, theoretical and experimental investigations have so far usually focused on mature networks [1], [16] sometimes including adaptive processes [18]–[20]. Little is known how developing networks can reach a final state of self-organized criticality [10], [17], [21]. In the current paper, we are therefore experimentally investigating the different stages of developing cortical cell cultures [22] to assess under which conditions these networks develop into a critical state. Specifically we are asking the following questions: 1) do the investigated cell cultures undergo a significant transition in their activity states and how is this related to self-organized criticality and 2) can specific predictions be made with respect to network activity and connectivity which would explain the observed behavior. To address the second aspect we are designing a model to simulate network development, which is based on activity-dependent axonal and dendritic growth leading to homeostasis in neuronal activity [23]–[28].
  In order to assess how self-organized criticality develops in cell cultures, we have monitored a total of 20 cultures and recorded their activity patterns between 13 and 95 days in vitro (DIV). In general, cultures start with about 500,000 dissociated cortical neurons, which develop over time into an interconnected network. To assess the different network states the activity at 59 electrodes was measured and analyzed at different DIV (see Methods). Figure 1 A shows 15 minutes of recorded activity for one typical culture at 42 DIV. At this temporal resolution individual bursts are visible as vertical dot-lines indicating activity at almost all electrodes, separated by rather long pauses which allow for robust separation of these bursts required for avalanche analysis. At a fine temporal resolution (Figure 1 B) one sees that the burst activity expresses certain patterns. Note, pauses have been graphically shortened in panels (B) and (C). Panel (C) shows the activity pattern that arises in our model, which at a first glance looks similar to that in the culture. Details about the model and an analysis which support similarity of model and real data, will be provided later. First we would like to describe the developmental stages in the cultures with respect to their avalanche distributions. In this work, avalanches are defined by the number of spikes between two windows without activity (see Methods). At early stages during development, usually before 13 DIV, connectivity is small and activity in the network very low. So, it is very difficult to obtain long enough recordings for plotting avalanche distributions. However, known from the literature [29], in this stage activity is best described by a Poisson like behavior. At about 13 DIV (see Figure 2), we receive the first distributions which develop towards criticality (Figure 3 A). Therefore, we call this state the initial state. The ideal power-law fit for each curve is shown by the dashed lines. If a distribution matches the power law line it can be called “critical” [6], [7]. A dominance of long avalanches is indicative of a supercritical state whereas a lack thereof is referred to as subcritical. This is measured by, which gives the quality of fit between ideal power law and actual distribution. For a system in a supercritical state is larger and for a subcritical state smaller than zero (see Methods). Values of are also shown in the different panels of Figure 3. For the cultures, we receive at (on average) 19 DIV values of in the interval from to. While this shows that the system develops towards criticality, we also observed that this behavior is very unstable. Quickly, within just (on average) one/two days, the distributions change shape and develop a substantial “bump” for larger avalanches. This indicates that at (on average) 22 DIV the network enters a supercritical regime (Figure 3 B). After (on average) 36 DIV network activity is curbed and it reaches a subcritical regime (Figure 3 C). This can be seen by the decrease of the distribution at larger avalanches. At (on average) 58 DIV the system becomes finally critical (Figure 3 d). Here we find that the deviation from a power law is nearly zero (for these examples). In general we find that the differences between all states are significant for the measured values of (ANOVA test). Figure 2 provides the data of all 20 cultures (see Methods) divided into the different states. All completely measured cultures undergo the same transitions from initial (black) to supercritical (red) to subcritical (green) and finally to a critical state (blue). The overlap between the first two states results from the very quick transition between them together with small differences in the speed of development of the different cultures. Average values of for these four states are, , , and (see Table in Figure 2). Differences are significant using the multiple comparison procedure with Bonferroni correction based on the one-way ANOVA test. Only the difference between the initial and critical state is not significant as in the initial state the network develops towards criticality until strong morphological changes set in (see Phase I). However, the activity given by the number of action potentials per minute is for the supercritical state significantly higher than for the initial, subcritical and critical state, which has the lowest mean activity. These were the only differences that were observed. In summary, these results show that there is a characteristic time course in the development of the avalanche distributions. The system starts with low activity and then enters a transitory initial state. Quickly it leaves this state and, passing supercritical and subcritical regimes, reaches the critical state. The first phase (Phase I) of the network development is characterized by dendritic growth to establish first synaptic contacts and to rise neuronal activities. At the beginning of the model development the dendritic acceptance increases (Figure 4 C). By this outgrowth the system creates synapses and forms a network. The distribution of the avalanches, the mean membrane potential, and the mean calcium concentration also changes (mean values over all neurons are given as upper case letters, while lower case letters indicate individual values). Similar to real cell cultures, all neurons at this phase are excitatory. With the help of a mean field approach it is possible to calculate average membrane potential and average calcium concentration during this phase. Different from real networks, where the activity is too small to render reliable measurements for very early developmental stages, in the model we can also analyze these. For this, the term, which determines the increase of the membrane potential according to the activity of the connected neurons in Equation 14 (see Methods), is simplified to a product of the mean membrane potential and an monotonous function dependent on the mean synaptic density (see below) and we get for the activity change: (1) The differential equation of the calcium concentration (Equation 15 in Methods) can be written as: (2) With these equations, we can now consider three different degrees of synaptic densities in the first phase; namely zero, small, and medium densities and for Phase II with a large density. Phase II of the network development is characterized by an overshoot in network activity. The membrane potential and calcium concentration (and) reach their maximum. This causes a phase transition in axonal and dendritic development: At that point, the dendritic acceptance begins to shrink and the axonal supply increases (see 4 C, D, Phase II). Moreover, during such transitions (accompanied with the formation of very many synapses) the action of the transmitter GABA switches from excitatory to inhibitory due to a change in the intracellular chloride concentration [30]. As we do not model changes in ion concentrations, we just change 20% of all neurons and assign them a negative value of, thereby making them inhibitory (is changed to in this second phase). To determine the influence of different degrees of inhibition, the ratio of to is chosen differently in different experiments (Figure 5 D–F). We can calculate the membrane potential as before with Equations 1 and 2 now with the constraint for a network without inhibition. As the membrane potential has by definition an upper limit of 1, the limit for to infinity during the phase transition (Phase II) is: (7) The calcium concentration has no upper limit and will theoretical rise to infinity (8) As the system remains only for a finite time in this second stage, will, however, remain finite. The mean membrane potential on the other hand reaches in the simulations indeed a value of 1 while the calcium concentration approaches (Figure 4 D). If the membrane potential is close to one, neurons theoretically fire at every time step. Due to the given refractory period of 4 time steps, however, only out of 100 neurons fire on average in one time step. Without inhibition too many neurons are active during this stage and distributions cannot be reasonably assessed because one will only measure one or two “endless” avalanches (Figure 5 D). Introducing inhibition changes this behavior substantially. The mean membrane potential decreases from to (Table 2) and the avalanche distribution shows now a measurable supercritical behavior (Figure 5 E, F). For measuring this avalanche distribution both excitatory and inhibitory neurons are considered. The membrane potential for the inhibitory neurons is larger than that for the excitatory neurons. This is due to their lower density (20% inhibitory as compared to 80% excitatory neurons). As in Phase I, the network will not reach a steady state in Phase II, either. However, by contrast to the first phase where activity and connectivity is slowly growing, in the second phase, connectivity and activity is quickly getting overly strong (Figure 4 B–D, Phase I and II). Therefore, the system remains supercritical for the whole second phase until pruning is reducing connectivity to the homeostatic value in Phase III. Note, that stronger inhibition dampens the membrane potential and the firing rate considerably but does not influence the supercritical behavior of the system; (Table 2) remains essentially the same across five orders of magnitude of increased inhibition (see also Figure 5 E, F). So far we have described the three development phases for our network model showing how criticality depends on network state, where the final state suggests some kind of fixed point behavior. In the following we will assess to what degree this process is characteristic for the system. To this end, we calculate its nullclines analytically [24] and compare these results with the simulations in Figure 8. For simplicity here we treat only a purely excitatory network. To be able to solve the problem analytically we assume that the change of the connectivity between neurons and their membrane potential is slow and derivatives can, thus, be set to zero. Furthermore, on longer time scales the differences between neurons are negligible and only the behavior of the means need to be considered. As a result one can calculate the nullcline of this system (see Supporting Material Text S2), which describes a hysteresis curve (Figure 8 A): (10) and are the mean values of and over all neurons. is a sigmoidal function as an approximation for the Heaviside function, which determines when an action potential is generated (see Supporting Material Text S2). In Figure 8 A we also plot the trajectories which belongs to this system and the other (trivial) nullcline, which describes the fact that the system develops into homeostasis. At fixed point development stops. In Figure 8 B we plot the actual development of and observed in the simulations. Ideally this curve should match one of the trajectories in panel A and one can see that this is essentially the case. The main deviation arises from the fact that, due to the required simplifications, the analytical solution in panel A shows during the phase transition (Phase II) infinite growth and this cannot be achieved in the simulation. This leads to a reduction in the rising slope of panel B and to the fact that the fixed point is shifted closer to the inflexion point of the isocline. When considering axons and dendrites separately, fixed point splits into a zone of many points, which correspond to the same connectivity, and hence lie on a hyperbola in Figure 8 C (dashed line). These fixed points form an omega-limit set in phase space and are represented by the equilibrium point in the --space. The approximate path of a trajectory from panels A and B is shown in Figure 8 C by the solid white line. Above we had stated that rates and membrane potentials are in Phase III fully invariant against system parameters and initial conditions, while connectivity is influenced by the level of inhibition. To this we can now add that the actual balance between axonal supply and dendritic acceptance (location of the different fixed points) remains dependent on the initial conditions (as well as on the inhibition) and should, therefore, be the most sensitive developmental parameter, e. g. much susceptible to pharmacological interference. Furthermore, as the rate essentially follows and and inversely, we can state that the isocline in Figure 8 A will, for larger inhibition, be shifted diagonally upwards away from the origin shifting the fixed point to a higher synaptic density. The dynamic behavior shown in Figure 8 is similar to that observed in the studies of Van Ooyen and Van Pelt [24] and our results show that the three development phases (Phase I, II and III) of this system are generic and independent of the chosen simulation parameters and confirm the existence of a strong phase transition. Figure 9 shows a comparison of the different criticality states between cell culture (top) and model data (bottom) summarizing some of the observations from above. Additionally, the exponent of the avalanche distribution and the time bins are given in Table 4 for each state in model and cell cultures. In the model, at the end of the transition from Poisson to power law (Figure 5 C), little connectivity in Phase I leads to an initial state similar to that observable in dissociated cell cultures (Figure 9 A). This is followed by strongly rising synaptic density in Phase II (B). Accompanying the overshoot in network activity and connectivity, the model network passes a transient phase of supercriticality (B, bottom) as do the cell cultures (B, top). Depending on the chosen strength of inhibition, we obtain in Phase III a subcritical state for the model (C, bottom,) similar to that in cell culture data (C, top). Thereafter, still in Phase III, we have gradually reduced the inhibitory strength to, hence balancing synaptic strength for inhibition and excitation (while keeping the number of inhibitory neurons constant). This leads to a final critical state in the model (D, bottom) similar to that found in cell culture data (D, top). Thus, this predicts that that developing inhibition is an important factor for the course of criticality in developing neuronal networks. Only if inhibition in the model is lowered in Phase III again, the network becomes critical. Therefore, it is likely that overall synaptic pruning in developing networks not only affects excitatory but also inhibitory synapses [31]. Moreover, neuronal networks seem to reach firing rate homeostasis earlier than the equilibrium for maturing inhibition (compare discussion in Figure 7). Additional inter-spike interval (ISI) and cross-correlation (CC) analyzes have been performed. ISIs and CCs are very similar between cultures and model across all stages but they do not contain interesting features (like oscillations) and therefore we do not show these diagrams here to save some space. The following predictions are derived from the model: These predictions are quite specific as they do not depend on the parameter choices in the model, which is one strength of this approach. Most predictions, if not all, can be tested in a straight forward way in future experiments, albeit requiring substantial and sometimes difficult experimental work which can only be addressed in future work.
  In general the chosen abstractions in the model appear to match the data description level quite well, but the question arises to what degree this still corresponds to reality in developing networks. Most importantly, the network stages described above must be related to the morphological development of dissociated neurons and their growing connectivity in culture which determines the activity pattern at every point in time [37]–[39]. It is well known [40] that the development of connectivity in cultures follows several phases. An initial phase (Phase I) is characterized by neuritic growth, followed (Phase II) by a structural overshoot and pruning followed by a maturation phase (Phase III) which finally leads to stable mean connectivity. Slowly growing connectivity in Phase I [41] leads over to the fast building of many synapses and a strong increase in activity in Phase II [38], [42], while pruning leads to Phase III with reduced number of synapses and lower activity [37]. Thereafter, firing behavior remains unchanged for two months [38], [42]. One may conclude that when synaptic pruning ceases, connectivity becomes stable and neuronal activities turn into homeostasis. Stable connectivity means that the sum of existing synapses does not vary much in time. The topology of the network can, however, not be predicted by the model as for this purpose a more detailed model of axons and dendrites would be required. However, neuronal development towards homeostasis substantially accelerates by increasing neuronal activities due to disinhibition by picrotoxin, a GABAergic synapse blocker [43]. Considering other transmitter studies, neuronal activity via increased glutamate release is likely to promote axonal outgrowth [44], [45] and therefore leads to a faster synapse formation and to an earlier maturation of the cell culture. Importantly, the behavior of dissociated neurons forming networks spontaneously occurs in any cell culture regardless of the original source of the plated neurons like cortex or hippocampus [46]. While certain simplifying assumptions had to be made to arrive at the basic differential equations (Equations 14–17) of our model, these experimental results clearly support the general dynamics assumed for our model. In our model, networks with about 20% inhibition where the only ones that reached a robust critical state. While this level of inhibition corresponds to that in real nets, the results is intriguing as homeostasis of the firing rate will also be reached with much different levels of inhibitory cells. As known from literature [47], [48], GABA changes during the development from an excitatory to an inhibitory transmitter. As this is a fast process, inhibition sets in rapidly in the overshoot Phase II [48], [49] and possibly with a too high level. As discussed above, the observed subcritical phase clearly suggests a pruning phase for the inhibition which lasts longer than the firing rate equilibration. An indication of the functional role of synaptic pruning of inhibitory synapses was recently obtained from the developing auditory system in gerbils [31]. Like others [24], [25], [50], [51], also our model assumes that the main determining force within a growing network is the attempt of the neurons to achieve on average activity homeostasis. Several existing studies indicate that neurons, which are too active, seek to reduce their firing [50], [52], whereas neurons that are too quiescent try to increase it [53], [54]. Activity reduction is achieved by a reduction of the inputs to the cell (for example dendritic withdrawal) and vice versa. At the same time, highly activated cells respond with axonal outgrowth [44], [45], [55], [56] as increased levels of intracellular calcium, as a second messenger, regulates growth cone motility and therefore affects neurite outgrowth [44], [56]–[60]. Self-organized criticality represents the situation that many systems of interconnected, nonlinear elements evolve over time into a critical state in which the probability distribution of avalanche sizes can be characterized by a power law. This process of evolution takes place without any external instructive signal. As analytically shown [7], an important feature of the power law is its scale invariance. This means that all neuronal avalanches regardless of their size (number of spikes) can be treated as physically equal [3]. Furthermore, avalanches remain stable in their spatial and temporal configuration for many hours, as already shown in cortical slices [15]. So, avalanches have optimal preconditions (equality and stability) to be a candidate for memory patterns. The stability of these effects is strongly supported by the way our model systems develop as will be discussed next. The current study shows that networks in cell cultures undergo a certain transition during their morphological development. Thus, this paper is in the tradition of a sequence of investigations [17], [40], [43], [45] that try to link cell culture activity and development to possible in vivo stages. Indications exist indeed that different activity states in cultures could be matched to in vivo states [61], but one needs to clearly state that culture and in vivo development also show clear differences. In vivo development is much more structured which will lead to differences in (ongoing) activity. As discussed above, dendritic and axonal fine structure and their spatial distribution, however, does not seem to critically affect the observed state-transitions. Hence, this supports that, at the level of avalanches, little difference might indeed exist between culture and in vivo. A study by Stewart and Plenz [21] suggests that avalanche frequency is correlated to the integrated amplitude of local field potentials, which grows until 25 DIV in their study. This indicates that also their networks had developed from a low-activity state into states that follow a power-law distribution. They show that distributions have in general an exponent of −1. 5, indicative of a branching parameter of 1 [1], and a closer look at their result suggests that transitory (sub- and supercritical) stages are also observable in this data set (see, e. g. , Fig. 3D in Stewart and Plenz 2008 [21]). A related study by Pasquale et al. [17] confirms this observation. It, thus, seems that the critical state represents the final state of the development, which – in the model – is reached together with firing rate homeostasis. This leads to a high degree of stability, which would be desirable also from a functional viewpoint. This is supported by the observation that in Phase III in the model sudden changes of the network structure (e. g. by a sudden change of inhibition) will only lead transiently to a stronger disruption of criticality. Indeed, the system soon find its way back into homeostasis and criticality is only little affected. Several previous studies [1], [17], [21], [32] focused on the exponent of the power law in the critical state. This is a characteristic parameter of the system and found to be around [1], [17], [21], [32]. We find that the exponent is in simulations and in cell cultures. Thus, the exponent matches previous results very well for the simulations. The difference in the experiments from the theoretical value of can occur from variations in the time bin, too harsh selection criteria, or a too small number of data points. Thus, deviations leading to the found value of fall into the tolerance range of these experiments. In addition, it is not clear if the theoretical value of gained from branching processes [62] can be applied to all self-organizing systems in the critical state (f. e. Bak et al. , 1987 [6]). Hence, it is equally well possible that the activity of the cultures is critical but does not exactly follow a branching process. In a previous study Beggs and Plenz [1] have shown, that the critical state is optimal for a neuronal system concerning information flow. If the system is subcritical information will die out. The opposite situation is an epileptic system with too many long avalanches (supercritical state). Thus, a neuronal network in the critical state has the maximal dynamical range to react to incoming (external) information arriving from complex interactions of the neural system with its environment. The experimental part of the current study shows that real networks will develop towards such a state and the model suggests that this state is rather stable and therefore computationally reliable. Follow-up investigations, hopefully triggered by this research, might shed a light on the structural and functional dynamics of self-organized criticality in real developing brains and possibly also contribute a better understanding of developmental pathologies.
  In order to investigate the relationship between network development and self-organized criticality, we extended the previous neurite outgrowth model by Van Ooyen and Van Pelt [10], [24], [25] by separate axons and dendrites. The model is essentially a two-dimensional recurrent neuronal network with uni-directional synapses. Model neurons are described by four equations; for activity, internal calcium concentration as well as dendritic acceptance, and axonal supply. The last two parameters determine the connectivity which is a generalisation of synaptic weights and the number of synapses between neurons. In line with previous experimental [27], [45], [50] and modelling studies [23], [24], [71], [72], the processes which determine the dynamics of this system can be summarized very briefly as: The activity of each neuron affects its calcium concentration. This, in turn, specifies the change of the dendritic and axonal offers, hence, the connectivity which will then gradually influence activity and so on. In the following we define parameters and equations. These equations are solved by the Euler method with an interval length of one simulated time step.
summary: Learning depends crucially on the synaptic distribution in a neural network. Therefore, investigating the development from which a certain distribution emerges is crucial for our understanding of network function. Morphological development is controlled by many different parameters, most importantly: neuronal activity, synapse formation, and the balance between excitation and inhibition, but it is largely unknown how these parameters interact on different time scales and how they influence the developing network structure. In our work, we consider the well-known concept of self-organized criticality. We have measured how real cell cultures change their activity patterns during the first 60 days of development traversing through different stages of criticality. With a dynamic model we can reproduce the observed developmental states and predict specific time-courses for the network parameters. For example, the model predicts a delayed, overshooting onset of inhibition with a longer time to reach maturation as compared to excitation. Furthermore, we suggest that the balance of dendrites and axons in the mature state is quite sensitive to the initial conditions of development. These and several more predictions are accessible by future experimental work and can help us to better understand neuronal networks and their parameters during development and also in the mature state.
keywords: neuroscience/neurodevelopment
  neuroscience/theoretical neuroscience
  neuroscience/neural homeostasis
title: Self-Organized Criticality in Developing Neuronal Networks