We model the neurons as leaky integrate-and-fire (LIF) units. The postsynaptic membrane potential of neuron i, when the neuron does not emit a spike, is given by the equation

where τ_m is the characteristic time of the membrane, C the membrane capacity, and I_i(t) the total current input to neuron i. The input is given by

where t_j are the spike times of neuron j, Q_ij is the total charge released at the synapse between neuron i and j, τ_s is the characteristic time of the synapse, t^i are the times of noise events releasing a random charge Q^i at some point of the membrane of neuron i, and the sum is extended to the spikes t_j and noise events t^i between the last spike t_i of neuron i, and the present time t. Defining J_ij = Q_ij/[C(1 − τ_s/τ_m)] and J^i=Q^i/[C(1−τs/τm)], we therefore have (Gerstner et al., 1993; Gerstner and Kistler, 2002)

where ϵ(t) = e^−t/τ_m − e^−t/τ_s. When the potential V_i(t) reaches the threshold value Θ_i, the neuron i emits a spike, and its potential is reset to the base value V_i = 0. In the present paper we set the same threshold Θ_i ≡ Θ for all the neurons, τ_m = 10 ms, τ_s = 5 ms, we extract the times t^i of noise events from a Poissonian distribution with a rate ρ = 1 ms⁻¹ for each neuron, and extract Ĵ_i from a Gaussian distribution with zero mean and standard deviation αρ∑jJij2. The constant α, which has the dimension of a rate, sets the “noise level” of the network.

The synapse strengths J_ij are held fixed during the simulation (no short term plasticity). They are set at the beginning with a “learning procedure” (Scarpetta et al., 2001, 2002; Scarpetta and Marinaro, 2005; Scarpetta and Giacco, 2012; Scarpetta et al., 2013), inspired to spike time dependent plasticity (STDP) (Markram et al., 1997, 2011). During this initial “learning procedure,” we store P patterns in the network connections. A pattern μ = 1, …, P is a phase-coded spike train of period T^μ, with one spike per neuron and per cycle, where the activity of neuron i is given by

The times t^μ_i are given by (ϕiμ2π)Tμ, where ϕ^μ_i are phases chosen randomly in [0, 2π) and then kept fixed, that give the “order of spiking” of neurons within pattern μ. Therefore, during the initial learning procedure, the network is forced to replay pattern μ, and the connections evolve due to STDP, so that in the interval [−T, 0] the change in the connection J_ij is given by

where H_i is a constant depending on the postsynaptic neuron i that sets the strength of the connections, TμT is a normalization factor, x_j(t) is the activity of the presynaptic neuron at time t, and x_i(t) the activity of the postsynaptic one. In STDP, the learning window A(τ) is the measure of the strength of the synaptic change when a time delay τ occurs between pre and post-synaptic spikes. The window A(τ) is the one introduced and motivated by Abarbanel et al. (2002), with the same parameters used in Abarbanel et al. (2002) to fit the experimental data of Bi and Poo (1998), see Figure 1. This function satisfies the balance condition ∫^∞_−∞ A(τ) dτ = 0. Notably, when A(τ) is used in Equation (5) to learn phase-coded patterns with uniformly distributed phases, then the balance condition assures that the sum of the connections on the single neuron ∑_j J_ij is of order 1/N, and therefore, it assures a balance between excitation and inhibition (Scarpetta et al., 2010). Note that, as we are studying a network of excitatory neurons, the negative connections have to be thought as connections mediated by fast inhibitory interneurons. When multiple phase-coded patterns are stored, the learned connections are simply the sum of the contributions from individual patterns, namely,

Throughout the paper we use a number of neurons N = 3000, a period T^μ = 333 ms, and a number of patterns P = 2. Moreover we use two values for the strengths H_i of the connections in Equation (5). A value H₀ for “normal” neurons, and a value H₁ = 3H₀ for “leader neurons,” that are chosen for each pattern μ as a fraction of 3% of the neurons that have consecutive phases, for example the lowest phases in the interval [0, 2π). Note that the values of H₀ are expressed in units of the threshold Θ of the neurons. The role of these few “leader” neurons, with higher incoming connection strenghts, is that of collect and amplify activity initiated by noise, and give rise to a cue able to initiate the short collective replay.

After the learning procedure, we perform a pruning procedure, by which only a fraction of the N(N − 1) connections J_ij survives. Namely, for each neuron i, we take all the incoming connections J_ij, and separately consider the positive (excitatory), and the negative (inhibitory) ones. As for the positive ones, we delete a fixed fraction f⁺_prune of them that have the lowest value. Then, we delete a fraction f^−,i_prune, that can depend on the neuron i, of the negative connections that have the lowest (absolute) value, choosing f^−,i_prune so that the sum of the incoming connections to neuron i at the end is as close as possible to zero. Throughout the paper we use f⁺_prune = 70%. As a consequence, at the end of the pruning process, about 12% of the N(N − 1) connections survived as positive connections, and 27% as negative connections, with statistical fluctuations of order of 1/N. After the learning and pruning procedure is applied, the dynamics of the network is studied with the connections J_ij fixed, that is we do not apply STDP nor short term depression.

