Stochastic oscillations and dragon king avalanches in self-organized quasi-critical systems

In the last decade, several models with network adaptive mechanisms (link deletion-creation, dynamic synapses, dynamic gains) have been proposed as examples of self-organized criticality (SOC) to explain neuronal avalanches. However, all these systems present stochastic oscillations hovering around the critical region that are incompatible with standard SOC. Here we make a linear stability analysis of the mean field fixed points of two self-organized quasi-critical systems: a fully connected network of discrete time stochastic spiking neurons with firing rate adaptation produced by dynamic neuronal gains and an excitable cellular automata with depressing synapses. We find that the fixed point corresponds to a stable focus that loses stability at criticality. We argue that when this focus is close to become indifferent, demographic noise can elicit stochastic oscillations that frequently fall into the absorbing state. This mechanism interrupts the oscillations, producing both power law avalanches and dragon king events, which appear as bands of synchronized firings in raster plots. Our approach differs from standard SOC models in that it predicts the coexistence of these different types of neuronal activity.

Notice that (ρ 0 = 0, Γ) is not a fixed point. The single fixed point is: To perform the fixed point stability analysis we calculate partial derivatives evaluated at the fixed point (denoted by ] * ): The eigenvalues of the system Jacobian matrix at the fixed point (ρ * , Γ * ) are determined by So, the characteristic equation is: That can be written as: with solutions: This equation has no meaning for τ < 2 since this would give ρ * = 1/τ > 1/2. When 2 < τ < 2 + √ 2 we have ∆ > 0, λ real and |λ | < 1, that is, the fixed point (ρ * , Γ * ) is a stable node. However, when τ > 2 + √ 2 = 3.4142 . . ., we have ∆ < 1 and two complex roots λ ± with |λ | = √ λ + λ − < 1, that is, we have a stable focus with: 1/6 which leads to Eq. (13) of the main Manuscript For small amplitude (harmonic) oscillations, the frequency is given by the angle of the eigenvalues written in polar coordinates. We can write it in terms of the determinant D and trace T of the Jacobian matrix J as ω ≡ |ω ± | = arctan 4D Case µ > 0 In the case µ > 0, instead of a single equation to describe ρ[t], we need some recurrence equations to describe the evolution of firing densities. As described in 1 and 2 , the density of firing neurons ρ[t] can be computed from the probability density p[t](V ) of potentials at time t: Suppose that at time t 0 every neuron in the network has already fired at least once. Then, we can describe ρ[t] for t > T 0 with respect to the evolution of densities of neurons according their firing age. In this case, Eq. (17) becomes: Where η k [t] is the proportion of neurons at time t > T 0 that fired k time steps ago, which all have potential U k [t]. The proportion of neurons with firing age k at t corresponds to fraction of the population with firing age k − 1 that did not fire at time t − 1: And from Eq. (2) of the main Manuscript we know that the potential evolution of these populations follow: The normalization condition ∑ ∞ k=0 η k [t] = 1 must be included explicitly. Although, we cannot obtain a simple 2d map for ρ and Γ as is the case for µ = 0, we can instead, solve Eqs (19-20) numerically for any µ at every time step as Γ evolves concomitantly according to Eq. (2).
In 2 we estimate from Eqs (19-20), approximating for small ρ, that in the stationary state: where (2) the value of ρ for the non trivial fixed point. Using its value in 21, we obtain that, for µ > 0: Figure 1 shows the numerical iteration for the simplified adaptive model with one parameter in two cases: for µ = 0 (Fig. 1a), which corresponds to the 2D map in Eqs (1 and 2), and for µ = 0.5 (Fig.˜1b) by the iterations of Eqs (2, 19 and 20).

Stability analysis of the stochastic neuron model with LHG gain dynamics
The MF map is given by: We have two fixed points. The absorbing state fixed point: and the non trivial fixed point (ρ * , Γ * ): Partial derivatives of F and G with respect to ρ and Γ are:

3/6
Evaluating at the trivial fixed point, we get: The eigenvalues are found from (a − λ )(d − λ ) − bc = 0, or: with solution: Notice that ∆ ≥ 0 so the eigenvalues are always real: the fixed point is an attractive or repulsive node. The final solution of Eq. 37 is: This means that the fixed point This analysis show that we can have a subcritical state (0, A), an attractive node, if we use A < 1. On the other hand, if we use A > 1, the absorbing state looses its stability and the non trivial fixed point (ρ * , Γ * ) is the unique stable solution. For A = 1, the point (ρ 0 = 0, Γ 0 = 1) is an indifferent node.
Evaluating the derivatives at the non trivial fixed point, we get: As before, the eigenvalues of the Jacobian matrix at the fixed point (ρ * , Γ * ) are determined by Eq. (8). After some algebra, we get: For large τ, we have:

4/6
which means that, for large τ, the map is very close to a Neimark-Sacker-like critical point. As before, the focus can be excited by fluctuations, generating the SO.
The frequency for small amplitude oscillations is ω ≡ |ω ± | = arctan 4D T 2 − 1, where D and T are respectively the determinant and trace of the Jacobian matrix. When τ → ∞, which is similar to the τ −1/2 behavior of the simple gain dynamics, Eq. (16).

Stability analysis of the Cellular Automata model with LHG synapses
Starting from Eqs (24 and 25) of the main Manuscript, the 2d MF map close to the stationary state is: As before, the stability of the map is given by the eigenvalues of the Jacobian matrix J(ρ, σ ) at the fixed points. The matrix elements are: Thus, we can plug in the solutions of Eqs (26 and 27) of the main Manuscript in Eqs (50-53) and find the eigenvalues of J(ρ, σ ).
Similarly to the previous model, the MF map of the CA with LHG synapses has two fixed points. The first one is the absorbing state fixed point and is given by The second one is nontrivial and must be solved numerically. We can compute the elements of the Jacobian matrix at the absorbing state fixed point which is the same result of Eqs (31-34). This matrix has two eigenvalues: 1 − 1/τ and A. The eigenvalue 1 − τ corresponds to the σ axis, which is a stable direction. The eigenvalue A corresponds to a unstable (stable) direction for A > 1 (A < 1).
If we assume that we are close to criticality, where ρ * is small, we can find an approximate solution for J(ρ, σ ) at the nontrivial fixed point. Expanding the binomial in Eq. (24) of the main Manuscript to first order and solving for ρ * , we find: