# On the temporal organization of neuronal avalanches

^{1}Institute of Computational Physics for Engineering Materials, ETH, Zurich, Switzerland^{2}Departamento de Física, Universitade Federal do Ceará, Fortaleza, Brazil^{3}Section on Critical Brain Dynamics, National Institute of Mental Health, National Institute of Health, Bethesda, MD, USA^{4}Department of Industrial and Information Engineering, Second University of Naples, National Institute for Nuclear Physics Gr. Coll. Salerno, Aversa, Italy

Spontaneous activity of cortex *in vitro* and *in vivo* has been shown to organize as neuronal avalanches. Avalanches are cascades of neuronal activity that exhibit a power law in their size and duration distribution, typical features of balanced systems in a critical state. Recently it has been shown that the distribution of quiet times between consecutive avalanches in rat cortex slice cultures displays a non-monotonic behavior with a power law decay at short time scales. This behavior has been attributed to the slow alternation between up and down-states. Here we further characterize the avalanche process and investigate how the functional behavior of the quiet time distribution depends on the fine structure of avalanche sequences. By systematically removing smaller avalanches from the experimental time series we show that size and quiet times are correlated and highlight that avalanche occurrence exhibits the characteristic periodicity of θ and β/γ oscillations, which jointly emerge in most of the analyzed samples. Furthermore, our analysis indicates that smaller avalanches tend to be associated with faster β/γ oscillations, whereas larger ones are associated with slower θ and 1–2 Hz oscillations. In particular, large avalanches corresponding to θ cycles trigger cascades of smaller ones, which occur at β/γ frequency. This temporal structure follows closely the one of nested θ − β/γ oscillations. Finally we demonstrate that, because of the multiple time scales characterizing avalanche dynamics, the distributions of quiet times between avalanches larger than a certain size do not collapse onto a unique function when rescaled by the average occurrence rate. However, when considered separately in the up-state and in the down-state, these distributions are solely controlled by the respective average rate and two different unique function can be identified.

## 1. Introduction

During sleep or under anesthesia, as well as *in vitro*, ongoing or spontaneous activity in cortex alternates between active periods with high probability of action potential firing and quiescent periods characterized by sparse firing (Plenz and Aertsen, 1996; Cossart et al., 2003; Cunningham et al., 2006; Hahn et al., 2006). These extracellular spiking dynamics correspond to so-called up and down-state fluctuations in the intracellular membrane potential of cortical neurons (Steriade et al., 1993; Plenz and Kitai, 1996; Wilson, 2008). During up- states, the intracellular membrane potential is close to firing threshold allowing neurons to fire action potentials in response to synaptic input. In contrast, the membrane potential is more hyperpolarized during the down-state leading to low probability of firing. The up-state is generally considered a cortical network property that arises from the propagation of activity among recurrently connected neurons (Plenz and Kitai, 1996; McCormick et al., 2003; Wilson, 2008; Millman et al., 2010). The resulting synaptic input depolarizes neurons beyond threshold supporting and prolonging the up-state. In that context, the up-state should be considered a metastable state, i.e., the membrane potential would rapidly decay to resting value, if network mechanisms prevented the required excitability or excitatory synaptic drive for individual neurons.

Conversely, down-states reflect relatively quiescent network periods during which the membrane potential of most neurons is close to or even lower than their resting value. Down-states generally result from disfacilitation, i.e., a substantial reduction or lack of excitatory drive in the network (Cowan and Wilson, 1994; Timofeev et al., 2001). Transitions to the down-state can be caused by various mechanisms such as synaptic depression at glutamatergic synapses (Stevens and Tsujimoto, 1995; Staley et al., 1998), an increase of a factor inhibiting glutamate release, such as nucleoside adenosine (Thompson et al., 1992), blockage of receptor channels by the presence, for instance, of external magnesium (Maeda et al., 1995), or spike adaptation, which arises from the intracellular accumulation of calcium entering during the action potential and opening potassium channels (Sanchez-Vives et al., 2000). Transitions to the up-state are generally thought to arise from non-linear amplification following recovery from disfacilitation. For example, spontaneous single action potentials, spontaneous miniature synaptic release, and recovery from synaptic vesicle depletion, i.e., synaptic depression, can cooperate to a non-linear amplification of small amplitude signals leading to the generation of larger depolarizing events rapidly transitioning the network to the up- state, as observed in cortical slabs (Timofeev et al., 2000) and slice cultures (Plenz and Aertsen, 1996).

During up-states, which usually last up to several hundreds of milliseconds, cortical neurons have been shown to fire irregularly often during nested oscillations (e.g., Plenz and Kitai, 1996). This highly variable firing pattern at short time scales of just a few milliseconds, over the last decade, has been found to reflect in fact a precise, scale-invariant organization of activity, so-called neuronal avalanches (Beggs and Plenz, 2003; Mazzoni et al., 2007; Gireesh and Plenz, 2008; Pasquale et al., 2008; Petermann et al., 2009; Shriki et al., 2013). Neuronal avalanches are intermittent bursts of activity cascades whose sizes and durations follow power law statistics, a typical feature of systems at criticality (Stanley, 1971). The statistics of time intervals separating successive avalanches has been recently studied in the spontaneous activity of rat cortex slice cultures (Lombardi et al., 2012). In Lombardi et al. (2012), these intervals are called waiting times and defined as the difference between the ending and starting time of consecutive avalanches. Here and in the following we will adopt a slightly different notation (Sanchez et al., 2002): We call quiet times the time intervals between the ending and starting time of consecutive avalanches, whereas we refer to waiting times as time intervals between starting times of consecutive avalanches.

The quiet time distribution, is widely used in the stochastic analysis of natural phenomena, such as earthquakes, solar flares (de Arcangelis et al., 2006a), and rock fractures, where it is usually called waiting time distribution. Indeed, for these phenomena the waiting times, do not differ from quiet times because event durations can be neglected and processes can be consistently treated as point processes. For neuronal avalanches this approximation is not always valid since the shortest quiet times are comparable with some avalanche durations, as we will show in the following. While numerous similarities between earthquakes and neuronal avalanches have been found (Plenz, 2012), the quiet time distribution has only been incompletely analyzed so far for avalanches. Of particular interest are the universal temporal scaling features observed for earthquakes. Distribution of earthquake waiting times, in which waiting times are restricted to earthquakes above a given magnitude threshold, depend on the threshold, but nevertheless collapse onto a universal, i.e., threshold independent, function when waiting times are rescaled by the average rate (Corral, 2004). This property reveals that seismicity has a complex organization in time with universal properties: the removal of small events by increasing the minimal detection threshold does not affect the fundamental organization of earthquake occurrence.

The quiet time distribution of neuronal avalanches is characterized by a peculiar non-monotonic behavior, with power law decay followed by a local minimum and a more or less pronounced peak at a characteristic slow time scale (Lombardi et al., 2012). Numerical simulations suggest that such a distribution reflects the alternation between up and down-states in the network, which acts as a homeostatic mechanism controlling network excitability (Lombardi et al., 2012). In the current work, we analyze the functional behavior of the quiet time distribution in relation to the structure of avalanche sequences. In particular, we examine the relationship between quiet times and avalanche sizes by studying the distributions *P*(Δ*t*; *s*_{c}) of quiet times between consecutive avalanches of sizes larger than a given threshold *s _{c}* and investigate whether the non-monotonic quiet time distribution identified in cortex cultures exhibits the universal scaling features reported for waiting time distributions of earthquakes. We first compare quiet and waiting time statistics for neuronal avalanches. Then we show that, (1) the avalanche process in the up-state is solely controlled by the average occurrence rate and the corresponding quiet time distribution has a universal, i.e., sample independent, power law decay. By systematically removing smaller avalanches from the experimental time series, (2) we then unveil correlations between sizes and quiet times and highlight that avalanche occurrence exhibits some of the characteristic periodicity of θ (4–15 Hz), β (15–30 Hz), and γ (30–100 Hz) oscillations. Indeed, in place of the original power law, we observe several peaks at short time scales when considering only avalanches with size

*s*above a given threshold