In Figure 2A we show the average spiking rate in Hz per neuron, as a function of the noise level α and the connection strengths H₀. While the average rate increases continuously as either of the parameters is increased, the distribution of the rates in a finite interval of time changes qualitatively. We bin the time in 1 ms intervals, evaluate the rate in Hz per neuron for every interval, and compute the distribution of the rates. For low connection strength, or low noise, the distribution is nearly exponential (see Figure 3A), with an average rate lower than 2 Hz. For high connection strength and high noise, the distribution is nearly Gaussian (see Figure 3C), with an average rate higher than 13 Hz. In an intermediate region the distribution is bimodal (see Figure 3B), and shows both peaks, one exponential at low rates, and one Gaussian at high rates with a minimum in the distribution. This three different regimes are shown in Figure 2B with different colors. The intermediate bimodal regime resembles the phase coexistence observed in a first order equilibrium phase-transition, even though in our case the transition is a non-equilibrium one.

The qualitative difference in the distribution of the spiking rates, corresponds to a different dynamical behavior. At low rates, when the distribution is nearly exponential with a maximum at zero rate, the dynamical behavior is dominated by noise. The potential of neurons is governed by a Ornstein–Uhlenbeck process, and with some probability crosses the threshold giving rise to a spike, that is not able, however, to generate a spreading activity in the network (Figure 4A). On the other hand, in the high rate regime, the noise triggers the replay of one of the patterns encoded in the network (Figure 4C). In this case, once the replay of the pattern has started, the noise is not able to stop it, so that the replay is permanent². The intermediate regime, corresponding to a bimodal distribution of the rates, is shown in Figure 4B. In this case the noise is able to start the replay of a pattern, but also to stop it, so that the activity is intermittent, and resembles the experimentally observed alternation of up and down states. The firing rate is high in the “up” state (during short replays) and is low during “down” states, therefore the distribution of the rates is bimodal. The corresponding region, shown with green dots in Figure 2B, separates the regimes of permanent collective replay of spatiotemporal pattern (red) and the region of quiescence with low activity (yellow). Such non-equilibrium phase transition has been recently studied in a similar model (Scarpetta and de Candia, 2013, 2014) showing that in the region where the order parameter, which measures the similarity between spontaneous dynamics and the stored dynamic patterns, passes from zero to one, the fluctuations of the order parameter are maximized.

Note that the replay of a pattern appears in the raster plot of Figure 4 as a sawtooth, when the neurons are sorted on the vertical axes in the order of the pattern that is being replayed. On the other hand, it appears as completely random when the neurons are sorted in another way, for example in the order of a pattern that is not being replayed. The alternation of states seen in the raster plot in the bimodal region (Figure 4B) resembles the alternation of up and down states observed to occur spontaneously. Notably, as reviewed in Luczak and MacLean (2012), there are experimental evidences that during the up states often neurons activate in a surprisingly similar sequential order, reproducing default spatiotemporal patterns.

In the intermediate regime, where the network alternates between up and down states, we observe that inside the periods of high firing rate (up states), at a finer level, the activity is made of a series of cascades or “avalanches,” separated by short drops in the rate. Cortical activity cascades that follow precise power laws, i.e., neural avalanches, have been observed experimentally during spontaneous cortical activity in vitro and in vivo (Plenz, 2012, and references therein).

Experimentally neural avalanches are defined in terms of local field potential recorded at electrodes, that average the activity of many neurons. In our model, we have to distinguish between the few spikes generated by noise, that we want to characterize as no activity, and the spikes generated when a collective pattern is replayed, that represent instead an activity in the network. Due to the separation we have seen on the global spiking rates, with rates lower than 2 Hz corresponding to no activity, and rates larger than 13 Hz representing the collective replay of a pattern, we identify “avalanches” as consecutive time bins with a rate higher than a threshold R_min = 7 Hz. Successive time bins are concatenated until an empty bin (rate lower then R_min) is reached, at which the concatenation process stops.

We define the size of an avalanche as the total number of spikes, that is the integral of the rates over the avalanche duration. In Figure 5 we show the distribution of the sizes (A) and durations (C) of the avalanches for H₀/Θ = 0.221, α = 0.06, that corresponds to point B in Figure 2B. Note that the distributions are well described by power laws, with exponent 3/2 for the sizes and 2 for the durations, as experimentally observed (Plenz, 2012). Such a behavior is quite robust, and is observed generically in the region of the non-equilibrium transition between the replay and non-replay of spatiotemporal patterns (Scarpetta and de Candia, 2013). It does not depend on the precise value of the R_min chosen. Indeed, as shown in Figure 5B, the size distribution follows approximately the same power law in a range of R_min from R_min = 5 Hz to R_min = 9 Hz.

We have computed the distribution P(Δt) of the waiting times between successive avalanches. In Figure 6A we show the distribution for synaptic strength H₀/Θ = 0.214, 0.221, and 0.228, and noise α = 0.06 ms⁻¹, in the region where the rate distribution is bimodal. The middle (red) curve at H₀/Θ = 0.221 corresponds to point B in Figure 2B. The distribution presents a regime between 10 and 50 ms characterized by a power law with exponent −3, preceded by a regime with a lower slope. For times larger than 50 ms, the distribution shows a broad plateau, that is longer the lower the noise, or the strength of the connections.

