PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 December; 6(12): e1001013.
Published online 2010 December 2. doi:  10.1371/journal.pcbi.1001013
PMCID: PMC2996321

Self-Organized Criticality in Developing Neuronal Networks

Karl J. Friston, Editor

Abstract

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.

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

Introduction

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

Results

Experimental approach

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

Figure 1
Raster plots at different temporal resolutions for experimental and model data.

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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e001.jpg, which gives the quality of fit between ideal power law and actual distribution. For a system in a supercritical state An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e002.jpg is larger and for a subcritical state smaller than zero (see Methods ). Values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e003.jpg are also shown in the different panels of Figure 3. For the cultures, we receive at (on average) 19 DIV values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e004.jpg in the interval from An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e005.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e006.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e007.jpg). In general we find that the differences between all states are significant for the measured values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e008.jpg (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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e009.jpg for these four states are An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e010.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e011.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e012.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e013.jpg (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.

Figure 2
Development of the deviation from a power law An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e014.jpg of cell cultures.
Figure 3
Avalanche distribution changes during morphological development of dissociated cell cultures.

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.

Modelling approach

Neurite growth and retraction towards firing rate homeostasis

The model uses two opposing mechanisms of axonal and dendritic growth and is driven by the goal to reach homeostasis of the mean firing rate. The first mechanism regulates dendritic growth probabilities inversely to neuronal activity and the second is the axonal outgrowth promoted by activity. Specific choices for the model are being discussed in the Discussion section, where we also summarize the different specific predictions made by the model and described in detail in the next sections.

As will be shown below, the model is capable of reproducing all different patterns of neuronal activity (Figure 3) based on the implemented rules for activity-dependent structural network formation. A neuron is represented by its membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e018.jpg and its inner calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e019.jpg (see Methods ) at the time point An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e020.jpg. After a disturbance, these variables will decay in time to the resting values An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e021.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e022.jpg for the membrane potential and the calcium concentration, respectively. Every time a neuron generates an action potential (see Equation 15 in Methods ), its calcium concentration increases by a constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e023.jpg.

Dependent on the difference between the current calcium concentration and a desired homeostatic value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e024.jpg, the neuron changes its input (dendritic acceptance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e025.jpg) and output (axonal supplies An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e026.jpg) by ways of a simulated growth or withdrawal process (see Methods ). The intersection between input and output of two neurons An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e027.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e028.jpg determines the synaptic density An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e029.jpg, and hence the connectivity, between them.

The difference between an inhibitory and excitatory neuron is defined by constants An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e030.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e031.jpg, which are prefactors of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e032.jpg.

In summary, the model comprises a negative feedback loop of the following kind (Figure 4 A): Neuronal activity (1) determines the calcium level (2) in the cell. This level leads to the simulated growth pattern of the neuron(s). The growth pattern determines the effective amount of axonal supplies and dendritic acceptances (3). Thus, growth of many neurons, influencing their respective neuritic offers, will lead to different synaptic densities (4) between neurons. We use this synaptic density as the simplest way to estimate the inputs (5) to any given cell. This input will then determine the cell's activity closing the loop at (1).

Figure 4
The development of the model shows three different phases (Phase I, II and III).

These interactions lead to the effect that the model development passes through three different morphological phases (Figure 4 B–D), which we will first describe qualitatively and in the following subsections also analyze mathematically as far as possible.

The initial supplies of the axons and acceptances of the dendrites are chosen such that no connections exist. As a consequence of the resulting too low activity the dendritic acceptance increases to build synapses and to enhance the activity in the first developmental phase I. It rises slowly and, at a certain point in time, increases explosively towards a maximum. Parallel to this increase in activity, the system undergoes a morphological transition (Phase II) until it reaches homeostasis (Phase III). As discussed later (see Discussion ), this is similar to the morphological development in such cultures (see inset in Figure 4 B). At the final stage the mean activity is equal to the homeostatic value (see Methods ) and changes of the axonal supplies and dendritic acceptances are negligible.

The three different phases in the above described development can be largely understood in an analytical way and we can also describe to what degree the system approaches criticality. The difficult recurrent processes, which drive the interactions within a network and lead to a specific avalanche distribution, however, defy analytical analysis and can only be obtained from simulations. Additionally, the effects of inhibition on the network dynamic in the different developmental phases are tested by simulations.

Phase I: First developmental phase An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e042.jpg

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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e043.jpg, and the mean calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e044.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e045.jpg and average calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e046.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e047.jpg, which determines the increase of the membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e048.jpg according to the activity of the connected neurons An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e049.jpg in Equation 14 (see Methods ), is simplified to a product of the mean membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e050.jpg and an monotonous function dependent on the mean synaptic density An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e051.jpg (see below) and we get for the activity change:

equation image
(1)

The differential equation of the calcium concentration (Equation 15 in Methods ) can be written as:

equation image
(2)

With these equations, we can now consider three different degrees of synaptic densities in the first phase An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e054.jpg; namely zero, small, and medium densities and for Phase II An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e055.jpg with a large density.

Network development before synapse formation An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e056.jpg

For the initial conditions of the model without connectivity, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e057.jpg is set to zero. Therefore, from Equation 1 one can obtain that the mean activity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e058.jpg reaches the resting potential:

equation image
(3)

If this solution is entered in Equation 2, we get:

equation image
(4)

Thus, also the mean calcium concentration reaches a constant value dependent on An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e061.jpg.

Taking the limit An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e062.jpg corresponds to letting the system under the given condition An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e063.jpg relax into its end state. Note however, that the actually ongoing development (Figure 4 A) will curtail this condition as eventually An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e064.jpg.

From Figure 5 A we can see that the avalanche distribution shows a poissonian form. This is also reflected by a large negative value for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e065.jpg (Table 1, first row). This changes as soon as the model begins to make the first connections between neurons as shown in the following.

Figure 5
Avalanche distribution of the model in Phase I and II.
Table 1
The mean synaptic density An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e071.jpg influences the membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e072.jpg, avalanche distribution, and mean firing rate per time step An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e073.jpg.

Network development with small and medium connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e094.jpg

As soon as the system has reached small connectivity, the behavior of the membrane potential, calcium concentration, and avalanche distribution changes. This corresponds to a situation where we have An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e095.jpg larger than zero but smaller than An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e096.jpg. So, the system is still in Phase I. It is easy to see that the dynamics change again if the density function becomes larger than An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e097.jpg and this is later discussed in Phase II.

We can solve the differential equations (Equation 14 and 15) for the mean variables with standard methods and get:

equation image
(5)

equation image
(6)

As the synaptic density function An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e100.jpg is in this case smaller than An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e101.jpg, the product An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e102.jpg is between zero and one (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e103.jpg). Therefore, the solutions for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e104.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e105.jpg with connectivity (Equations 5 and 6) are larger than without (Equations 3 and 4) but remain bounded. As a consequence, the membrane potential and the rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e106.jpg in the system rise slowly in time (Figure 4 B–D, Phase I) dependent on the density which still growths as homeostasis is not yet reached (Figure 4 B).

Also the avalanche distribution changes slowly with rising connectivity and activity from a Poisson to a power law distribution (see transition in Figure 5 A–C, Table 1).

In the whole Phase I, the network never attains steady state. Hence connectivity and activity continue to change. Criticality essentially follows these changes. The transition from small to medium synaptic density only leads to a qualitative change in the distribution (Figure 5 B,C), which now becomes very similar to the ones measured around 18 DIV in the real cell cultures (Figure 3 A).

Phase II: Second developmental phase An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e107.jpg

Phase II of the network development is characterized by an overshoot in network activity. The membrane potential and calcium concentration (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e108.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e109.jpg) 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e110.jpg, thereby making them inhibitory (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e111.jpg is changed to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e112.jpg in this second phase). To determine the influence of different degrees of inhibition, the ratio of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e113.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e114.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e115.jpg for a network without inhibition. As the membrane potential has by definition an upper limit of 1, the limit for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e116.jpg to infinity during the phase transition (Phase II) is:

equation image
(7)

The calcium concentration has no upper limit and will theoretical rise to infinity

equation image
(8)

As the system remains only for a finite time in this second stage, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e119.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e120.jpg (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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e121.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e122.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e123.jpg (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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e124.jpg is larger than that for the excitatory neurons An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e125.jpg. This is due to their lower density (20% inhibitory as compared to 80% excitatory neurons).

Table 2
The system is overly active without inhibition so that it is not possible to determine the avalanche statistics (there exist one or two large avalanches across the whole second phase).

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; An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e149.jpg (Table 2) remains essentially the same across five orders of magnitude of increased inhibition (see also Figure 5 E,F).

Phase III: Third developmental phase

Firing rates become independent from parameter settings

Phase III is that of morphological homeostasis of the network and the network has now equilibrated reaching a steady state, where firing rate is stable in the mean.

It is obvious that the average steady state rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e150.jpg (the asterisk An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e151.jpg indicates steady state values) follows the averages of potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e152.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e153.jpg) and synaptic density An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e154.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e155.jpg), while it is inversely related to inhibition An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e156.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e157.jpg).

Let us first consider the system without inhibition. Also in this case in Phase III we receive a stable rate with constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e158.jpg. As a consequence An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e159.jpg should be constant, too. The top row for each fixed point (FP) 1–3 in Table 3 demonstrates that this is indeed the case. (The meaning of the different fixed points will be discussed in the next section. This can for now still be ignored.)

Table 3
In the homeostatic state (Phase III), the membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e160.jpg is independent of the attained fixed point (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e161.jpg) and the inhibition (ratio of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e162.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e163.jpg).

With different levels of inhibition the steady state connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e236.jpg changes. Larger inhibition leads to larger connectivity and vice versa. This is due to the effect that inhibition tries to lower the rate. As a consequence of the system being homeostatic (Equations 14–17) connectivity will increase to keep the rate constant. Because of the constant rate and the co-variation of inhibition and connectivity, we expect again that the membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e237.jpg should be constant. Table 3 shows this, too. For each ratio of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e238.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e239.jpg the membrane potential and number of spikes (firing rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e240.jpg) remain the same.

As a central conclusion we observe that rates An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e241.jpg and membrane potentials An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e242.jpg are in Phase III fully invariant against system parameters and initial conditions. Connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e243.jpg, however, is influenced by the level of inhibition.

Analytical approximation of the firing rate in the steady state

As the firing rate is the most accessible variable in cell cultures, we are now showing how to compute the firing rate in the model analytically. When the network is in a homeostatic equilibrium, the calcium concentration for each neuron on average equals the target value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e244.jpg. With this, and assuming that action potentials are uniformly distributed in time (see Supporting Information Text S1), it is possible to calculate the firing rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e245.jpg:

equation image
(9)

This solution quite accurately approximates the values for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e247.jpg obtained by the simulation (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e248.jpg see Table 3). A more detailed analysis shows that the remaining small difference arises from the discrete sampling in the numerics (not shown).

Homeostasis criticality is influenced by inhibition

Above we observed that inhibition influences the final connectivity that gives rise to network homeostasis. Here we find that also the avalanche distribution is dependent on inhibtion (Figure 6). Without inhibitory neurons the distribution is slightly supercritical. With 20% inhibitory neurons with the same synaptic weighting as the excitatory neurons (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e249.jpg), we obtain a critical distribution. Further increase of the inhibition to a ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e250.jpg of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e251.jpg drives the system significantly into a subcritical regime. A further increase to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e252.jpg does not significantly increase subcriticality anymore.

Figure 6
In the homeostatic equilibrium (Phase III), the degree of inhibition determines whether the network finally reaches a critical state or remains sub- or supercritical.

This demonstrates that criticality after equilibration, hence on the long run, depends on connectivity but neither on the mean membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e261.jpg nor on the resulting average firing rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e262.jpg.

Criticality is subject to acute changes in inhibition

This can be nicely demonstrated by disturbing an equilibrated system by a sudden change of inhibition (Figure 7). After such a jump, connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e263.jpg (and, thus, criticality, see Figure 7) changes, but mean membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e264.jpg and rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e265.jpg will relax back to their previous values. This long term process is initiated by the activity change that follows the artificially induced change of inhibition. Panels C and D in Figure 7 show that immediately after the jump, criticality changes to relatively high(low) values for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e266.jpg for the super(sub)-critical case (left vs right columns in Figure 7). A similar experiment has been performed by Beggs and Plenz [1] (see inset in Figure 7 A) where reduced inhibition also led to a supercritical system. While in our system activity fully builds back, super(sub)-criticality does not.

Figure 7
A sudden change of the inhibition in Phase III destabilizes the system.

A comparison between panel B in Figure 6, which represents the fully relaxed case, with panels C and D in Figure 7, which represent the situation immediately after the jump, shows this clearly. Hence, while the activity change leads to an immediate change in criticality, it is the lasting change of connectivity that leads to the fact that also the changed criticality persists albeit on a reduced level.

Thus, the model predicts that sudden activity changes should affect criticality in Phase III, but in a reversible way. Lasting changes of inhibition, on the other hand, should also lead to lasting small changes in the criticality without affecting the mean firing rate in the network.

Dynamic network behaviour: Isoclines and fixed points

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.

Figure 8
Development of the network in phase space.

To be able to solve the problem analytically we assume that the change of the connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e280.jpg between neurons and their membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e281.jpg 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):

equation image
(10)

An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e283.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e284.jpg are the mean values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e285.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e286.jpg over all neurons. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e287.jpg is a sigmoidal function as an approximation for the Heaviside function An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e288.jpg, 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e289.jpg, which describes the fact that the system develops into homeostasis. At fixed point An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e290.jpg development stops. In Figure 8 B we plot the actual development of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e291.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e292.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e293.jpg is shifted closer to the inflexion point of the isocline.

When considering axons and dendrites separately, fixed point An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e294.jpg splits into a zone of many points, which correspond to the same connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e295.jpg, 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e296.jpg in the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e297.jpg-An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e298.jpg-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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e299.jpg and membrane potentials An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e300.jpg are in Phase III fully invariant against system parameters and initial conditions, while connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e301.jpg is influenced by the level of inhibition. To this we can now add that the actual balance between axonal supply An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e302.jpg and dendritic acceptance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e303.jpg (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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e304.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e305.jpg and inversely An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e306.jpg, 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.

Comparison between cell culture and model development

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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e307.jpg 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, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e308.jpg) similar to that in cell culture data (C, top). Thereafter, still in Phase III, we have gradually reduced the inhibitory strength to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e309.jpg, 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).

Figure 9
Comparing real data (top) with model (bottom).
Table 4
The exponent of the avalanche distributions and the time bins An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e310.jpg are given for the model (top) and cell culture (bottom) data in the different states.

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.

Predictions of the model

The following predictions are derived from the model:

  • Criticality at the end of development is optimally reached with 20% inhibition with a strength equal to that of the excitation. This observation does not depend crucially on the distribution of the inhibitory neurons and corresponds approximately to the normal degree of inhibition in cortical networks and cultured networks, respectively. A higher degree leads to sub- and a lower degree to supercritical behavior. A subcritical state is observed in cell cultures before they reach a mature state. This predicts that the onset of inhibition in the cultures must be strong with a connectivity much larger than in the end. Also the time-course of reaching firing rate homeostasis appears to be shorter than the curbing of this overly strong inhibition which takes longer and thus leads to the subcritical state observed in the interim. These are model predictions, because there is currently no data existing about the temporal development of inhibition. This data would be required to extend our model by implementing some time-course (dynamic coupling function) of the inhibitory development. Due to the lack of this data, this seems not useful at the moment because there is no way to constrain such a model extension. In general, the average homeostatic firing rate is independent of the level of inhibition and will in Phase III be reached regardless. All these predictions could be tested by measuring the degree of inhibition in the developing cultures and by pharmacologically interfering with the normal development forcing cultures to develop stronger (or weaker) inhibitory networks.
  • Average rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e329.jpg and membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e330.jpg are at the end of the development fully invariant against system parameters and initial conditions with which development had started. Connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e331.jpg within the network, however, is influenced by the level of inhibition. Strong inhibition leads to more and weak inhibition to fewer connections. The latter could be assessed in parallel with the first prediction performing histological analyzes.
  • For a network that has reached homeostasis, criticality can probably still be influenced by a sudden, prolonged change of inhibition. Following the second prediction, we expect some lasting connectivity changes to take place in these cases leading to a mildly changed criticality. Remarkably, first experimental evidence exists that an acute change in inhibition in fact alters criticality in mature cultured networks [1]. It further remains to be tested whether the mean firing rate will not be affected by a change in inhibition. It should quickly relax to its old value. Another recent experiment from Shew et al. (2009) [32] shows that AP5-DNQX, which blocks excitation, acts – at a first glance – like increasing inhibition in our model. However, it is not clear whether this experimental result and the model predictions can be compared at this stage, because theory suggests that inhibition acts dissipating in branching processes [33]. Thus, with respect to criticality decrease of excitation does not necessarily correspond to increase in inhibition. Nevertheless, it would be of great interest to assess whether this prediction holds and especially how long a homeostatic culture will still remain susceptible to such interference. Additional long term mechanisms, not captured by our model, might in reality terminate this effect after some time.
  • The model further predicts that the actual balance between axonal supply An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e332.jpg and dendritic acceptance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e333.jpg is quite sensitively depending on the initial conditions under which a cell culture starts its individual development. Thus, any histological analysis of connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e334.jpg should best be performed by carefully assessing dendritic and axonal parameters in the different cultures. Even with very similar initial conditions, we expect those to vary widely across cultures, while the total connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e335.jpg should be very similar. In this study this is reflected in the behavior of the PKC (Protein Kinase C)-inhibited cultures, which do not show any visible differences in their avalanche and firing rate characteristics as compared to non-treated cultures. The PKC-inhibited cultures express much richer dendritic structures [34][36]. However, the lack of difference in activity patterns found here suggests that the final synaptic density in these cultures does not substantially differ from that in the controls. The model predicts that this should go hand in hand with a shift of the fixed point along the hyperbola in Figure 8 C, which, however, does not alter the found avalanches. The prediction of a fix-point shift could be tested by detailed and complex histological analyzes of the existing synapses in controls and PKC-inhibited cultures, which goes beyond the scope of this study.

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.

Discussion

In the current study we have investigated how the activity patterns in developing cell cultures can be measured and modeled in terms of self-organized criticality. We have shown that the activity distributions in real cultures undergo a transition from a stage with little activity to a supercritical and then a subcritical state and finally to critical behavior. These transitions were significant for the cell cultures analyzed.

We used an extended version of the neurite outgrowth model by Van Ooyen and co-workers [10], [24], [25] with separate axonal and dendritic fields. The axonal and dendritic growth is driven by the goal to reach firing rate homeostasis as modeled in previous papers by Butz and co-workers [26], [27]. The model was able to reproduce the different developmental phases and several interesting predictions have been made.

Relating criticality to morphological development

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 in neuronal networks

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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e336.jpg [1], [17], [21], [32]. We find that the exponent is An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e337.jpg in simulations and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e338.jpg in cell cultures. Thus, the exponent matches previous results very well for the simulations. The difference in the experiments from the theoretical value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e339.jpg 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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e340.jpg fall into the tolerance range of these experiments. In addition, it is not clear if the theoretical value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e341.jpg 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.

Methods

Experimental approach and data evaluation

Preparation of the cell cultures

Primary cortical cell cultures were prepared as described previously [63], [64]. Cells were derived from cortices of neonatal wistar rats by mechanical (chopping with scalpel, trituration) and enzymatical (0.05%, Trypsin, 15min at An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e342.jpgC.) dissociation and plated at densities (CASY cell counter, Innovatis) of 500,000 cells per An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e343.jpg onto polyethyleneimine-coated micro-electrode arrays (59 TiN electrodes, 200/500An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e344.jpgm electrode pitch, Multi Channel Systems). Cultures developed in 1ml growth medium, minimum essential medium (Gibco) supplemented with heat-inactivated horse serum (5%), L-glutamine (0.5mM), glucose (20mM) and gentamycin (10An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e345.jpgg/ml). One third of medium was exchanged twice per week. Cultures were maintained at 5% An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e346.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e347.jpgC. In a subset of cultures PKC (Protein Kinase C) was chronically inhibited by addition of a PKC antagonist (Gö6976, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e348.jpg, Calbiochem) at the first exchange of the culture medium at An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e349.jpg DIV. This different treatment has no significant influence on different parameters of the cell cultures (see Table 5) and, therefore, on the results of this paper.

Table 5
Cell cultures with PKC and without PKC (untreated) are compared.

Electrophysiology

Electrophysiological recordings were performed on the different DIVs at the same time for one hour under culture conditions with a MEA1060-BC system amplifier (Multi Channel Systems) [65]. Raw electrode signals were digitally high-pass filtered at 200Hz and action potentials were detected by voltage threshold (3 times of standard deviation from the mean) using MC-Rack software (Multi Channel Systems).

Selection criteria

Clustering of neurons [66] at few electrodes can distort the avalanche statistics as clustering is a culture phenomenon and not seen in-vivo. To avoid clustering induced effects, we demand that activity has to be nearly uniformly distributed across all electrodes. So, all electrodes are excluded from the statistics which have measured more activity than two times the standard deviation from the mean activity per electrode. As shown in the literature [1], [3], [17] a too small number of measuring electrodes distorts the avalanche distribution (fewer long avalanches) and decreases classification of different states. This is avoided by choosing only samples with at least 50 active electrodes.

A total of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e366.jpg cultures has been originally considered in this study. Of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e367.jpg cultures An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e368.jpg were controls and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e369.jpg were PKC inhibited. A total of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e370.jpg cultures did not obey the tight prior selection criteria for allowing rigorous criticality analysis and had to be excluded. The final number of analyzed cultures was, thus, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e371.jpg controls and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e372.jpg PKC inhibited. This shows that both groups were equally affected by the selection criteria. Interestingly, and also predicted by the model (see “Predictions of the Model”), no differences with respect to the activity analyzes of the current study were found between controls and PKC-inhibited cultures. Thus results were pooled.

Definition of avalanches