*s*. Therefore, close in time smaller avalanches are crucial for the power law in the quiet time distribution of up-states to emerge. We observe that these avalanches tend to be related to short quiet times and fast β/γ oscillations, while larger avalanches are associated with slower θ and 1–2 Hz oscillations. In particular, we notice a sort of hierarchical structure in avalanche sequences: In the up-states, large avalanches occurring with θ frequency trigger cascades of smaller avalanches corresponding to faster oscillations. Finally we demonstrate (3) that the distributions

_{c}*P*(Δ

*t*;

*s*

_{c}) of quiet times between avalanches with size

*s*above a given threshold

*s*do not collapse if quiet times are rescaled by the average rate

_{c}*r*= 1/〈Δ

*t*〉. However, when the different temporal scales that govern up and down-states are taken into account, a proper collapse can be obtained. Specifically, the distributions

*P*(Δ

*t*;

*s*

_{c}) in the up-state and in the down-state show the same functional behavior if quite times are rescaled by the respective average avalanche rate.

## 2. Materials and Methods

### 2.1. Experimental Setup

Coronal slices from rat dorsolateral cortex (postnatal day 0–2; 350 μm thick) are attached to a poly-D-lysine coated 60-microelectrode array (MEA; Multichannelsystems, Germany) and grown at 35.5 C in normal atmosphere in standard culture medium without antibiotics for 4–6 weeks before recording. Avalanche activity is measured from cortex-striatum-substantia nigra triple cultures or single cortex cultures as reported previously (Beggs and Plenz, 2003). In short, spontaneous avalanche activity is recorded outside the incubator in standard artificial cerebrospinal fluid (ACSF; laminar flow of 1 ml/min) under stationary conditions for up to 10 h. The spontaneous local field potential (LFP) is sampled continuously at 1 kHz at each electrode and low-pass filtered at 50 Hz. Negative deflections in the LFP (nLFP) were detected by crossing a noise threshold of −3 *SD* followed by negative peak detection within 20 ms. nLFP times and nLFP amplitudes were extracted. Neuronal avalanches are defined as spatio-temporal clusters of nLFPs on the MEA (Beggs and Plenz, 2003). A neuronal avalanche consists of a consecutive series of time bins of width ϵ that contain at least one nLFP on any of the electrodes. Each avalanche is preceded and ended by at least one time bin with no activity. Without loss of generality, the present analysis is done with width ϵ individually estimated for each culture from the average inter nLFP interval on the array at which the power law in avalanche sizes s, *P*(*s*) ~ *s*^{−α}, yields α = 3/2. ϵ ranged between 3 and 6 ms for all cultures. Avalanche size is defined as the sum of absolute nLFP amplitudes (μV) on active electrodes or simply the number of active electrodes. Size distributions are obtained using logarithmic binning for sizes expressed in μV. A quiet time Δ*t* is defined as the time interval between the ending time of an avalanche *t ^{f}_{j}* and the starting time

*t*

^{i}

_{j + 1}of the following one, namely Δ

*t*=

_{j}*t*

^{i}

_{j + 1}−

*t*. A waiting time δ

^{f}_{j}*t*is defined as the time interval between the starting time of an avalanche

*t*and the starting time

^{i}_{j}*t*

^{i}_{j + 1}of the following one, namely δ

*t*=

_{j}*t*

^{i}_{j + 1}−

*t*. Quiet (waiting) time distributions are obtained using logarithmic binning for quiet (waiting) times expressed in ms.

^{i}_{j}#### 2.1.1. Up and down-state

The following procedure is used to discriminate between up and down-states. An up-state consists of a consecutive series of avalanches separated by a quiet time Δ*t* shorter than Δ*t**, where Δ*t** is defined as the local minimum between the initial power law regime and the local peak observed between 500 and 1000 ms. Conversely, every quiet time longer than Δ*t** belongs to a down-state and a consecutive series of avalanches separated by quiet times longer than Δ*t** is considered a down-state. The mean rate in the up-state is defined as *r _{up}* = 1/〈Δ

*t*〉

*, whereas the mean rate in the down-state is defined as*

_{up}*r*= 1/〈Δ

_{dw}*t*〉

*.*

_{dw}### 2.2. Numerical Model

#### 2.2.1. Network and dynamics

We consider *N* = 64000 neurons at random positions, characterized by their potential *v _{i}*. Neurons are connected by a scale-free network of synapses. More precisely to each neuron

*i*we assign an out-going connectivity degree,

*k*∈ [2, 100], according to the degree distribution

_{outi}*P*(

*k*) ∝

*k*

^{−2}of the functional network measured in Eguiluz et al. (2005). Choosing different network topologies, the model exhibits the same scaling behavior of avalanche size and duration distributions (de Arcangelis et al., 2006b; Pellegrini et al., 2007; de Arcangelis and Herrmann, 2012). The universality class of the neuronal avalanche process is the one of the mean field branching process (Zapperi et al., 1995; Lauritsen et al., 1996). To each synaptic connection we assign an initial random strength

*g*∈ [0.15, 0.3] and to each neuron an excitatory or inhibitory character. Outgoing synapses are excitatory if they belong to excitatory neurons, inhibitory otherwise. The network has a fraction

_{ij}*p*of inhibitory synapses, which is fixed. Each synapse is directed, meaning that it can be used by neuron

_{in}*i*to send a signal to neuron

*j*but not viceversa. As a consequence

*g*≠

_{ij}*g*and in general out-degree and in-degree of a neuron do not coincide. Therefore, once the network of output connections is established, we identify the resulting degree of in-connections,

_{ji}*k*, for each neuron

_{inj}*j*, namely we identify the number of synapses directed to each neuron

*j*. The number

*k*of in-going synapses can be considered as the dentritic tree of neuron

_{inj}*j*. We then assume that each neuron

*j*has a soma whose surface is proportional to

*k*.

_{inj}Whenever at time *t* the value of the potential in neuron *i* is above a certain threshold, *v _{i}* ≥

*v*, the neuron fires and its potential

_{max}*v*arrives at each of the

_{i}*k*connected neurons. In our simulations we use

_{outi}*v*= 6. However, as in every SOC-like model, this parameter is not relevant and results are independent of this particular choice.

_{max}For real neurons the production of neurotransmitters at the presynaptic terminals, and then the charge entering the postsynaptic neuron, is controlled by the frequency of action potentials, which depends on the integrated stimulation received by the neuron. Here the integrated stimulation is given by *v _{i}*, the membrane potential of the firing neuron. Therefore, we assume that the total charge

*q*that can enter into connected neurons is proportional to

_{i}*v*·

_{i}*k*. The change in the intracellular membrane potential of the postsynaptic neuron

_{outi}*j*is proportional to the relative synaptic strength

*g*/∑

_{ij}*,*

_{l}gilIn Equation 1 it is assumed that the received charge is distributed over the surface *k _{inj}* of the soma of the post-synaptic neuron. The plus or minus sign is for excitatory or inhibitory synapses respectively. After firing, the neuron is set in a refractory state lasting

*t*= 1 time step, during which it is unable to receive or transmit any charge, and its membrane potential is set to

_{ref}*v*= 0.

_{rest}#### 2.2.2. Avalanche activity

When a neuron fires, it may bring to threshold some of the connected neurons thus generating an avalanche, a cascade of activity which propagates through the network involving a variable number of neurons. During an avalanche there is no further external stimulation. As soon as no more neurons are able to fire, the avalanche ends and size is recorded as the number of firing neurons *s*, or, alternatively, as the sum *s*_{ΔV} of all positive potential variations (depolarizations) δ*v ^{+}_{i}* occurred in the network, namely

*s*

_{ΔV}= ∑

_{i}δ

*v*

^{+}

_{i}. By definition a single neuron firing does not constitute an avalanche. Avalanches are also characterized by their duration

*T*, which is defined as the number of iterations taken by the activity propagation. The numerical time step for each iteration corresponds to the real time between the triggering of an action potential in the presynaptic neuron and the change of the membrane potential in the postsynaptic neuron, therefore it is of the order of 4–6 ms. After an avalanche ends, an external stimulus triggers further activity in the system. Distributions of sizes and durations are shown in Supplementary Figure 1.