A power law regime in the waiting times between avalanches has been observed also experimentally, for example in freely behaving rats (Ribeiro et al., 2010), or in cortical slices (Lombardi et al., 2012, in preparation). The second regime, corresponding to large waiting times, is also observed in Lombardi et al. (2012, in preparation).

The power law regime corresponds to waiting times between successive avalanches within the same up state. The power law in the distribution indicates temporal correlation, i.e., that consecutive avalanches belonging to the same up state are correlated. Indeed in our model the up state is the result of the system falling in one of the many metastable spatiotemporal pattern attractors, corresponding to a collective replay activity. The large bump at long times in the waiting time distribution is related on the other hand to down states, that is intervals in which the network does not replay any of the encoded patterns.

We have plotted data of Figure 6A also in an alternative way. While P(Δt)δt is the probability of observing a waiting time between Δt and Δt + δt, in Figure 6B we plot P(Δt), where P(Δt)δλ is the probability of observing a waiting time between Δt and Δt(1 + δλ). Note that P(Δt) = P(Δt)Δt. With this alternative definition, the distribution becomes non-monotonic, with a pronounced maximum at high values of Δt, and the exponent of the initial power law becomes −2.

In Figure 7 we show the distribution P(Δt) in the region of high activity, marked in red in Figure 2B. In this case the replay of the pattern becomes continuous, and therefore the plateau at long times, in the distribution of waiting times, disappears. The distribution therefore shows only the power law regime, as shown in Figure 7, corresponding to the point in phase space marked with letter C in Figure 2B. On the other hand, in the region of low activity, marked in yellow in Figure 2B, avalanches become very sparse, so that the distribution of waiting times is different from zero only for very long times, and in this case the initial power law regime disappears.

As in Lombardi et al. (2012), we define the up states as periods of high activity characterized by the concatenation of consecutive avalanches with waiting times lower than T_max = 50 ms, the maximum time falling inside the power law regime of waiting times. Successive avalanches are concatenated until a waiting time larger than T_max is reached, at which the concatenation process stops. Similarly, down states are defined as a concatenation of waiting times larger than T_max. An isolated avalanche preceded and followed by a waiting time larger than T_max does not stop the down state.

In Figure 8 we show the distribution of the durations of down and up states, for the same parameters of Figure 6. The distribution of durations of down states (Figure 8A) is well fitted by an exponential (continuous lines in the figure), showing that the transition from down to up states is controlled by a Poissonian probability, due to the noise focusing that triggers a replay of one of the patterns encoded in the network. On the other hand, the distribution of durations of up states cannot be fitted by an exponential as well as the one of down states, except for a narrow interval at large times. Indeed, one observes an excess of durations around 100 ms, that corresponds in the model to the duration of one period of the replayed pattern. Moreover, when one approaches the region of parameter space where the replay of the patterns becomes continuous, the distributions show significant deviations from the exponential also for large times, and could be better fitted by a stretched exponential. This is apparent going from blue to red and green curves in Figure 8B, that are all in the bimodal region, but get closer and closer to the region of self-sustained replay. Note that, when one goes deep inside the region of self-sustained replay, no down states are observed in practice, and up states last for a time of the order of the experimental time.

It has been conjectured that the alternation between up and down states depends crucially on the short term synaptic depression (STSD). As we have shown, in our model this instability between up and down states is present even in absence of short term depression, and is due instead to the particular structure of connections, that are far from being random. Such structure determines in the network a large transition region of phase space, where there is a co-presence of both dynamical attractor states, corresponding to the replay of the patterns encoded, and the attractor corresponding to quiescence of the network.

However, as short term depression is present in real synapses in the brain, we show here that it does not invalidate the behavior displayed by the model considered here, but changes only the parameters, such as the strength of connections, where the transition region appears. We have added STSD in the model, implementing a dynamics on the connections J_ij according to the equation

where J⁰_ij are the connections given by Equation (6), and τ_r = 10 ms is the recovery time of synapses. Moreover, we depress J_ij by a factor f_stsd = 0.5, every time the presynaptic neuron j fires a spike.

In Figure 9A, we show the distribution of the rates, evaluated as in Figure 3, in absence of STSD for H₀/Θ = 0.221 (point B in Figure 2), and in presence of STSD for H₀/Θ = 0.221 and 0.231. Note that, for the same synaptic strength at rest, the effect of STSD is to lower the fraction of time in which the network is in the up state. However, for a slightly higher H₀/Θ, the distribution is very similar to the one without STSD. In Figure 9B, we show the raster plot in the case of H₀/Θ = 0.231 and τ_r = 10 ms, showing a behavior very similar to the one displayed in Figure 4B, with the alternation of up and down states, and of different patterns replayed in the up states. Therefore, the effect of STSD is to slightly shift the region of the phase space in which the transition is observed.

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.