order to assess the distribution of avalanches in cell cultures and in the model in the same manner, we search for the beginning and the end of an avalanche by a gliding time bin of a fixed size. Whenever the system is silent (no spikes) for at least the duration of the time bin, an avalanche ends and a new one starts with the next spike. The time bin is the mean time interval between two spikes in the system. Too long time intervals are sorted out by first calculating the mean cross-correlation of all electrode signals and secondly by getting the time value for which 99% of the integration area is under this mean cross-correlation curve. This time value is the maximum time interval which is taken into account of the mean time interval [1]. This way we ensure that bursts on a longer timescales do not distort the statistics. This definition of avalanches can be used for cell culture and model data.

Measuring the deviation from a power law

To distinguish between the different states of an SOC system, the measure An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e373.jpg is defined. For this, the theoretical power law distribution An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e374.jpg is calculated by a linear regression of the linear start on the left side (at approximately An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e375.jpg) of the distribution An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e376.jpg in the log-log-plot with MATLAB to the end of the linear behaviour. Now, for each data point An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e377.jpg the regression An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e378.jpg is subtracted from An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e379.jpg. The mean of the differences of all data points is the measure An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e380.jpg.

equation image
(11)

We define the following relation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e382.jpg and the state of the system:

equation image

equation image

equation image

These thresholds are heuristic as there is no theoretical background for this. In general they correspond well to state classification if made by human inspection. Note, the results of this paper are not crucial dependent on narrow threshold margins.

Additional tests for criticality and data evaluation

Measuring a power law for the avalanche distribution is not sufficient to conclusively show that a system is in the critical state, because a power law can also result from the summation of two exponentials [7]. Therefore, the critical state in model and cell cultures has to be analysed by additional tests showing the scale-free behavior. Several tests were performed and results are shown here. First, the avalanche distribution has to show a power law even with less neurons (model) or electrodes (cell culture), indicating that the system is spatially scale-free. Figure 10 A,B demonstrates this. Furthermore, a system in the critical state has also to be temporally scale-free. To show this, different time bins are used for analyzing the avalanche distribution. Also these distributions show a power law relation (Figure 10 C,D) for model and cell cultures.

Figure 10
Additional tests for criticality used for cell cultures and model.

A third informative test is to assess the scale-free behavior of the inter-avalanche intervals (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e390.jpg). For this a minimum event size An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e391.jpg has to be introduced (minimum number of spikes in one avalanche) and then the time interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e392.jpg between the occurrence of two avalanches larger as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e393.jpg can be measured. Thus, for all values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e394.jpg the probability distribution An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e395.jpg of getting an inter-avalanche interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e396.jpg given An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e397.jpg is assessed and can be re-scaled by the rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e398.jpg of having an avalanche larger than An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e399.jpg per time unit (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e400.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e401.jpg). If this is done for different An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e402.jpg and all distributions form a single curve, the system is scale-free and the curve is the scaling function An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e403.jpg [67].

equation image
(12)

This is done for our model and cell cultures in the critical state (Figure 10 E,F). Note, the actual shape of the different scaling functions is of less importance. The scale-free property is confirmed as long as all functions for a given system collapse onto the same function [67]. This is indeed observed in Figure 10 E,F.

Finally, a fourth test for criticality is the Fano Factor [68][70], for which the number of spikes An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e405.jpg in a time window from An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e406.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e407.jpg has to be considered using following equation

equation image
(13)

The Fano Factor An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e409.jpg assumes a point process of events (spikes) and relates the clustering of these events to a Poisson process for which An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e410.jpg. When the Fano Factor is below one, it indicates that the point process is more orderly than a Poisson process, and a Fano Factor above one indicates increased clustering at the given time scale An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e411.jpg [69]. For a scale-free point process (e.g.; a system in the critical state), the Fano Factor needs to be a power law with the form An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e412.jpg. The exponent is an approximation of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e413.jpg exponent which is related to the exponent of the avalanche distribution [9].

For the critical state in our model and cell cultures the Fano Factor shows a power law behavior for a wide range of time windows An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e414.jpg (Figure 10 G,H). There are no large differences between model and cell culture exponents (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e415.jpg for model and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e416.jpg for cell cultures). Only at large values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e417.jpg, around An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e418.jpg steps or An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e419.jpg, model and cell culture data do not show a power law behavior anymore and start to differ. However, the important range for the avalanche analysis is at smaller An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e420.jpg-values (compare time bin An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e421.jpg in Figure 10 G,H and Table 4).

Computational modelling approach

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 An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e422.jpg, internal calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e423.jpg as well as dendritic acceptance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e424.jpg, and axonal supply An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e425.jpg. The last two parameters determine the connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e426.jpg 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.

Membrane potential

As in the main text in the following, mean values are given as upper case letters, while lower case letters indicate individual values. For a fixed connectivity, given by the synaptic density between dendritic and axonal offers, each neuron has a certain activity. In accordance to the definition of the neuron model in the work by Abbott and Rohrkemper [10], the activity of the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e427.jpg-th neuron at the time point An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e428.jpg is given by a membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e429.jpg, limited by a hard bound to 1, which decays in time exponentially with time constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e430.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e431.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e432.jpg is the resting membrane potential. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e433.jpg increases proportionally to the connectivity An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e434.jpg if a neighboring neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e435.jpg generates an action potential (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e436.jpg; An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e437.jpg is a uniformly distributed number between 0 and 1. This relation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e438.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e439.jpg is obtained by the Heaviside-function An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e440.jpg).

equation image
(14)