#### 2.2.3. Synaptic plasticity

We implement a Hebbian-like plasticity rule at the end of each avalanche. The strength *g _{ij}* of the used connections is increased proportionally to the membrane potential variation |δ

*v*| of the postsynaptic neuron

_{j}*j*induced by the presynaptic neuron

*i*during the avalanche,

whereas the strength of all inactive synapses is reduced by the average strength increase per bond

where *N _{B}* is the number of bonds. We set a minimum and a maximum value for the synaptic strength

*g*= 0.0001 and

_{ij}, g_{min}*g*= 1.0. Whenever

_{max}*g*<

_{ij}*g*, synapse

_{min}*g*is pruned. Since cortical plasticity such as long-term potentiation acts on time scales of seconds to minutes, which is much longer than the duration of avalanches, we apply the plasticity protocol for a certain number of stimulations and then study avalanche activity without further changing synaptic strengths. Specifically, since we don't want to alter the scale-free connectivity of the initial network, we apply plasticity rules until the first few synapses are pruned. After this plastic adaptation the

_{ij}*g*are distributed as shown in Supplementary Figure 2.

_{ij}#### 2.2.4. Up-down state dynamics

Alternation between the up and down-state was simulated on the basis of two concepts. First, the transition from one state to the other has a high degree of synchronization. Second, a down-state occurs when activity in the up-state reaches a level at which the up-state can't sustain itself anymore. Such a decrease in activity can result from either the exhaustion of available synaptic vesicles (Staley et al., 1998) or the increase of factors inhibiting their release (Thompson et al., 1992). For simplicity, we assume that the transition happens after a sufficiently large avalanche, which causes a lack of available neurotrasmitters and a sufficiently strong network inhibition.

Accordingly, at the end of each avalanche we measure its size in terms of the sum of depolarizations δ*v*^{+}_{i} of all neurons, *s*_{ΔV}. As soon as avalanche is larger than a threshold *s ^{min}*

_{ΔV},

*s*

_{ΔV}>

*s*

^{min}_{ΔV}, the system transitions into a down-state and neurons become hyperpolarized proportionally to their previous activity; namely, we reset

This rule models the local inhibition experienced by a neuron, due to spike adaptation (Sanchez-Vives et al., 2000), adenosine accumulation (Thompson et al., 1992), synaptic vesicles depletion (Staley et al., 1998) or blockade of receptor channels by the presence of external magnesium (Maeda et al., 1995). The down-state ends whenever a new avalanche occurs, namely the system transitions in an up-state. When in the up-state, all neurons firing in the previous avalanche of size *s*_{ΔV} are set to the depolarized value

This equation states that the neuron's intracellular membrane potential depends on the response of the whole network via *s*_{ΔV} and implements an homeostatic mechanism at the single neuron level: When avalanche sizes *s*_{ΔV} are close to the threshold *s ^{min}*

_{ΔV}, the ratio

*s*

_{ΔV}/

*s*

^{min}_{ΔV}is close to 1 and membrane potentials are reset closer to a zero resting value, thus avoiding an explosive growth of the following avalanche. Conversely, the network does sustain the depolarized state of the single neuron and the membrane potential stays closer to the firing threshold. We wish to stress that this mechanism is driven by the whole network activity, following the idea that the up-state in the cortex is a cooperative network state (Wilson, 2008). Furthermore, it is in agreement with measurements of the neuronal membrane potential, which remains significantly depolarized in the up-state (Wilson, 2008), and, at the same, keeps activity balanced. Through Equation (5), the threshold

*s*

^{min}_{ΔV}controls the level of excitability of the system.

At the network level, the high activity in up-states is sustained by a stimulation which has a random value in the interval *d _{u}* = [0,

*s*

^{min}_{ΔV}/

*s*

_{ΔV}): After an avalanche, at each time step we randomly choose a neuron and increase its membrane potential by

*rad*·

*d*, where

_{u}*rad*is a random number in the interval [0, 1). We notice that the amplitude of

*d*depends on past network activity through the size of the previous avalanche

_{u}*s*

_{ΔV}. As for Equation 5, the stimulation in the up-state is based on an homeostatic principle: The larger the previous avalanche the smaller

*d*and viceversa.

_{u}Conversely, during the down-state, the system experiences a general disfacilitation mimicked by weak random stimulation: At each time step we randomly choose a neuron and increase its membrane potential by a small constant quantity (30–40 smaller than *v _{max}*). This drive reproduces the effect of the small depolarizations due to miniature potentials (minis) from spontaneous synaptic release observed in the down-state (Timofeev et al., 2001). The drive slowly brings the system back in an up-state not correlated to past activity (Lombardi et al., 2012).

During the avalanche propagation the drive is stopped, as in usual SOC models. This procedure implements the separation of time scales between fast avalanche propagation and slow neuron stimulation.

Equations 4 and 5 each depend on a single parameter, *h* and *s ^{min}*

_{ΔV}, which introduce a memory effect at the level of single neuron activity and the entire system, respectively. In order to reproduce the experimentally observed behavior we only need to control the ratio

*R*=

*h/s*

^{min}_{ΔV}, as shown in Lombardi et al. (2012).

## 3. Results

### 3.1. Waiting Time and Quiet Time Distributions

The definition of quiet time and waiting time is sketched in Figure 1A and can be summarized in the following equality:

that is the *j*th waiting time is obtained summing up the *j*th quiet time and the duration *T _{j}* of the

*j*th avalanche. It follows that δ

*t*≃ Δ

*t*if the relation

*T*≪ Δ

*t*holds. In case of neuronal avalanches durations

*T*range from a few to few tens of milliseconds (Figure 1B) and are then comparable with the shortest Δ

*t*s (Lombardi et al., 2012). As a consequence, we expect quiet time and waiting time distribution to differ at short time scales. In Figure 1 we show the distribution

*P*(Δ

*t*) of quiet times between successive avalanches in six different cortex slice cultures (Lombardi et al., 2012) and compare them to the corresponding distributions

*P*(δ

*t*) of waiting times. The quiet time distribution has been extensively discussed in Lombardi et al. (2012), where it was called waiting time distribution. Here we briefly recall its main features, namely the power law behavior at short time scales, from few to 200–300 ms, and a local maximum situated at longer time scales, which leads to a peculiar non-monotonic behavior. The initial power law decay indicates that avalanches are temporally correlated if sufficiently close in time, which requires a sustained synaptic and firing activity in the network, namely an up-state. Conversely longer quiet times correspond to down-states and sparse synaptic activity in the network.

**Figure 1. Distributions of duration T, quiet times Δt and waiting times δt for six cortex slice cultures**.

**(A)**Definition of avalanche, quiet time and waiting time. nLFPs in the same time bin ϵ or consecutive bins define an avalanche. Avalanche duration

*T*is given by the number

*n*of consecutive non-empty bins times the bin amplitude ϵ, namely

*T*=

*n*· ϵ. A quiet time Δ

*t*is the time interval between the end of an avalanche and the start of the following one. A waiting time δ

*t*is the time interval between the start of an avalanche and the start of the following one. The following equality holds: δ

*t*= Δ

*t*+

*T*.

**(B)**Duration distributions. For better comparison duration

*T*is expressed in multiples of ∈. The initial power law regime extends for about one order of magnitude and is followed by an exponential cutoff.

**(C)**Distribution of quiet times: All curves show an initial power law regime with an exponent μ ranging between −2.0 and −2.5. For larger Δ

*t*, distributions are characterized by a local minimum followed by a more or less pronounced maximum at Δ

*t*≃ 1 − 2 s. Upper inset: Distributions of waiting times. Lower inset: illustrative comparison between quiet (cyan) and waiting (blue) time distribution for the blue curve in the main panel. The two distributions only differ at short time scales where durations are comparable to quiet times.

