Self-Organized Supercriticality and Oscillations in Networks of Stochastic Spiking Neurons

Networks of stochastic spiking neurons are interesting models in the area of Theoretical Neuroscience, presenting both continuous and discontinuous phase transitions. Here we study fully connected networks analytically, numerically and by computational simulations. The neurons have dynamic gains that enable the network to converge to a stationary slightly supercritical state (self-organized supercriticality or SOSC) in the presence of the continuous transition. We show that SOSC, which presents power laws for neuronal avalanches plus some large events, is robust as a function of the main parameter of the neuronal gain dynamics. We discuss the possible applications of the idea of SOSC to biological phenomena like epilepsy and dragon king avalanches. We also find that neuronal gains can produce collective oscillations that coexists with neuronal avalanches, with frequencies compatible with characteristic brain rhythms.


Introduction
Neuronal network models are extended dynamical systems that may present different collective behaviors or phases characterized by order parameters. The separation regions between phases can be described as bifurcations in the order parameters, or phase transitions. In several models of neuronal activity, the relevant phase change is a continuous transition from an absorbing silent state to an active state [1][2][3]. In such continuous transition we have a critical point (in general, a critical surface) where concepts of universality classes and critical exponents (among others) are valid. At criticality, we observe avalanches of activity described by power laws for their size and duration. Also, the avalanche profile shows fractal scaling. Since the landmark findings of Beggs and Plenz in 2003 [2], these behaviors have been reported also in biological networks, see reviews [4][5][6].
The motivation for the idea that criticality is important to understand neuronal activity is not only empirical. Several works have shown that there are advantages for a network to operate at the critical state [3,[7][8][9]. However, it is not clear how biological networks tune themselves to the critical region that, in a parametric space with P free parameters, has a maximum of P − 1 dimensions.
An important idea discussed in several papers is that, since criticality depends on the strength of the synapses (links) between the neurons, a homeostatic mechanism for dynamic synapses tunes the network toward the critical region. There are two main paradigms: self-organization of Hebbian synapses [10][11][12][13][14][15], and self-organization of dynamic synapses [16][17][18][19][20] following Tsodyks and Markran [21,22]. With these synaptic mechanisms it is possible to achieve, or at least to approximate, a self-organized critical (SOC) state.
With a different approach, we have shown recently that dynamic neuronal gains, biophysically linked to firing dependent excitability of the Axonal Initial Segment (AIS) [23][24][25][26], can also lead to self-organized criticality [27]. This new mechanism is simpler than dynamic synapses because, for biological networks with N neurons, we have of the order of 10 4 N synapses [28] -and thus 10 4 N dynamic equations -but only N equations for the neuronal gains.
It has also been observed in Brochini et al. [27] that, for achieving exact SOC, all papers in the literature used a time scale τ for synaptic recovery proportional to N, with N → ∞. The use of this non-local information (N), and a diverging recovery time τ, is not plausible biologically. A similar time scale τ is present in the neuronal gain recovery. When we use a biological range for τ that does not scale with N, we observe that the network turns out (slightly) supercritical, a phenomenon which we called self-organized supercriticality (SOSC). That is, both dynamic synapses and dynamic gains with fixed τ, which seems to be the reasonable biological assumption, presents SOSC instead of SOC.
Here we report for the first time an extensive study of the neuronal gain mechanism and SOSC. First, we present new mean-field results for phase transitions in a fully connected model of integrate-and-fire stochastic neurons with fixed gains. We find both continuous and discontinuous phase transitions. Then, we introduce a simplified gain dynamics depending only on the τ parameter, which also has a simple mean-field solution in the case of the continuous transition and presents SOSC. We compare this solution with extensive simulations for different system sizes N and values for τ. Surprisingly, we found collective oscillations produced by the gain dynamics that coexist with neuronal avalanches.
If a presynaptic neuron j fires at discrete time t, then X j [t] = 1. This event increments by W ij the potential of every postsynaptic neuron i that have not fired at time t. The potential of a non-firing neuron may also integrate an external stimulus I i [t]. Apart from these increments, the potential of a non-firing neuron decays at each time step towards zero by a factor µ ∈ [0, 1], which models the effect of a current leakage.
The neuron membrane potentials evolve as: This is a special case of the general model from [31] with the filter function g(t − t s ) = µ t−t s , where t s is the time of the last firing of neuron i [27]. In contrast to standard IF neurons, the firing is not deterministic above a threshold, but stochastic. We also have X i [t + 1] = 0 if X i [t] = 1 (refractory period of one time step).
The firing function 0 ≤ Φ(V) ≤ 1 is sigmoidal, that is, monotonically increasing. We also assume that Φ(V) is zero up to some threshold potential V T . If Φ is the shifted Heaviside step function Φ(V) = Θ(V − V T ), we have a deterministic discrete-time leaky integrate-and-fire (LIF) neuron. Any other choice for Φ(V) gives a stochastic neuron.
In Brochini et al. [27] we have studied a linear saturating function with neuronal gain Γ similar to that used in [32]. Here we study the so called rational function that does not have a saturating potential, see Fig. 1a: Notice that we recover the deterministic LIF model The use of the rational instead of the linear saturating function is convenient and gives some theoretical advantages, for example to avoid the anomalous cycles-2 observed in [27].

Mean-field calculations
The network's activity is measured by the fraction ρ[t] of firing neurons (or density of active sites): The density of active neurons ρ[t] can be computed from the probability density p[t](V) of potentials at time t: where p[t](V) dV is the fraction of neurons with potential in the range [V, V + dV] at time t. Neurons that fire between t and t + 1 have their potential reset to zero. They contribute to p[t + 1](V) a Dirac impulse at potential V = 0, with amplitude ρ[t] given by equation (4). The potentials of all neurons also evolve according to equation (1). This process modifies p[t](V) also for V = 0.
In the mean-field limit, we assume that the synaptic weights W ij follow a distribution with average W = W ij and finite variance. By disregarding correlations, the term in Eq. (1) corresponding to the sum of all presynaptic inputs simplifies to Wρ[t].
If the external input is constant, I i [t] = I, a stationary state is achieved, which depends only on the average synaptic weight W, the leakage parameter µ and the parameters that define the function Φ(V), that is, Γ and V T . In Brochini et al. [27] it is shown that the stationary p(V) is composed of delta peaks with height η k situated at voltages U k given by: for all k ≥ 1. Here U k corresponds to the potential value of the population of neurons that have firing age k. The firing age is the amount of time steps since the neuron fired for the last time. The normalization condition ∑ ∞ k=0 η k = 1 must be included explicitly. Equations (6)(7)(8) can be solved numerically for any firing function Φ, so this result is very general.

Phase transitions for the rational Φ(V)
In terms of non-equilibrium statistical physics, ρ is the order parameter, I is an uniform external field and Γ and W are the main control parameters. The activity ρ also depends on V T and µ. 4.1.1. The case with µ > 0, I = 0, V T = 0 By using Eqs. (5)-(8), we obtain numerically ρ(W, Γ) for several values of µ > 0, for the case with I = 0, V T = 0 (Fig. 1b). Only the first 100 peaks (U k , η k ) were considered, since, for the given µ and Φ, there was no significant probability density beyond that point. The same numerical method can be used for studying the cases I = 0, V T = 0.
We also obtained an analytic approximation (see Appendix) for small ρ: where Γ C = (1 − µ)/W defines the critical line and ∆ Γ = (Γ − Γ C )/Γ is the reduced control parameter. So, the critical exponent for the order parameter near criticality is β = 1, characteristic of the mean-field directed percolation (DP) universality class [37]. We also compare Eq. (9) with the numerical results for ρ(Γ, µ) in Fig. 1c.

Analytic results for µ = 0
In the case µ = 0 it is possible to do a simple mean-field analysis valid for N → ∞. This case is illustrative because it presents all phase transitions that occur with µ > 0.
When µ = 0 and I i [t] = I (uniform constant input), the stationary density p(V) consists of only two Dirac peaks at potentials U 0 = 0 and U 1 = I + Wρ. Eq. (8) simplifies to: since η 0 = ρ and η 1 = 1 − ρ. By inserting the function Eq. (2) in Equation (10), and remembering that Φ(0) = 0, we get: with solutions: 4.1.3. The case with I = 0, V T = 0: continuous transition where the phase transition line is and the critical exponent is β = 1. This corresponds to a standard mean-field continuous (second order) absorbing state phase transition ( Fig. 2a and 2b with V T = I = 0).

The case with I < V T : discontinuous transition
For I < V T , we have discontinuous (first order) phase transitions when ∆ = 0, see Eq. (13): which, after some algebra, leads to the phase transition lines: which have the correct limit, Eq. (16), when V T → 0, I → 0. The transition discontinuity is: In Fig. 2a we show examples of the phase transitions, which occurs when the unstable point ρ − collapses with the stable point ρ + . It is important to notice that, for any V T > 0, the unstable point never touches the absorbing point ρ 0 = 0, so the zero solution is always stable. Only for the case V T = 0 the solution ρ 0 looses stability and ρ + is the unique solution above the critical line Γ C = 1/W C . Fig. 2b gives the phase diagram Γ × W for some values of V T with I = 0, see Eq. (19). Finally, we give the phase diagram for the variables ΓW versus Γ(V − I), see Fig. 3 and Eq. (18).
This line is a first order phase transition which terminates at the second order critical point 4.2. Self-organized supercriticality (SOSC) through dynamic gains with µ = 0, If we fine tune the model to some point in the critical line Γ C = 1/W, we can observe perfect neuronal avalanches with size distribution P S (s) ∝ s −3/2 and duration distribution P D (d) ∝ d −2 [27]. As expected, these are mean-field exponents fully compatible with the experimental results [2,6].
This fine tuning, however, is not plausible biologically. What we need is some homeostatic mechanism that makes the critical region an attractor of some self-organization dynamics. In the literature, a well studied mechanism is dynamic synapses W ij [t] [16][17][18]. For example, in discrete time [19,20]: where τ is a synaptic recovery time, A is an asymptotic value and u ∈ [0, 1] is the fraction of depletion of neurotransmitter vesicles when the presynaptic neuron fires. In Brochini et al. [27], we proposed a new self-organization mechanism based in dynamic neuronal gains Γ i [t] while keeping the synapses W ij fixed [27]. The idea is to create a feedback loop based only in the local activity X i [t] of the neuron, reducing the gain when the neuron fires, and recovering slowly after that. The biological motivation for dynamic gains is spike frequency adaptation, a well-known phenomenon that depends on the decrease (and recovery) of sodium ion channels density at the axon initial segment (AIS) when the neuron fires [24,25].
The dynamics for the neuronal gains studied in [27] has a form similar to that used in [16,[18][19][20] for synapses: The advantage of neuronal gains is that now we have only N dynamical equations (notice the term X i [t] that refers to the activity of the postsynaptic neuron, not of the presynaptic one as in Eq. (21)). For dynamic synapses, we need to simulate N(N − 1) equations for the fully connected graph model, and 10 4 N for a biologically realistic network, and this is computationally very costly for large N.
A problem with this dynamics, however, also present in dynamic synapses, is that we have a three dimensional parameter space (τ ∈ [1, ∞], A ∈ [1/W, ∞], u ∈ [0, 1]) that must be fully explored to characterize the stationary value Γ * (τ, A, u, N). Here we propose a new simplified dynamics with a single free parameter, the gain recovery time τ: The self-organization mechanism can be viewed in Fig. 4. So, we reduce our parametric study to determine the curves Γ * (1/τ, 1/N), see Figs. 5a and 5b. The fluctuations measured by the standard deviation SD of the Γ[t] time series, after the transient, diminishes for increasing τ (Fig. 5c) and probably goes to zero for τ → ∞, in accord to Campos et al. [20]. However, in contrast to this idealized τ → ∞ limit, as discussed in the references [18,20], the fluctuations do not converge to zero for finite τ in the thermodynamic limit N → ∞ (see Fig. 5d). This occurs because, for low τ, the adaptation mechanism produces oscillations of Γ[t] around the value Γ * (τ) We can do a mean-field analysis of Eq. (23) to find the value Γ * (τ). Denote the average gain as Averaging over the sites, we have : A solution is Γ * = 0, but this is unstable, see Eq. (23). Another solution is obtained by inserting Eq. (15), ρ * = (Γ * − Γ C )/(2Γ * ), in Eq. (25): Notice that this is valid only when Eq. (15) is valid, that is, for Γ * ≥ 1/W. Also, Eq. (26) presumes that ρ * is a stable fixed point, which can not be true for some interval of values of τ, see below.
A first order approximation leads to: This mean-field calculation shows that, if τ → ∞, we obtain an exact SOC state Γ * → Γ C . Or, for finite networks, it would be required a scaling τ = O(N a ) with an exponent a > 0, as done previously for dynamic synapses [16,[18][19][20]. However, this scaling for τ cannot be justified biologically. So, biology requires a finite recovery time τ which always leads to supercriticality, see Eq. (26) or (27). This supercriticality is self-organized in the sense that it is achieved and maintained by the gain dynamics Eq. (23). We call this phenomena self-organized supercriticality (SOSC).
The deviation from criticality can be small. For example, if τ = 1000 ms (assuming 1 time step equals 1 ms in the model): Even a more conservative value τ = 100 ms gives Γ * ≈ 1.02 Γ C . Although not perfect SOC [5], this result is sufficient to explain a power law with exponent 3/2 for small (s < 1000) neuronal avalanches plus a supercritical bump (Fig. 6). By using Eq. (26) in Eq. (15) we also obtain: showing that the network presents supercritical activity for any finite τ. This result, however, is valid only in the infinite size limit. For finite size networks, fluctuations interrupt this constant ρ * = 1/τ activity, leading the system to the absorbing state and defining the end of the avalanches. Simulations reveal that these fixed points (Γ * , ρ * ) correspond only to mean values around which both Γ[t] and ρ[t] oscillate, see Figs. 7a, 7b and 7c. These global oscillations are unexpected since the model has been devised to produce avalanches, not oscillations. The finite size fluctuations and oscillations drive the network to the absorbing zero state, generating the avalanches. What we see in the histogram of Fig. 6 is a combination of power law avalanches in some range plus large events (superavalanches or Dragon-kings due to the Γ[t] oscillations).

Discussion
We examined a network of stochastic spiking neurons with a rational firing function Φ that has not been studied previously. We obtained numeric and analytic results that show the presence of continuous and discontinuous absorbing phase transitions. Classic SOC is possible only at the continuous transition, which means that we need to use zero firing thresholds (V T = 0). In some sense, this is a kind of fine tuning of the Φ function, but not the usual one where the synaptic strength W and the neuronal gain Γ are the main control parameters.
The presence of a well behaved absorbing phase transition in the directed percolation class enables the use of an homeostatic mechanism for the neuronal gains that tunes the network to the critical region. The dynamics on the gains is biologically plausible, and can be related to a decrease and recovery, due to the neuron activity, of the firing probability at the Axon Initial Segment (AIS) [23]. Our dynamic Γ i [t] mimics the well known phenomenon of spike frequency adaptation [24,25] and is a one-parameter simplification of the three-parameter dynamics studied by us in [27]. We observe that this gain dynamics is equivalent to approach the critical line with fixed W and variable Γ[t], that is, performing vertical movements in Fig. 2b. Previous literature approaches the critical point W C by dynamic [16][17][18][19][20] or Hebbian [10][11][12][13][14][15] synapses. This corresponds to fix Γ and allow changes in W ij [t] along the horizontal axis, see Fig. 2b.
The two homeostatic strategies are similar, but we stress that we have only N equations for the gains Γ i [t] instead of N(N − 1) equations for the synapses W ij [t], so that our approach imply a huge computational advantage. Indeed, previous literature as [13,16] report system sizes on the range of N = 1, 000 − 4, 000, to be compared to our maximal size of N = 160, 000.
We found that the fixed point Γ * predicted by a mean-field calculation is not exactly critical, but instead supercritical, and that the distance from criticality depends on the gain recovery time τ. Previous claims about achieving exact SOC by using dynamic synapses are based on the erroneous assumption that we can use a synaptic recovery time τ ∝ N a → ∞ [16,[18][19][20]. But if we use a finite τ, which is not only plausible but biologically necessary, we obtain SOSC, not SOC [20,27]. However, we found that for large but plausible values of τ, the system is only slightly supercritical and presents power law avalanches (plus small supercritical bumps) compatible with the biological data.
SOSC enables us to explore supercritical networks that are robust, that is, the stationary state, with or without oscillations, is achieved from any initial condition and recovers from perturbations. So, the question now is: are there self-organized supercritical (SOSC) oscillating neuronal networks in the brain? A first evidence would be a supercritical bump in the P(s) distributions. Indeed, we found several papers where such bumps seem to be present, see for example the first plot in Fig. 2 of Friedman et al. [38] and Fig. 4 of Scott et al. [39]. It seems to us that, since the main paradigm for neuronal avalanches is exact SOC, with pure power laws, it is possible that researchers report what is expected and do not comment or emphasize small supercritical bumps, even if they are present in their published data. Therefore, we suggest that experimental researchers reevaluate their data in search for small supercritical bumps. The presence of supercritical bumps can also be masked by the phenomenon of subsampling [40][41][42], so the analysis must be done with some care.
Supercriticality, in the form of the so called Dragon-king avalanches [13,43,44], has been conjectured to be at the basis of hyperexcitability in epilepsy [6,45,46]. Also, networks can be put artificially in the hyperexcitable state and show bimodal distributions P S (s) with large supercritical bumps [47]. The SOSC phenomenon seems to be a natural explanation for such hyperexcitability. In [47], the supercritical bumps are fitted by a supercritical branching process but are not explained in mechanistic terms as, in our case, due to different values of the biophysical τ recovery time.
The unexpected oscillations in Γ[t] around Γ * have amplitudes that depend on τ and vanish for large τ (Fig. 5d and Fig. 7a). These oscillations in Γ[t] induce oscillations in the activity ρ[t] (Fig. 7b  and 7c). In our model, the discrete time interval ∆t is postulated as describing the width of a spike, that is, ∆t ≈ 1 − 2 ms. From our simulation data, with these values for ∆t, we obtain frequencies f ≈ 0.5 − 16 Hz, depending on the τ value.
Interestingly, this frequency range includes Delta, Theta and Alpha rhythms. Coexistence of Theta waves and neuronal avalanches has been observed experimentally [48]. Also, some theoretical work recently discussed coexistence of oscillations and avalanches [49].
The presence of oscillations can mean that the fixed point (Γ * , ρ * ) is unstable below some biffurcation point τ b (even for N → ∞) or that it is stable but has a very small negative Lyapunov exponent, such that finite-size fluctuations drive Γ[t], ρ[t] away from equilibrium, producing excursions (oscillations) in the (Γ, ρ) plane. At this point, without further study, we can not decide what is the correct scenario. Notice that similar oscillations for W[t] were also observed for dynamic synapses [16,18,19], although these authors have not studied in detail such phenomenon.
Finally, from a conceptual point of view, the observed subcriticality in some of our simulations (see Fig. 5b and 5d) is less important than supercriticality (SOSC), because it is a finite-size effect for small N. Our largest networks have N = 160, 000, which is small compared to real biological networks that have at least one or two orders of magnitude more neurons.
Nevertheless, there is in the literature claims that subcritical states are present in certain experimental conditions [50][51][52]. How we can conciliate these findings? Here we offer an answer based on the findings of Priesemann et al. [52]. These authors found that, in order of explain in vivo experiments with awake animals, they need three ingredients: subsampling [40], increased input (violating the standard separation of scales of SOC models) and small subcriticality of the networks. If we increase the inputs in our network, by a Poisson process on the variable I[t] for example, the overall result is that the homeostatic mechanism turns out our network subcritical. This occurs because increased forced firing imply an overall depression of the gains Γ i [t] in Eq. (23), so that a new equilibrium is achieved with Γ * < Γ c .
Then, under external input like in awake animals, our adaptive networks turns out subcritical and returns to criticality or supercriticality for spontaneous activity without external input. We obtained preliminary simulation results confirming this scenario and a comprehensive study of the effect of external input shall be done in a next paper.

Materials and Methods
All numerical calculations were done by using MATLAB. Simulation codes were made in Fortran90.
In the study of neuronal avalanches, we simulate the evolution of finite networks with N neurons, uniform synaptic strengths W ij = W (W ii = 0) and Φ(V) rational with V T = 0. The avalanche statistics were obtained after the transient of the neuronal gains self-organization. A silent instant when X i [t] = 0 for all i defines the end of an avalanche. We start a new avalanche by forcing the firing of a single random neuron i, setting V i [t + 1] to a value high enough for the neuron spikes.

Conclusion
We have shown in this paper that dynamic neuronal gains lead naturally to self-organized supercriticality (SOSC) and not SOC. The same occurs with dynamic synapses [20]. So, we propose that neuronal avalanches are related to SOSC instead of exact SOC. This opens an opportunity for reevaluation of the accumulated experimental data.
SOSC suggests that neuronal tissues could be more prone to Dragon-king avalanches [43] and hyperexcitability than one would expect from simple power laws. This prediction of larger and increased instability due to supercriticality may be important for studies in epilepsy [13].
Finally, the emergence of oscillations coexisting with neuronal avalanches seems to unify in a single formalism two theoretical approaches and two different research communities: those that emphasize critical behavior and avalanches, and those that emphasize oscillations and synchronized activity.
In a future work, we intend to study with more care the mechanism that generates these oscillations and how to relate them to EEG data. In order to simulate more biological networks, we also intend to study the cases V T > 0, I > 0 and µ > 0.
Notice that ∑ ∞ k=1 η k−1 = ∑ ∞ k=0 η k = 1 by normalization. Using this fact and also Eq. (30) (using Φ(U 0 ) = 0), after some rearrangement we obtain: With respect to the last term, we use Eq. (34) and (30) to obtain which is valid for µ < 1. Here, the sum ρ ∑ ∞ k=1 η k−1 φ(U k−1 )µ k is composed of terms in ρ 3 than can be dismissed. Using this approximation in Eq. 39 we obtain two solutions. One is the absorbing state ρ = 0. The other solution is where we considered Γ c = (1 − µ)/W. Moreover, in the critical region we can approximate ΓW ≈ 1 − µ, leading to: We compare this analytical approximation with numerical solutions for ρ(ΓW, µ) near the critical point, see Fig. 1c.