An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e442.jpg defines if a presynaptic neuron is inhibitory (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e443.jpg) or excitatory (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e444.jpg). In the beginning of the simulation, all neurons are excitatory comparable to the very early development of biological neuronal networks [48]. At some point during simulation, a certain subset of neurons (20% of all) is converted into inhibitory neurons (see subsection Phase II of the Results section). We further define a refractory period of four time steps.

Calcium concentration

We model the calcium dynamics in our neuron model related to the work by Abbott and Rohrkemper [10]. The membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e445.jpg affects the calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e446.jpg which has a slower exponential time constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e447.jpg. If a neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e448.jpg is active, it receives an influx of calcium and the concentration increases by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e449.jpg.

equation image
(15)

An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e451.jpg determines the change of the synaptic density An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e452.jpg.

Dendritic acceptance, axonal supply and connectivity

The development of dendrites and axons depends indirectly, by ways of the calcium concentration, on the activity. Lipton and Kater [45] showed that the deviation of the calcium concentration of a neuron from a certain target value determines its outgrowth. In line with previous modelling studies [10], [23], [24], [72], we defined a homeostatic value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e453.jpg for which the axonal and dendritic offers (and with them the synaptic density) remains unchanged if An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e454.jpg equals An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e455.jpg. The rules for growth and retraction of continuous axonal supply and dendritic acceptance were taken from the modelling approach by Dammasch and Butz [23], [73] and can be described as follows: If the neuron has a too high membrane potential An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e456.jpg, hence a too high calcium concentration An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e457.jpg, the dendritic acceptance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e458.jpg shrinks proportionally to a constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e459.jpg leading to a decrease of its input. If the membrane potential is too low, the dendritic acceptance increases to receive more inputs raising the calcium concentration by an increasing membrane potential to the homeostatic value. The axonal supply An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e460.jpg behaves equally but with inverted sign and different growth rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e461.jpg.

equation image
(16)

equation image
(17)

Finally, we define the connectivity between neurons An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e464.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e465.jpg by:

equation image
(18)

with

equation image

equation image

and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e469.jpg as distance between neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e470.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001013.e471.jpg. This essentially represents the overlap of the axonal and dendritic probability zones which can be understood as an abstract representation for the probability of synapse formation. In the simulations the neurons are distributed on a open grid to avoid developmental unsteadinesses followed from irregularities in the neuronal density. However, the overall behavior is not effected by this.

Parameter settings

The following table ( Table 6 ) shows all standard parameter values. Exceptions are indicated in the text.

Table 6
Standard parameters used for the simulations.

Supporting Information

Text S1

Derivation of the steady state firing rate R*.

(0.03 MB PDF)

Text S2

Derivation of the nullcline of the model.

(0.03 MB PDF)

Acknowledgments

We thank Anna Levina for helpful discussions about avalanche statistics.

Footnotes

The authors have declared that no competing interests exist.