This non-monotonic behavior, with the same general features, can be still observed in the waiting time distributions (Figure 1C, upper inset). However, the power law exponent is generally slightly lower than the one measured for *P*(Δ*t*), as shown in the lower inset of Figure 1C. On the other hand, for time intervals larger than 200–300 ms, which are related to down-states, the two distributions basically coincide (Figure 1C, lower inset), meaning that, for this range of values, *T* ≪ Δ*t* and δ*t* ≃ Δ*t* is a good approximation. From 6 it follows that the waiting time distribution *P*(δ*t*) results from the combination of two quantities, quiet times and durations. While for long time scales *P*(δ*t*) is dominated by the former, at short time scales both of them contribute to its functional behavior. In this range of values both Δ*t* and *T* are power law distributed and add up to give again a power law: Short durations significantly couple with short quiet times and, due to lack of characteristic values, the net results is a power law with a larger slope. Is this power law carrying the same information as the statistics of time intervals without activity, i.e., quiet times? Evidently it does not, for the following reason: Durations, which are power law distributed, are not negligible and concluding that avalanches are temporally correlated from the power laws in waiting time distribution would be misleading (Sanchez et al., 2002). Indeed in Lombardi et al. (2012) only quiet time statistics has been considered. Nevertheless, some specific information can be extracted from waiting time distributions, as we will discuss in the following.

### 3.2. Temporal Features of Up and Down-State

In Lombardi et al. (2012) we have used numerical simulations to investigate the origin of the non-monotonic behavior in the quiet time distribution and concluded that it arises from the slow alternation of up and down-state. Accordingly, in this section we systematically isolate each contribution to the overall quiet time distributions (see Materials and Methods) and further investigate the temporal features of these two network states.

In Figure 2A we show the experimental distributions of quiet times between consecutive avalanches in the up-states (panel a) (see Materials and Methods): After rescaling Δ*t* by the mean rate *r _{up}* in the up-state, distributions collapse onto a unique power law with exponent μ ≃ −2.2. This implies that the avalanche process in the up-state is solely controlled by the average occurrence rate and the corresponding quiet time distribution has a universal, i.e., sample independent, power law decay (Figure 2A). On the other hand, down-states produce long quiet times mostly contributing to the tail of the overall

*P*(Δ

*t*), exhibiting a distribution with a characteristic value τ

_{d}, as found numerically (Lombardi et al., 2012). This behavior has a simple interpretation: The recurrence of up-states has a more or less pronounced characteristic time. If the distribution of quiet times in the down-state is peaked around a particular value τ

_{d}and is sufficiently narrow, then a non-monotonic behavior can be observed in the quiet time distribution of the entire avalanche activity. Although distributions of quiet times in the down-states do exhibit common features across samples, they do not generally collapse onto a unique function after rescaling δ

*t*by

*r*, the mean rate in the down-state (Figure 2B).

_{dw}**Figure 2. Experimental distribution of quiet times in the up-state and in the down-state**. **(A)** Distribution of quiet times between successive avalanches occurring in the up-state. The curves, rescaled by the mean rate *r _{up}*, show a universal power law scaling. The dashed line represents a power law with exponent −2.2.

**(B)**Distribution of quiet times between successive avalanches occurring in the down-state. In this case, rescaling by the mean rate

*r*does not lead to a universal behavior.

_{dw}To complete the investigation of up and down-state temporal features, we consider the distributions *P*(*T _{up}*) and

*P*(

*T*) of up and down-state durations, respectively (Figure 3). Numerical curves are over plotted with experimental results. We notice here that, both numerically and experimentally, the two states are characterized by time scales that differ by about one order of magnitude. Moreover, their respective duration distributions exhibit a distinct functional behavior. On average, the durations of down-states are distributed around

_{dw}*T*≃ 2000 ms and the tail of the distribution is well fitted by an exponential (Millman et al., 2010). This property characterizes most of the analyzed samples (Supplementary Figure 3). Conversely, the distribution

_{dw}*P*(

*T*) exhibits a tail compatible with a power law. However, in this case, the power law behavior arises by averaging over many cultures and does not necessarily characterize the up-state duration in each culture (Supplementary Figure 4).

_{up}**Figure 3. Distribution of durations of up-states (A) and down-states (B) averaged over 100 configurations of a network of N = 64000 neurons with p_{in} = 0.1 (black symbols)**. Experimental data are averaged over all samples (green curves).

### 3.3. Temporal Structure of Avalanche Process

We have shown that different quiet time distributions of distinct experimental samples show a qualitatively similar behavior. In particular, at short time scales, the distributions of quiet times are all characterized by the same power law (Figure 2), a general and robust feature of up-states. Here we go further in the characterization of the avalanche process and question how the functional behavior of the distribution *P*(Δ*t*) depends on the fine structure of avalanche sequences. In order to do that, we study the distributions *P*(Δ*t*; *s*_{c}) of quiet times between consecutive avalanches of size larger than a given threshold *s _{c}*. In this way we remove smaller avalanches from the time series and analyze how the distribution changes as a function of

*s*. If the different distributions

_{c}*P*(Δ

*t*;

*s*

_{c}) collapse onto a unique function, then the temporal properties of the avalanche process are invariant under the aforesaid removal procedure. This specific point will be addressed in the next section.

In Figure 4 we show the distribution *P*(Δ*t*; *s*_{c}) for different values of *s _{c}*. By removing avalanches we are making the time series sparser and, as a result, we would expect the distributions

*P*(Δ

*t*;

*s*

_{c}) to become broader and broader as we increase

*s*. Indeed this effect is observed but it is minor for a wide range of

_{c}*s*values, which suggests that large quiet times tend to separate large avalanches. On the contrary, as a main effect, we observe that the distributions

_{c}*P*(Δ

*t*;

*s*

_{c}) show peaks that were not present in the original

*P*(Δ

*t*). These peaks become pronounced for values of

*s*larger than 40 μV and are either located at the time scales within the power law regime or at very long quiet times. The first peak appears at Δ

_{c}*t*≃ 40−60 ms (Δ

*t*

_{β}in the following) and can be related to the period of β oscillations (Figure 4C). The second one arises at Δ

*t*∈ [80, 250] ms (Δ

*t*

_{θ}in the following) and corresponds to the period of θ oscillations. This peak is visible in all samples. In particular it is very pronounced in Figures 4B–D,F. Quiet times around Δ

*t*

_{θ}seem to play a special role with respect to our removal process: While the probability increases with

*s*for Δ

_{c}*t*longer than Δ

*t*

_{θ}and decreases for the shorter ones, it stays nearly constant in a neighborhood of Δ

*t*

_{θ}(Figures 4B–D,F). This means that the ratio

*N*(Δ

*t*

_{θ};

*s*)/

_{c}*N*(

*s*) ≃

_{c}*const*, namely the number of quiet times corresponding to θ period scales with the total number

*N*(Δ

*t*;

*s*). Since the number of avalanches larger than

_{c}*s*is simply given by the number of quiet times plus one, then the number

_{c}*N*

_{θ}(

*s*) of avalanches related to θ oscillations scales with the total number

_{c}*N*(

*s*) of avalanches, namely it decreases proportionally to

_{c}*N*(

*s*) for increasing values of

_{c}*s*. On the other hand, the number of avalanches separated by longer and shorter quiet times decreases slower and faster than

_{c}*N*(

*s*), respectively. This point can be understood as follows. If, for a given Δ

_{c}*t*, the probability

*P*(Δ

*t*;

*s*) increases (decreases) with

_{c}*s*(Figures 4B–D,F), then the numerator of the ratio

_{c}*N*(Δ

*t*;

*s*)/

_{c}*N*(

*s*) decreases slower (faster) than the denominator and so does the corresponding number of avalanches. Alternatively, one can look at the quantity

_{c}*N*(Δ

*t*;

*s*

_{c}), which we show in the Supplementary Figure 5, and notice that it decreases faster for small than for large Δ

*t*s. Therefore, long quiet times tend to occur between large avalanches whereas shorter quiet times tend to separate the smaller ones. From Figure 4 we notice that, whenever the peak around Δ

*t*

_{θ}is not pronounced (Figures 4A,E), the Δ

*t*characteristic of slow 1 Hz oscillations between up and down-states plays the role of fixed point. Finally a further peak appears at Δ

*t*≃400-500 ms, which corresponds to a ≈ 2 Hz oscillation (Figure 4B). This peak behaves as the one at Δ

*t*

_{θ}, namely it behaves as a fixed point for our removal procedure.

**Figure 4. Experimental quiet time distributions for different values of the threshold s_{c} on avalanche size**. Already for

*s*= 80 μV, distributions clearly exhibit one or more additional peaks. Beside the one at large time scales, Δ

_{c}*t*≃ 1000–2000 ms, which is related with the characteristic time of up-state recurrence, at least one peak between 60 and 250 ms is always visible on time scales originally characterized by the power law decay and corresponds to the period of β or θ oscillations. This is particularly pronounced in

**(B,C,D,F)**, less in

**(A,E)**. The distributions in

**(B,C)**exhibit one more peak around 500 ms, related to 2 Hz oscillations. It is worth to notice that the probability

*P*(Δ

*t*) for Δ

*t*corresponding to the θ

**(B,C,D,F)**and 1–2 Hz oscillations

**(A,B)**is nearly a fixed point for this transformation. Insets: Experimental distributions of waiting times δ

*t*, for different values of

*s*. In this case one more peak appears around 20–30 ms, which corresponds to γ oscillations.

_{c}Since avalanche durations and periods of fast oscillations are of the same order of magnitude, in order to capture their relation with avalanche sizes we have considered the distributions *P*(δ*t*; *s*_{c}) of waiting times between consecutive avalanches of size larger *s _{c}*, which are shown in the insets of Figure 4. The picture emerging from the analysis of the quantity δ

*t*is basically the same we have drawn looking at the quiet times, except for a peak corresponding to the faster γ oscillations, which can be now clearly observed in the insets of Figures 4A,B,D–F of Figure 4. The probability associated with this peak, which is situated at very short δ

*t*, decreases with

*s*whenever it coexists with very pronounced θ peaks (Figures 4B–D,F), indicating that, at least in this particular case, faster oscillations tend to be associated with smaller avalanches.

_{c}To summarize, our removal procedure uncovers a rich temporal structure hidden behind the scale free behavior in the quiet time distribution: Beside the characteristic time associated with down-state duration, avalanche occurrence keeps the temporal features of θ and β/γ oscillations. They jointly emerge in most of the analyzed experimental samples (Figures 4B–D,F). While short quiet times and fast β/γ oscillations tend to be associated with smaller avalanches, slower oscillations are in general related to larger avalanches, but without any characteristic size. Indeed, varying the threshold *s _{c}* in a range of values within the power law regime of the size distribution

*P*(

*s*), typically between 30 μV and 400 μV (Figure 5B), the probability

*P*(Δ

*t*;

*s*

_{c}) of Δ

*t*associated with θ (Figures 4B–D,F) or slower oscillations (Figures 4A,E) remains nearly unchanged. In particular, the θ peak coexists with a faster decrease of the probability of γ period, thus suggesting that a sort of hierarchical structure for avalanche sequences, which follows closely the temporal organization of nested θ − β/γ oscillations: Within up-states, large avalanches occur with θ frequency and trigger smaller ones in a faster γ cycle (Figure 5A). Remarkably, within γ cycles the quiet times have no characteristic value. Indeed the quiet time distributions do not show peaks at very short time scales. Then, quiet times and durations, which are both power law distributed, show a peculiar coupling in the γ cycles. δ

*t*s corresponding to these oscillations are short, which implies that both

*T*and Δ

*t*are short. Considering the scaling relation between duration

*T*and

*s*(Friedman et al., 2012), this is the same that saying small avalanches are associated with short quiet times.

**Figure 5. Neuronal avalanches organize into a hierarchical structure corresponding to temporal organization of nested θ − β/γ oscillations**. **(A)** During up-states large avalanches (blue bars) occur with θ frequency and trigger smaller avalanches related to faster γ oscillations (green bars). Here bar widths indicate durations: Avalanche start is at the right side of the bar. Bar heights indicate sizes. Spacing between blue bars corresponds to a θ period. Spacing between the starting points of green bars corresponds to γ period. γ cycles do not show characteristic quiet times. Sizes *s* of avalanches related to θ cycles tend to fall within the blue region of the size distribution *P*(*s*) plotted in **(B)**, whereas the ones corresponding to nested γ oscillations fall within the green region. Therefore, the relationship between avalanches and oscillations does not imply characteristic sizes. In particular, for *s _{c}* ≥ 80 μV, the number of avalanches

*N*

_{θ}related to θ cycles scales with

*s*as the total number

_{c}*N*of avalanches, namely

*N*

_{θ}/

*N*≃

*const*.

**(B)**Distributions of avalanche sizes for the experimental samples in Figure 1. Same color is used here for each sample.

Figure 4 indicates that quiet times and avalanche sizes are correlated. The analysis of the scatter plots between Δ*t* and the relative previous and following avalanche also provides some evidence that correlations exist (Supplementary Figures 6, 7). In order to further validate this result, we reshuffle avalanche sizes while keeping the sequence of starting and ending times fixed. More precisely, we reassign to each avalanche a size taken at random from the measured size distribution. Then, we apply the same procedure described above. If no correlations existed between sizes and waiting times, then we should still observe the same peaks in the distributions *P*(Δ*t*; *s*_{c}). As shown in Figure 6, in this case no peaks emerge in the power law regime, which implies that, in the up-state, waiting times are strongly correlated with sizes. In particular, periods of θ, β, and γ oscillations are correlated with sizes of corresponding avalanches. On the other hand, for longer waiting times we observe the same qualitative behavior discussed for the original time series. Therefore, we can state that correlations with avalanche sizes are weak, but a more quantitative analysis is needed to exclude that they are significant.

**Figure 6. Quiet time distributions evaluated for the reshuffled avalanche time series and for different values of the threshold s_{c} on avalanche size**. In this case no additional peaks arise at short time scales. Distributions still exhibit peaks at longer time scales, as the ones shown in Figure 4.

#### 3.3.1. Up and down-state

From the analysis performed above it is evident that the functional behavior of the quiet time distribution arises from the superposition of many dynamic mechanisms. In Section 3.2 we have argued that non-monotonicity results from the alternation between up and down-state, which implies already two different mechanisms governing avalanche activity at the short and large time scales. Then we have shown that also the characteristic times of θ, β, and γ oscillations enter in the process. As a consequence, we do not expect the distribution *P*(Δ*t*; *s*_{c}) being controlled by a single parameter, as observed in other stochastic time series (Corral, 2004). Indeed, rescaling the quiet times by the mean avalanche rate *r* = 1/〈Δ*t*〉, does not lead to a collapse of the curves onto a single one (not shown).

However, one can apply the same removal procedure separately to up and down-states and then rescaling quiet times by the respective average occurrence rate, *r _{up}* and

*r*, in order to find universal features for each of the two network states. We start considering the distributions

_{dw}*P*(Δ

*t*;

*s*

_{c}) in the down-state and we rescale them by

*r*= 〈Δ

_{dw}*t*〉

*. As shown in Figure 7, distributions collapse onto a unique function, which shows a characteristic value and an exponential tail. This functional behavior is common to all samples except the one in Figure 7E, whose departure from an exponential could be interpreted as an effect of the very sharp peak at Δ*

_{dw}*t*≃ 1 s and not as a result of a different dynamics in the down-state. The existence of a universal function implies that the quiet time distribution in the down-states is uniquely controlled by

*r*. On the other hand, following the same procedure for the up-state does not provide a good data collapse (not shown). Peaks that emerge at short Δ

_{dw}*t*after the removal of smaller avalanches, tell us that avalanche occurrence in the up-state is not solely controlled by one time constant, that is 1/

*r*. Nevertheless, here we show that the distributions of quiet times shorter than Δ

_{up}*t*

_{θ}are solely controlled by

*r*

^{θ}

*= 1/〈Δ*

_{up}*t*〉

_{Δt < Δtθ}, where 〈·〉