The project was supported by BCCN grants from the German ministry for research and education (BMBF) via the Bernstein Center for Computational Neuroscience (BCCN) Göttingen as well as Freiburg. F.W. acknowledges BMBF BFNT funding (Project 3a), and M.B. from the Computational Life Sciences programme grant (635.100.017) of the Netherlands Organization for Research (NWO). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Beggs J, Plenz D. Neuronal avalanches in neocortical circuits. J Neurosci. 2003;23(35):11167–11177. [PubMed]
2. Chialvo D. Critical brain networks. Physica A. 2004;340:756–765.
3. Plenz D, Thiagarajan T. The organizing principles of neuronal avalanches: Cell assemblies in the cortex. Trends Neurosci. 2007;30:101–110. [PubMed]
4. Gutenberg B, Richter C. Seismicity of the earth. Princeton: Princeton University Press; 1956. 273
5. Harris T. The theory of branching processes. New York: Dover Publications; 1989. 230
6. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of 1/f noise. Phys Rev Lett. 1987;59:381–384. [PubMed]
7. Newman M. Power laws, pareto distributions and zipf's law. Contemp Phys. 2005;46:323–351.
8. Paczuski M, Maslov S, Bak P. Avalanche dynamics in evolution, growth and depinning models. Phys Rev E. 1996;53:414–443. [PubMed]
9. Bak P, Tang C, Wiesenfeld K. Self-organized criticality. Phys Rev A. 1988;38:364–374. [PubMed]
10. Abbott L, Rohrkemper R. A single growth model constructs critical avalanche networks. Prog Brain Res. 2007;165:9–13. [PubMed]
11. Chen DM, Wu S, Guo A, Yang Z. Self-organized criticality in a cellular automaton model of pulse-coupled integrate-and-fire neurons. J Physiol A Math Gen. 1995;28:5177–5182.
12. Corral A, Perez C, Diaz-Guilera A, A A. Self-organized criticality and synchronization in a lattice model of integrate-and-fire oscillators. Phys Rev Lett. 1995;74:118–121. [PubMed]
13. Eurich C, Herrmann J, Ernst U. Finite-size effects of avalanche dynamics. Phys Rev E. 2002;66:066137. [PubMed]
14. Stassinopoulos D, Bak P. Democratic reinforcement: A principle for brain function. Phys Rev E. 1995;51:5033–5039. [PubMed]
15. Beggs J, Plenz D. Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci. 2004;24:5216–5229. [PubMed]
16. Mazzoni A, Broccard F, Garcia-Perez E, Bonifazi P, Ruaro M, et al. On the dynamics of the spontaneous activity in neuronal networks. PLoS ONE. 2007;2(5):e439. doi: 10.1371/journal.pone.0000439. [PMC free article] [PubMed]
17. Pasquale V, Massobrio P, Bologna L, Chiappalone M, Martinoia S. Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience. 2008 doi: 10.1016/j.neuroscience.2008.03.050. [PubMed]
18. Levina A, Herrmann J, Geisel T. Dynamical synapses causing self-organized criticality in neural networks. Nat Phys. 2007;3:857–860.
19. Bornholdt S, Röhl T. Self-organized critical neural networks. Phys Rev E. 2003;67:066118. [PubMed]
20. Meisel C, Gross T. 2009. Self-organized criticality in a realtistic model of adaptive neural networks arXiv: 0903.2987v1.
21. Stewart C, Plenz D. Homeostasis of neuronal avalanches during postnatal cortex development in vitro. J Neurosci Methods. 2008;169:405–416. [PMC free article] [PubMed]
22. Gürel T, Egert U, Kandler S, De Raedt L, Rotter S. Predicting spike activity in neuronal cultures. 2007. pp. 2942–2947. In: IEEE International Joint Conference on Neural Networks.
23. Dammasch IE, Wagner GP, Wolff JR. Self-stabilization of neuronal networks. i. the compensation algorithm for synaptogenesis. Biol Cybern. 1986;54:211–222. [PubMed]
24. Van Ooyen A, Van Pelt J. Activity-dependent outgrowth of neurons and overshoot phenomena in developing neural networks. J Theor Biol. 1994;167:27–43.
25. Van Ooyen A, Van Pelt J, Corner M. Implications of activity-dependent neurite outgrowth for neuronal morphology and network development. J Theor Biol. 1995;172:63–82. [PubMed]
26. Butz M, Teuchert-Noodt G, Grafen K, Van Ooyen A. Inverse relationship between adult hippocampal cell proliferation and synaptic rewiring in the dentate gyrus. Hippocampus. 2008;18(9):879–898. [PubMed]
27. Butz M, Van Ooyen A, Wörgötter F. A model for cortical rewiring following deafferentation and focal stroke. Front Comp Neurosci. 2009;3:10. [PMC free article] [PubMed]
28. Helias M, Rotter S, Gewaltig MO, Diesmann M. Structural plasticity controlled by calcium based correlation detection. Front Comp Neurosci. 2008;10.3389:1–21. [PMC free article] [PubMed]
29. Pospischil M, Piwkowska Z, Bal T, Destexhe A. Characterizing neuronal activity by describing the membrane potential as a stochastic process. J Physiol Paris. 2009;1–2:98–106. [PubMed]
30. Chen B, Boukamel K, Kao J, Roerig B. Spatial distribution of inhibitory synaptic connections during development of ferret primary visual cortex. Exp Brain Res. 2005;160:496–509. [PubMed]
31. Magnusson AK, Kapfer C, Grothe B, Koch U. Maturation of glycinergic inhibition in the gerbil medial superior olive after hearing onset. J Physiol (Lond) 2005;568:497–512. [PubMed]
32. Shew WL, Yang H, Petermann T, Roy R, Plenz D. Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J Neurosci. 2009;29:15595–15600. [PMC free article] [PubMed]
33. Lauritsen K, Zapperi S, Stanley H. Self-organized branching process: Avalanche models with dissipation. Phys Rev E. 1996;54:2483–2488. [PubMed]
34. Egert U, Schrenk K, Kapfhammer J, Metzger F. Electrophysiological characterization of long and short-term drug effects in acute slices and organotypic cultures of the cerebellum. In: Eur J Neurosci. 2000;12:141.
35. Metzger F, Kapfhammer J. Protein kinase c activity modulates dendritic differentiation of rat purkinje cells in cerebellar slice cultures. Eur J Neurosci. 2000;12:1993–2005. [PubMed]
36. Okujeni S, Kandler S, Egert U. Homeostatic regulation of activity in cortical networks. 2008. pp. 74–75. In: Proc. 6th Int. Meeting on Substrate-Integrated Microelectrodes.
37. Van Huizen F, Romijn H, Habets A. Synaptogenesis in rat cerebral cortex cultures is affected during chronic blockade of spontaneous bioelectric activity by tetrodotoxin. Brain Res. 1985;19:67–80. [PubMed]
38. Habets AM, Van Dongen AM, Van Huizen F, Corner MA. Spontaneous neuronal firing patterns in fetal rat cortical networks during development in vitro: a quantitative analysis. Exp Brain Res. 1987;69:43–52. [PubMed]
39. Muramoto K, Ichikawa M, Kawahara M, Kobayashi K, Kuroda Y. Frequency of synchronous oscillations of neuronal activity increases during development and is correlated to the number of synapses in cultured cortical neuron networks. Neurosci Lett. 1993;163:163–165. [PubMed]
40. Corner M, Ramakers G. Spontaneous firing as an epigenetic factor in brain development: Physiological consequences of chronic tetrodotoxin and picrotoxin exposure on cultured rat neocortex neurons. Brain Res. 1992;65:57–64. [PubMed]
41. Nakanishi K, Kukita F. Functional synapses in synchronized bursting of neocortical neurons in culture. Brain Res. 1998;795:137–146. [PubMed]
42. Kamioka H, Maeda E, Jimbo Y, Robinson HP, Kawana A. Spontaneous periodic synchronized bursting during formation of mature patterns of connections in cortical cultures. Neurosci Lett. 1996;206:109–112. [PubMed]
43. Van Huizen F, Romijn H, Habets A, Van den Hooff P. Accelerated neural network formation in rat cerebral cortex cultures chronically disinhibited with picrotoxin. Experimental Neurology. 1987;97:280–288. [PubMed]
44. Mattson M. Neurotransmitters in the regulation of neuronal cytoarchitecture. Brain Res. 1988;13:179–212. [PubMed]
45. Lipton S, Kater S. Neurotransmitter regulation of neuronal outgrowth, plasticity and survival. Trends Neurosci. 1989;12:265–270. [PubMed]
46. van den Pol AN, Strecker GJ, Dudek FE. Excitatory and inhibitory amino acids and synaptic transmission in the suprachiasmatic nucleus. Prog Brain Res. 1996;111:41–56. [PubMed]
47. Ganguly K, Schinder A, Wong S, Poo MM. Gaba itself promotes the developmental switch of neuronal gabaergic responses from excitation to inhibition. Cell. 2001;105:521–532. [PubMed]
48. Jiang B, Huang Z, Morales B, Kirkwood A. Maturation of gabaergic transmission and the timing of plasticity in visual cortex. Brain Res Rev. 2005;50:126–133. [PubMed]
49. Kobayashi M, Hamada T, Kogo M, Yanagawa Y, Obata K, et al. Developmental profile of gabaa-mediated synaptic transmission in pyramidal cells of the somatosensory cortex. Eur J Neurosci. 2008;28:849–861. [PubMed]
50. Wolff J, Wagner G. Synergetics of the brain. Springer; 1983. Selforganization in synaptogenesis: interaction between the formation of excitatory and inhibitory synapses. pp. 50–59.
51. Turrigiano GG, Nelson SB. Homeostatic plasticity in the developing nervous system. Nat Rev Neurosci. 2004;5:97–107. [PubMed]
52. Pratt KG, Watt AJ, Griffith LC, Nelson SB, Turrigiano GG. Activity-dependent remodeling of presynaptic inputs by postsynaptic expression of activated camkii. Neuron. 2003;39:269–281. [PubMed]
53. Kirov SA, Harris KM. Dendrites are more spiny on mature hippocampal neurons when synapses are inactivated. Nat Neurosci. 1999;2:878–883. [PubMed]
54. Kirov SA, Goddard CA, Harris KM. Age-dependence in the homeostatic upregulation of hippocampal dendritic spine number during blocked synaptic transmission. Neuropharmacology. 2004;47:640–648. [PubMed]
55. Rekart JL, Sandoval CJ, Routtenberg A. Learning-induced axonal remodeling: evolutionary divergence and conservation of two components of the mossy fiber system within rodentia. Neurobiol Learn Mem. 2007;87:225–235. [PubMed]
56. Hutchins BI, Kalil K. Differential outgrowth of axons and their branches is regulated by localized calcium transients. J Neurosci. 2008;28:143–153. [PMC free article] [PubMed]
57. Mattson MP, Kater SB. Calcium regulation of neurite elongation and growth cone motility. J Neurosci. 1987;7:4034–4043. [PubMed]
58. Kater SB, Mattson MP, Cohan C, Connor J. Calcium regulation of the neuronal growth cone. Trends Neurosci. 1988;11:315–321. [PubMed]
59. Kater SB, Mattson MP, Guthrie PB. Calcium-induced neuronal degeneration: a normal growth cone regulating signal gone awry (?). Ann N Y Acad Sci. 1989;568:252–261. [PubMed]
60. Korkotian E, Segal M. Morphological constraints on calcium dependent glutamate receptor trafficking into individual dendritic spine. Cell Calcium. 2007;42:41–57. [PubMed]
61. Stett A, Egert U, Guenther E, Hofmann F, Meyer T, et al. Biological application of microelectrode arrays in drug discovery and basic research. Anal Bioanal Chem. 2003;377:486–495. [PubMed]
62. Zapperi S, Lauritsen K, Stanley H. Self-organized branching processes: Mean-field theory for avalanches. Phys Rev Lett. 1995;75:4071–4074. [PubMed]
63. Shahaf G, Marom S. Learning in networks of cortical neurons. J Neurosci. 2001;21(22):8782–8788. [PubMed]
64. Marom S, Shahaf G. Development, learning and memory in large random networks of cortical neurons: lessons beyond anatomy. Q Rev Biophys. 2002;35:63–87. [PubMed]
65. Egert U, Knott T, Schwarz C, Nawrot M, Brandt A, et al. Mea-tools: an open source toolbox for the analysis of multi-electrode data with matlab. J Neurosci Methods. 2002;117:33–42. [PubMed]
66. Soussou W, Yoon G, Brinton R, Berger T. Neuronal network morphology and electrophysiology of hippocampal neurons cultured on surface-treated multielectrode arrays. IEEE Trans Biomed Eng. 2007;54:1309–1320. [PubMed]
67. Corral A, Telesca L, Lasaponara R. Scaling and correlations in the dynamics of forest-fire occurence. Phys Rev E. 2008;77:016101. [PubMed]
68. Lowen SB, Teich MC. The periodogram and allan variance reveal fractal exponents greater than unity in auditory-nerve spike trains. J Acoust Soc Am. 1996;99:3585–3591. [PubMed]
69. Lowen S, Ozaki T, Kaplan E, Saleh B, Teich M. Fractal features of dark, maintained, and driven neural discharges in the cat visual system. Methods. 2001;24:377–394. [PubMed]
70. García-Marín A, Jiménez-Hornero F, Ayuso J. Applying multifractality and the self-organized criticality theory to describe the temporal rainfall regimes in Andalusia (southern Spain). Hydrol Process. 2008;22:295–308.
71. Van Ooyen A, editor. Modeling Neural Development. Cambridge: MIT Press; 2003. 336
72. Butz M, Wörgötter F, Van Ooyen A. Activity-dependent structural plasticity. Brain Res Rev. 2009;60:287–305. [PubMed]
73. Butz M, Lehmann K, Dammasch I, Teuchert-Noodt G. A theoretical network model to analyse neurogenesis and synaptogenesis in the dentate gyrus. Neural Netw. 2006;19:1490–1505. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science