_{Δt < Δtθ}indicates the average over Δ

*t*< Δ

*t*

_{θ}. Indeed, rescaling them by

*r*

^{θ}

*leads to a data collapse onto a unique function which follows a power law with an exponent μ ≃ −2 (Figure 8). This collapse is particularly good in samples that show a clear power law behavior for quiet times shorter than Δ*

_{up}*t*corresponding to θ and 1 Hz oscillation period (Figures 8B–E,F). Conversely, curves do not collapse whenever a further, shorter characteristic time is present (Figure 8C).

**Figure 7. Distributions of quiet times P(Δt; s_{c}) in the down-state for the experimental data samples of Figure 1C and the numerical samples reproducing blue squares and red diamonds curve of Figure 1C**. Distributions are rescaled by the mean rate

*r*in the down-state. In five of the analyzed samples the tail of the distribution is well fitted by an exponential (black dashed line in

_{dw}**A–D,F**). Numerical data are shown in

**(E,F)**together with the corresponding experimental curves of Figure 1C and shifted by 1 order of magnitude to the left, for clarity. Numerical distributions are averaged over 100 configurations of a network of

*N*= 64000 neurons with

*p*= 0.1.

_{in}**Figure 8. Distributions P(Δt < Δt_{θ}; s_{c}) of quiet times shorter than Δt_{θ} in the up-state for the experimental data samples of Figure 1C and the numerical samples reproducing blue squares and red diamonds curve of Figure 1C**. Δ

*t*

_{θ}is sample dependent and its value varies in the interval [80, 250] ms. Distributions are rescaled by the mean rate

*r*

^{θ}

*. Numerical data are shown in*

_{up}**(E,F)**together with the corresponding experimental curves of Figure 1C and shifted by 1 order of magnitude to the left, for clarity. Numerical distributions are averaged over 100 configurations of a network of

*N*= 64000 neurons with

*p*= 0.1 and are rescaled by

_{in}*r*. The dashed line is a power law with exponent −2.2.

_{up}We obtain a similar result for numerical distributions. However, in this case removing avalanches according to their size does not lead to many peaks at short quiet times, which implies that there are only two characteristic time scales for numerical avalanches. In this case, we just need to consider separately up and down-state and rescale the quiet times by respective average occurrence rate, *r _{up}* and

*r*. As shown in Figures 7E,F, 8E,F we obtain a good data collapse in both cases.

_{dw}## 4. Discussion

The distribution of quiet times between consecutive avalanches in cortex slice cultures displays a power law decay at short time scales, namely from few to 200–300 ms, and is generally characterized by a local maximum at longer quiet times, which leads to a non-monotonic behavior. Numerical simulations show that this non-monotonic distribution results from the slow alternation between up and down-states (Lombardi et al., 2012). The model suggests that in the up-state, where neurons mutually sustain their spiking activity, network mechanisms act as a form of short-term memory, which produces clusters of correlated avalanches and thus gives rise to the initial power law regime in the quiet time distribution. On the other hand, the synaptic activity during down-state can be modeled as a random process that slowly brings the system back into the up-state, with no memory of past activity. Indeed numerical distributions exhibit an exponential tail similar to the ones observed experimentally (Lombardi et al., 2012).

Accordingly, here we have defined as up-state (down-state) a consecutive series of avalanches separated by Δ*t* shorter (longer) than the longest Δ*t* falling within the power law regime of *P*(Δ*t*) and systematically evaluated the quiet time distribution for up and down-state. We have shown that, while a power law with exponent μ ≃ −2 is a property of up-states in all analyzed samples, the recurrence of up-states has a characteristic time τ_{d} which is sample dependent (≃ 1 s on average). Indeed, the lasting times of down-states, which are simply quiet times between successive up-states, are distributed around a certain value 1 s < *T* < 2 s, the tail of the distribution being well fitted by an exponential. Since the exponential behavior is characteristic of Poisson processes, we conclude that consecutive up-states are basically not correlated. Moreover, from the properties of Poisson processes it follows that, given the sequence of quiet times Δ*t* between successive up-states, the jumps, i.e., the differences between two consecutive Δ*t*s, are also exponentially distributed. The distribution of jumps is commonly used to characterize stochastic processes. It has been analyzed for burst sequences in spontaneous activity of dissociated cultures of cortical neurons (Segev et al., 2002) and has been approximated with a symmetric Levy distribution. While Levy is indicative of self similarity in the process, spectral analysis was consistent with long range temporal correlation. Beside differences with cultures considered here, discrepancies can be also due to the definition of burst adopted in Segev et al. (2002), which substantially differs from our definition of up-states.

We have shown that beside the characteristic recurrence time τ* _{d}* between consecutive up-states, the analysis of quiet time distributions is able to capture the presence of θ, β, and γ oscillations in avalanche occurrence. The connection between nested oscillations and neuronal avalanches has been pointed out in Gireesh and Plenz (2008). Investigation of spontaneous neuronal activity in the rat cortex layer 2/3 has revealed that, during the second week postnatal, bursts develop a temporal organization of higher frequency oscillations, β and γ, nested into lower frequencies θ oscillations, while the spatio-temporal organization of LFPs is characterized by the scaling behavior of neuronal avalanches. Here we have further enlightened the relation between avalanche sizes and the temporal structure of the avalanche process. When avalanches of all sizes are considered, the distribution of quiet times in the up-state is scale free. On the contrary, disregarding avalanches smaller than ≃ 80μV, peaks corresponding to oscillations in θ, β, and γ frequency bands are clearly visible. Smaller avalanches (60–160 μV) tend to be associated with shorter quiet times and faster β/γ oscillations, larger ones to longer quiet times and slower θ or 1–2 Hz oscillations. Of considerable interest is the behavior of the θ and 1 Hz peaks under the removal procedure, which are nearly independent of the threshold

*s*on avalanche sizes: It doesn't matter how many avalanches are removed, the probability for quiet times around the period of θ or 1 Hz oscillation does not change for a large range of

_{c}*s*values. Equivalently, avalanches corresponding to these frequency bands are a constant fraction of the total number, which implies that they have no characteristic size. This suggests a special role in the temporal organization of spontaneous activity. In particular, we have noticed that large avalanches occurring with θ frequency trigger cascades of smaller avalanches corresponding to the higher frequency oscillations, in a sort of hierarchy which is reminiscent of the temporal organization of nested θ − β/γ oscillations (Gireesh and Plenz, 2008; He et al., 2010).

_{c}These results indicate that correlations between quiet times and avalanche sizes could be relevant and deserve further investigation. This point is intimately related to the existence of a universal scaling function for the distributions *P*(Δ*t*; *s*_{c}). A stochastic process for which such a universal function exists is a fixed point of the transformation which has been illustrated and performed in Section 3 (Corral, 2007). It can be shown that the only process without correlations which is invariant under this transformation is the Poisson process (Daley and Vere-Jones, 1988). More precisely, if sizes are independent of any other variable, the removal of events is equivalent to a so called random thinning and, under certain conditions, the resulting process converges to a Poisson process. Here we have demonstrated that the distributions *P*(Δ*t*; *s*_{c}) do not collapse onto a unique function when Δ*t* is rescaled by the average occurrence rate *r*. This is because of the multiple time scales in avalanche dynamics, which result from different mechanisms governing avalanche triggering during up and down-states. Indeed distributions *P*(Δ*t*; *s*_{c}) for the down-state are simply controlled by the respective average rate: When Δ*t* is rescaled by *r _{dw}*, the distributions

*P*(Δ

*t*;

*s*

_{c}) for the down-state collapse onto the same curve with an approximately exponential tail, which therefore implies that sizes of avalanches separated by large quiet times are either independent or weakly correlated, as well as sizes and quiet times. On the other hand, in the up-state we observe that the peak associated with period of θ oscillations and those corresponding to the β/γ scale differently with

*s*and therefore cannot be controlled by the same time scale,

_{c}*r*. In other words oscillations introduce additional characteristic times in the up-state. However, we have shown that the power law for short quiet times is universal and controlled by 〈Δ

_{up}*t*〉

_{Δt<Δtθ}. A similar analysis has been recently performed for spike avalanches in freely behaving (FB) and anesthetized rats (AR) (Ribeiro et al., 2010), where the quiet time distributions show consistently a monotonically decreasing behavior. Universal scaling features are observed for FB rats when quiet times are rescaled by the average occurrence rate, whereas curves for AR do not collapse onto a unique function. Our analysis suggests that the different behavior between anesthetized and freely behaving rats could be due to different dynamic mechanisms characterizing spontaneous activity in AR.

Waiting time distribution and its universal features have been widely investigated for earthquakes (Corral, 2004; de Arcangelis et al., 2006a). In this case the distribution is not exponential, but monotonic and solely controlled by *r*, except for corrections at short waiting times (Bottiglieri et al., 2010). On the other hand, many similarities between neuronal avalanches and earthquakes can be recognized, which have suggested a common interpretation in term of self organized criticality (SOC). SOC was originally proposed as an explanation for long range correlations emerging in processes far from equilibrium (Tang et al., 1988) and has rapidly become a useful interpretative scheme for many stochastic natural phenomena that exhibit scale free statistics. As for neuronal avalanches and earthquakes, in many cases, e.g., solar flares (Boffetta et al., 1999), waiting time distributions are not exponential. Conversely, in the original sand pile model introduced by Bak, Tang and Wiesenfeld (BTW) to exemplify SOC idea, waiting times are exponentially distributed (Boffetta et al., 1999) and this fact was used to question SOC as an interpretation for solar flares (Boffetta et al., 1999) and earthquakes (Yang et al., 2004). However, Paczuski et al. (2005) have argued that an experimental sequence of bursts can arise from a single avalanche observed at a finite detection threshold, which would give rise to a power law in the waiting time distribution of the BTW model. In addition, several different models have been proposed in order to show that SOC-like dynamics can provide temporal correlations among avalanches (Rios and Zhang, 1999; Baiesi and Maes, 2006) and a non-exponential distribution of waiting times (Sanchez et al., 2002; Lippiello et al., 2005; Baiesi and Maes, 2006). In particular, it has been shown that in the so called running sand pile (Sanchez et al., 2002), waiting times between avalanches with size above a large enough threshold are power law rather than exponentially distributed. Non-exponential waiting time distributions also arise if avalanches are triggered on the basis of the entire history of local stimulations (Lippiello et al., 2005). Here we have shown that our model, inspired in SOC, is able to capture the peculiar, non-exponential and non-monotonic behavior of the waiting time distribution for neuronal avalanches recorded in cortex slice cultures (Lombardi et al., 2012). Moreover, numerically generated up and down-states, exhibit the same universal features found experimentally. This point is particularly important because it indicates that the lack of universality in the waiting time distribution for spike avalanches in anesthetized rats (Ribeiro et al., 2010) could be due to the coexistence of different dynamic mechanisms, each one controlling ongoing activity at different temporal scales. Indeed, in freely behaving rats, where no down-states are observed, the waiting time distribution is controlled by the average occurrence rate (Ribeiro et al., 2010), which, for our model, is equivalent to *r _{up}*. From our simulations it emerges that the crucial features of this temporal evolution are (1) the different single neuron behavior in the two phases, namely the ability to oscillate between a very depolarized and hyperpolarized state, (2) the homeostatic mechanism driving activity in the up-state and (3) the network disfacilitation following up-states. The good agreement with experimental data indicates that the transition from an up-state to a down-state has a high degree of synchronization, whereas the onset of up-states is usually more gradual. According to our numerical results, the alternation between up and down-states is the expression of an homeostatic regulation which, during a burst, is activated to control the excitability of the system and avoid pathological behavior.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We thank the SNF for funding within project 205321-138074.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fnsys.2014.00204/abstract

## References

Baiesi, M., and Maes, C. (2006). Realistic time correlations in sandpiles. *Europhys. Lett*. 75, 413–419. doi: 10.1209/epl/i2006-10137-2

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. *J. Neurosci*. 23, 11167–11177.

Boffetta, G., Carbone, V., Giuliani, P., Veltri, P., and Vulpiani, A. (1999). Power laws in solar flares: self-organized criticality or turbulence? *Phys. Rev. Lett*. 83, 4662. doi: 10.1103/PhysRevLett.83.4662

Bottiglieri, M., de Arcangelis, L., Godano, C., and Lippiello, E. (2010). Multiple-time scaling and universal behavior of the earthquake interevent time distribution. *Phys. Rev. Lett*. 104, 158501. doi: 10.1103/PhysRevLett.104.158501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Corral, A. (2004). Long-term clustering, scaling and universality in the temporal occurrence of earthquakes. *Phys. Rev. Lett*. 92, 108501. doi: 10.1103/PhysRevLett.92.108501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Corral, A. (2007). Statistical features of earthquake temporal occurrence. *Lect. Notes Phys*. 705, 191–221. doi: 10.1007/3-540-35375-5-8

Cossart, R., Aronov, D., and Yuste, R. (2003). Attractor dynamics of network up states in the neocortex. *Nature* 423, 283–288. doi: 10.1038/nature01614

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cowan, R. L., and Wilson, C. J. (1994). Spontaneous firing patterns and axonal projections of single corticostriatal neurons in the rat medial agranular cortex. *J. Neurophysiol*. 71, 17–32.

Cunningham, M. O., Pervouchine, D. D., Racca, C., Kopell, N. J., Davies, C. H., Jones, R. S. G., et al. (2006). Neuronal metabolism governs cortical network response state. *Proc. Natl. Acad. Sci. U.S.A*. 103, 5597–5601. doi: 10.1073/pnas.0600604103

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Daley, D. J., and Vere-Jones, D. (1988.) *An Introduction to the Theory of Point Processes*. New York, NY: Springer.

de Arcangelis, L., Godano, C., Lippiello, E., and Nicodemi, M. (2006a). Universality in solar flare and earthquake occurrence. *Phys. Rev. Lett*. 96:051102. doi: 10.1103/PhysRevLett.96.051102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

de Arcangelis, L., and Herrmann, H. J. (2012). Activity dependent neuronal model on complex networks. *Front. Physiol*. 3:62. doi: 10.3389/fphys.2012.00062

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

de Arcangelis, L., Perrone-Capano, C., and Herrmann, H. J. (2006b). Self-organized criticality model for brain plasticity. *Phys. Rev. Lett*. 96:028107. doi: 10.1103/PhysRevLett.96.028107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Eguiluz, V. M., Chialvo, D., Cecchi, G. A., Baliki, M., and Apkarian, A. V. (2005). Scale-free brain functional networks. *Phys. Rev. Lett*. 94:018102. doi: 10.1103/PhysRevLett.94.018102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Friedman, N., Ito, S., Brinkman, B. A. W., Shimono, M., Lee DeVille, R. E., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data. *Phys. Rev. Lett*. 108:208102. doi: 10.1103/PhysRevLett.108.208102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gireesh, D. E., and Plenz, D. (2008). Neuronal avalanches organized as nested theta-and beta/gamma-oscillations during development of cortical layer 2/3. *Proc. Natl. Acad. Sci. U.S.A*. 105, 7576–7581. doi: 10.1073/pnas.0800537105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hahn, T. T. G., Sakmann, B., and Mehta, M. R. (2006). Phase-locking of hippocampal interneurons membrane potential to neocortical updown states. *Nat. Neurosci*. 9, 1359–1361. doi: 10.1038/nn1788

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

He, B. J., Zempel, J. M., Snyder, A. Z., and Raichle, M. E. (2010). The temporal structures and functional significance of scale free brain activity. *Neuron* 66, 353–369. doi: 10.1016/j.neuron.2010.04.020

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lauritsen, K. B., Zapperi, S., and Stanley, E. (1996). Self-organized branching processes: avalanche models with dissipation. *Phys. Rev. E* 54, 2483–2488. doi: 10.1103/PhysRevE.54.2483

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lippiello, E., de Arcangelis, L., and Godano, C. (2005). Memory in self-organized criticality. *Europhys. Lett*. 72:678. doi: 10.1209/epl/i2005-10292-x

Lombardi, F., Herrmann, H. J., Perrone-Capano, C., Plenz, D., and de Arcangelis, L. (2012). Balance between excitation and inhibition controls the temporal organization of neuronal avalanches. *Phys. Rev. Lett*. 108, 228703. doi: 10.1103/PhysRevLett.108.228703

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Maeda, E., Robinson, H. P., and Kawana, A. (1995). The mechanisms of generation and propagation of synchronized bursting in developing networks of cortical neurons. *J. Neurosci*. 15, 6834–6845.

Mazzoni, A., Broccard, F. D., Garcia-Perez, E., Bonifazi, P., Ruaro, M. E., and Torre, V. (2007). On the dynamics of the spontaneous activity in neuronal networks. *PLoS ONE* 2:e439. doi: 10.1371/journal.pone.0000439

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

McCormick, D. A., Shu, Y., Hasenstaub, A., Sanchez-Vives, M., Badoual, M., and Bal, T. (2003). Persistent cortical activity: mechanism of generation and effects on neuronal excitability. *Cereb. Cortex* 13, 1219–1231. doi: 10.1093/cercor/bhg104

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Millman, D., Mihalas, S., Kirkwood, A., and Niebur, E. (2010). Self-organized criticality occurs in non-conservative neuronal networks during ’up’ states. *Nat. Phys*. 6, 801–805. doi: 10.1038/nphys1757

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Paczuski, M., Boettcher, S., and Baiesi, M. (2005). Interoccurence times in the B.T.W. sandpile model: a comparison with the observed ststistics of solar flares. *Phys. Rev. Lett*. 95:181102. doi: 10.1103/PhysRevLett.95.181102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pasquale, V., Massobrio, P., Bologna, L. L., Chiappalone, M., and Martinoia, S. (2008). Self-organization and neuronal avalanches in networks of dissociated cortical neurons. *J. Neurosci*. 153, 1354–1369. doi: 10.1016/j.neuroscience.2008.03.050

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pellegrini, G. L., de Arcangelis, L., Herrmann, H. J., and Perrone-Capano, C. (2007). Activity-dependent neural network model on scale-free networks. *Phys. Rev. E* 76:028107. doi: 10.1103/PhysRevE.76.016107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Petermann, T., Thiagarajan, T. C., Lebedev, M., Nicolelis, M., Chialvo, D. R., and Plenz, D. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. *Proc. Natl. Acad. Sci. U.S.A*. 106, 15921–15926. doi: 10.1073/pnas.0904089106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Plenz, D. (2012). Neuronal avalanches and coherence potentials. *Eur. Phys. J. Special Top*. 205, 259–301. doi: 10.1140/epjst/e2012-01575-5

Plenz, D., and Aertsen, A. (1996). Neural dynamics in cortex-striatum co-cultures-ii. spatiotemporal characteristics of neuronal activity. *Neuroscience* 70, 893–924. doi: 10.1016/0306-4522(95)00405-X

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Plenz, D., and Kitai, S. T. (1996). Generation of high-frequency oscillations in local circuits of rat somatosensory cortex cultures. *J. Neurophysiol*. 76, 4180–4184.

Ribeiro, T. L., Copelli, M., Caixeta, F., Belchior, H., and Chialvo, D. R. (2010). Spike avalanches exhibit universal dynamics across the sleep-wake cycle. *PLoS ONE* 5:e14129. doi: 10.1371/journal.pone.0014129

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Rios, P. D. L., and Zhang, Y.-C. (1999). Universal 1/f noise from dissipative self-organized criticality model. *Phys. Rev. Lett*. 82, 472. doi: 10.1103/PhysRevLett.82.472

Sanchez, R., Newman, D. E., and Carreras, D. A. (2002). Waiting-time statistics of Self-Organized-Criticality Systems. *Phys. Rev. Lett*. 88:068302. doi: 10.1103/PhysRevLett.88.068302

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., Novak, L. G., and McCormick, D. A. (2000). Cellular mechanisms of long-lasting adaptation in visual cortical neurons. *J. Neurosci*. 20:4286.

Segev, R., Morris, B., Hulata, E., Cohen, N., Palevski, A., Kapon, E., et al. (2002). Long term behavior of lithographically prepared *in vitro* neuronal networks. *Phys. Rev. Lett*. 88:118102. doi: 10.1103/PhysRevLett.88.118102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Shriki, O., Alstott, J., Carver, F., Holroyd, T., Hanson, R. N., Smith, M. L., et al. (2013). Neuronal avalanches in the resting meg of the human brain. *J. Neurosci*. 33, 7079–7090. doi: 10.1523/JNEUROSCI.4286-12.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Stanley, H. E. (1971.) *Introduction to Phase Transition and Critical Phenomena*. New York, NY: Oxford University Press.

Staley, K. J., Longacher, M., Bains, J. S., and Yee, A. (1998). Presynaptic modulation of CA3 network activity. *Nat. Neurosci*. 1, 201–209. doi: 10.1038/651

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Steriade, M., Nunez, A., and Amzica, F. (1993). A novel slow (< 1Hz) oscillation of neocortical neurons *in vivo*: depolarizing and hyperpolarizing components. *J. Neurosci*. 13, 3252–3265.

Stevens, C., and Tsujimoto, T. (1995). Estimates for the pool size of releasable quanta at a single central synapse and for the time required to refill the pool. *Proc. Natl. Acad. Sci. U.S.A*. 92, 846–849. doi: 10.1073/pnas.92.3.846

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Tang, C., Bak, P., and Wiesenfeld, K. (1988). Self-organized criticality. *Phys. Rev. A* 38:364. doi: 10.1103/PhysRevA.38.364

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Thompson, S. M., Haas, H. L., and Ghwiler, B. H. (1992). Comparison of the actions of adenosine at pre and postsynaptic receptors in the rat hippocampus *in vitro*. *J. Physiol*. 451, 347–363.

Timofeev, I., Grenier, F., Bazhenov, M., Sejnowski, T. J., and Steriade, M. (2000). Origin of slow cortical oscillation in deafferented cortical slabs. *Cereb. Cortex* 10, 1185–1199. doi: 10.1093/cercor/10.12.1185

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Timofeev, I., Grenier, F., and Steriade, M. (2001). Disfacilitation and active inhibition in the neocortex during the natural sleep-wake cycle: an intracellular study. *Proc. Natl. Acad. Sci. U.S.A*. 98, 1924–1929. doi: 10.1073/pnas.98.4.1924

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wilson, C. (2008). Up and down states. *Scholarpedia* 3, 1410. doi: 10.4249/scholarpedia.1410

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Yang, X., Du, S., and Ma, J. (2004). Do earthquakes exhibit self-organized criticality? *Phys. Rev. Lett*. 92:228501. doi: 10.1103/PhysRevLett.92.228501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zapperi, S., Lauritsen, K. B., and Stanley, E. (1995). Self-organized branching process: a mean-field theory for avalanches. *Phys. Rev. Lett*. 75, 4071. doi: 10.1103/PhysRevLett.75.4071

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keywords: neuronal avalanches, oscillations, rat cortex, waiting times, neuronal networks, criticality

Citation: Lombardi F, Herrmann HJ, Plenz D and De Arcangelis L (2014) On the temporal organization of neuronal avalanches. *Front. Syst. Neurosci*. **8**:204. doi: 10.3389/fnsys.2014.00204

Received: 26 June 2014; Accepted: 01 October 2014;

Published online: 28 October 2014.

Edited by:

Paolo Massobrio, University of Genova, ItalyReviewed by:

James A. Roberts, Queensland Institute of Medical Research Berghofer Medical Research Institute, AustraliaLeonardo Paulo Maia, Universidade de São Paulo, Brazil

Copyright © 2014 Lombardi, Herrmann, Plenz and De Arcangelis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Fabrizio Lombardi, Institute of Computational Physics for Engineering Materials, ETH, Schafmattstr. 6, Zürich, 8093, Switzerland e-mail: fabrizio.lombardi@ifb.baug.ethz.ch