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

TL;DR

This paper analyzes how stochastic oscillations in self-organized quasi-critical systems can produce power-law avalanches and rare dragon king events, explaining neuronal activity patterns that differ from standard self-organized criticality models.

Contribution

It provides a linear stability analysis of SOqC systems, revealing how demographic noise near a stable focus leads to diverse neuronal avalanches and dragon king events.

Findings

Stability analysis shows fixed points are stable foci losing stability at criticality.

Demographic noise near the focus induces stochastic oscillations that can produce avalanches.

Power law avalanches and dragon king events coexist due to noise-induced oscillations.

Abstract

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. This phenomenology has been called self-organized quasi-criticality (SOqC). Here we make a linear stability analysis of the mean field fixed points of two SOqC 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.

Equations

$$P(X_i[t]=1|V_i[t])=\Phi (V_i[t])=\frac{\Gamma V_i[t]}{1+\Gamma V_i[t]}\Theta (V_i[t]),$$

$$V_i[t+1]=\{\begin{array}{ll}\mu V_i[t]+I_i[t]+\frac{1}{N}{\sum }_{j=1}^N{W}_{ij}{X}_j[t]& \text{mbox}if{X}_i[t]=0,\\ 0& \text{mbox}if{X}_i[t]=1\end{array}$$

$$\rho [t]=\int \Phi (V)p(V)[t]dV,$$

$$\rho [t+1]=\frac{\Gamma W\rho [t](1-\rho [t])}{1+\Gamma W\rho [t]}\Theta (W\rho [t]).$$

$${\rho }^{\ast }=\frac{\Gamma W-1}{2\Gamma W},$$

$$W_{ij}[t+1]=W_{ij}[t]+\frac{\Delta t}{\tau }\left(A-W_{ij}[t]\right)-u{W}_{ij}[t]X_j[t].$$

$${\Gamma }_i[t+1]={\Gamma }_i[t]+\frac{1}{\tau }\left(A-{\Gamma }_i[t]\right)-u\Gamma [t]X_i[t],$$

$${\Gamma }_i[t+1]={\Gamma }_i[t]+\frac{1}{\tau }{\Gamma }_i[t]-{\Gamma }_i[t]X_i[t]=\left(1+\frac{1}{\tau }-X_i[t]\right){\Gamma }_i[t].$$

$$\rho [t+1]$$

$$\Gamma [t+1]$$

$${\rho }^{\ast }=\frac{1}{\tau },{\Gamma }^{\ast }=\frac{{\Gamma }_c}{1-2/\tau },$$

$${\rho }^{\ast }=\frac{{\Gamma }^{\ast }W-1}{2{\Gamma }^{\ast }W}=\frac{{\Gamma }^{\ast }-{\Gamma }_c}{2{\Gamma }^{\ast }},$$

$$|{\lambda }^{\pm }|=\sqrt{1-\frac{\tau +2}{\tau (\tau -1)}}.$$

$$\omega =\arctan \frac{\sqrt{\tau +2/\tau -4}}{\tau -2}.$$

$$\rho [t+1]$$

$$\Gamma [t+1]$$

$${\rho }^{\ast }$$

$${\Gamma }^{\ast }$$

$${\rho }^{\ast }=\frac{{\Gamma }^{\ast }-1}{2{\Gamma }^{\ast }},$$

$${\lambda }^{+}{\lambda }^{-}=\left(1-\frac{1}{\tau }\right)\left(1-\frac{2(A-1)}{A+u\tau +1}\right)+\frac{u(A-1{)}^2}{(A+u\tau +1)(2A+u\tau )}.$$

$$\left|{\lambda }^{\pm }\right|=\sqrt{{\lambda }^{+}{\lambda }^{-}}=1-\frac{2(A-1)+u}{2u\tau }+O({\tau }^{-2}),$$

$$P(X_i[t+1]=1)=(1-X_i[t])\left[1-\prod _j^{K_i}(1-P_{ij}{X}_j)\right],$$

$$P_{ij}[t+1]=P_{ij}[t]+\frac{1}{\tau }\left(\frac{A}{K}-P_{ij}[t]\right)-u{P}_{ij}[t]X_j[t].$$

$$\rho [t+1]$$

$$\sigma [t+1]$$

$${\rho }^{\ast }$$

$${\sigma }^{\ast }$$

$$|{\lambda }^{\pm }|=1-\left(\frac{(A-1)(2K-1)}{2uK}+\frac{1}{2}\right)\frac{1}{\tau }+O({\tau }^{-2}).$$

$$|{\lambda }^{\pm }|=1-\frac{A-1+u/2}{u\tau }+O({\tau }^{-2}),$$

$$\rho [t+1]$$

$$\Gamma [t+1]$$

$${\rho }^{\ast }=\frac{1}{\tau },{\Gamma }^{\ast }=\frac{{\Gamma }_c}{1-2/\tau }.$$

$$a={\frac{\partial F}{\partial \rho }]}_{\ast }$$

Full text

1 Introduction

Conservative self-organized critical (SOC) systems are by now well understood in the framework of out of equilibrium absorbing phase transitions [22, 30, 21, 50]. But, since natural systems that present power law avalanches (earthquakes, forest fires etc.) are dissipative in the bulk, compensation (drive) mechanisms have been proposed to make the models at least conservative on average. However, it seems that these mechanisms do not work so well: they produce “dirty criticality”, in which the system hovers around the critical region with stochastic oscillations (SO) [5, 4].

The literature reserves the acronym aSOC (adaptive SOC) to models that have an explicit dynamics in topology or network parameters (link deletion and creation, adaptive synapses, adaptive gains) [6, 39, 40, 23, 7, 16]. In the area of neuronal avalanches [3, 12, 2, 13], a well known aSOC system that presents SO is the Levina, Herrmann and Geisel (LHG) model [36, 37, 4] which uses continuous time integrate-and-fire neurons and dynamic synapses with short-term depression inspired by the work of Tsodyks and Markram [53, 54]. After that, similar models have been studied, for example excitable cellular automata [31] with LHG synapses [14, 9] and discrete time stochastic neurons with dynamic neuronal gains [7, 16, 15].

Not all SOqC systems (say, forest-fire models) are aSOC (which always have adaptive network parameters) but it seems that all aSOC systems are examples of SOqC. The exact origin of SO in aSOC systems is a bit unclear [5, 4, 16, 15]. Here, we examine some representative discrete time aSOC models at the mean-field (MF) level. We find that they evolve as 2d MF maps whose fixed point is a stable focus very close to a Neimark-Sacker-like bifurcation, which defines the critical point.

This kind of stochastic oscillation is known in the literature, sometimes called quasicycles [43, 38, 51, 56, 1, 11, 45] but, to produce them, one ordinarily needs to fine tune the system close to the bifurcation point. In contrast, for aSOC systems, there is a self-organization dynamics that tunes the system very close to the critical point [36, 4, 14, 9, 7, 16].

Although aSOC models show no exact criticality, they are very interesting because they can explain the coexistence of power law distributed avalanches, very large events (“dragon kings” [17, 44, 59]) and stochastic oscillations [16, 15]. Also, the adaptive mechanisms are biologically plausible and local, that is, they do not use non-local information to tune the system toward the critical region as occurs in other models [19, 18].

2 Network model with stochastic neurons

Our basic elements are discrete time stochastic integrate-and-fire neurons [25, 27, 26, 24, 34]. They enable simple and transparent analytic results [7, 16] but have not been intensively studied. We consider a fully connected topology with $i=1,\ldots,N$ neurons. Let $X_{i}$ be a firing indicator: $X_{i}[t]=1$ means that neuron $i$ spiked at time $t$ and $X_{i}[t]=0$ indicates that neuron $i$ was silent at time $t$. Each neuron spikes with a firing probability function $\Phi(V_{i}[t])$ that depends on a real valued variable $V_{i}[t]$ (the membrane potential of neuron $i$ at time $t$). Notice that, although the firing indicator is binary, the model is not a binary cellular automaton but corresponds to a stochastic version of leaky integrate-and-fire neurons.

The firing function can be any general monotonically increasing function $0\leq\Phi(V)\leq 1$. For mathematical convenience, we use the so-called rational function [16]:

[TABLE]

where $\Gamma$ is the neuronal gain (the derivative $d\Phi/dV$ for small $V$) and $\Theta$ is the step Heaviside function. This firing function is shown in Fig. 1a.

In a general case the membrane voltage evolves as:

[TABLE]

where $\mu$ is a leakage parameter and $I_{i}[t]$ is an external input. The synaptic weights $W_{ij}$ ($W_{ii}=0$) are real valued with average $W$ and finite variance. In the present case we study only excitatory neurons ($W_{ij}>0$) but the model can be fully generalized to an excitatory-inhibitory network with a fraction $p$ of excitatory and $q=1-p$ inhibitory neurons [32]. The voltage is reset to zero after a firing event in the previous time step. Since $\Phi(0)=0$, this means that two consecutive firings are not allowed (the neuron has a refractory period of one time step).

In this paper we are interested in the second order absorbing phase transition that occurs when the external fields $I_{i}$ are zero. Also, the universality class of the phase transition is the same for any value of $\mu$ [16], so we focus our attention to the simplest case $\mu=0$. In the MF approximation, we substitute $X_{i}$ by its mean value $\rho=\left\langle X_{i}\right\rangle$, which is our order parameter (density of firing neurons or density of active sites). Then, equation (2) reads $V[t+1]=W\rho[t]$, where $W=\left\langle W_{ij}\right\rangle$. From the definition of the firing function Eq. (1), we have:

[TABLE]

where $p(V)[t]$ is the voltage density at time $t$.

To proceed with the stability analysis of fixed points, it suffices to obtain the map for $\rho[t]$ close to stationarity. After a transient where all neurons spike at least one time, equations (2) lead to a voltage density that has two Dirac peaks, $p(V)[t+1]=\rho[t]\delta(V)+(1-\rho[t])\delta(V-W\rho[t])$. Inserting $p(V)$ in equation (3), we finally get:

[TABLE]

The $1-\rho[t]$ fraction corresponds to the density of silent neurons in the previous time step. We call this the static model because $W$ and $\Gamma$ are fixed control parameters.

The order parameter $\rho$ has two fixed points: the absorbing state $\rho^{0}=0$, which is stable (unstable) for $\Gamma W<1$, ($>1$) and a non-trivial firing state:

[TABLE]

which is stable for $\Gamma W>1$. This means that $\Gamma W=1$ is a critical line of a continuous absorbing state transition (transcritical bifurcation), see Fig. 1b and 1c. If parameters are put at this critical line one observes well behaved power laws for avalanche sizes and durations with the mean-field exponents $-3/2$ and $2$ respectively [7]. However, such fine tuning should not be allowed for systems that are intended to be self-organized in criticality.

2.1 Stochastic neurons model with dynamic synapses

Adaptive SOC models try to turn the critical point into an attractive fixed point of some homeostatic dynamics. For example, in the LHG model, the spike of the presynaptic neuron produces depression of the synapse, which recovers within some large time scale [36, 4].

A model with stochastic neurons and LHG synapses uses the same equations (1 and 2),but now the synapses change with time [7]:

[TABLE]

Here, $\tau$ is the recovery time toward a baseline level $A$ and $0<u<1$ is the fraction of synaptic strength that is lost when the presynaptic neuron fires ($X_{j}=1$). From now, we always use $\Delta t=1$ ms, the typical width of a spike.

2.2 Stochastic neurons model with dynamic neuronal gains

Instead of adapting synapses toward the critical line $W_{c}=1/\Gamma$, we can adapt the gains toward the critical condition $\Gamma_{c}=1/W$, see Fig. 1c. This can be modeled as individual dynamic neuronal gains $\Gamma_{i}[t]$ ($i=1,\ldots,N$) that decrease by a factor $u$ if the neuron fires (diminishing the probability of subsequent firings) with a recovery time $1/\tau$ toward a baseline level $A$ [7]:

[TABLE]

which is very similar to the LHG dynamics. Notice, however, that here we have only $N$ equations for the neuronal gains instead of $N(N-1)$ equations for dynamic synapses, which allows the simulation of much larger systems. Also, the neuronal gain depression occurs due to the firing of the neuron $i$ (that is, $X_{i}[t]=1$) instead of the firing of the presynaptic neuron $j$. The biological location is also different: adaption of neuronal gains (that produces firing rate adaptation) is a process that occurs at the axonal initial segment (AIS) [33, 49] instead of dendritic synapses.

Like the LHG model, this dynamics inconveniently has three parameters ($\tau,A$ and $u$). Recently, we proposed a simpler dynamics with only one parameter [16, 15]:

[TABLE]

Averaging over the sites, we obtain the 2d MF map:

[TABLE]

2.3 Stability analysis for stochastic neurons with simplified neuronal gains

This case with a single-parameter dynamics ($\tau$) for the neuronal gains has the simplest analytic results, so it will be presented first and with more detail. For finite $\tau$, $\rho^{0}=0$ is no longer a solution, see equation (10), and the 2d map has a single fixed point $(\rho^{*},\Gamma^{*})$:

[TABLE]

where $\Gamma_{c}=1/W$. The relation between $\rho^{*}$ and $\Gamma^{*}$ is:

[TABLE]

which resembles the expression for the transcritical phase transition in the static system, see equation (5). Here, however, $\Gamma^{*}$ is no longer a parameter to be tuned but rather a fixed point of the 2d map to which the system dynamically converges. Notice that the critical point of the static model can be approximated for large $\tau$, with $(\rho^{*},\Gamma^{*})\rightarrow_{\tau\rightarrow\infty}(0,\Gamma_{c})$.

Performing a linear stability analysis of the fixed point (see Supplementary Information), we find that it corresponds to stable focus. The modulus of the complex eigenvalues is:

[TABLE]

For large $\tau$ we have $|\lambda^{\pm}|=1-O(1/\tau)$, with a Neimark-Sacker-like critical point occurring when $|\lambda^{\pm}|=1$ where the focus turns out indifferent. For example, we have $|\lambda^{\pm}|\approx 0.990$ for $\tau=100$, $|\lambda^{\pm}|\approx 0.998$ for $\tau=500$ and $|\lambda^{\pm}|\approx 0.999$ for $\tau=1,000$. Since, due to biological motivations, $\tau$ is in the interval of $100-1000$ ms, we see that the focus is at the border of losing their stability. We call this point Neimark-Sacker-like because, in contrast to usual Neimark-Sacker one, the other side of the bifurcation, with $|\lambda^{\pm}|>1$, does not exist.

2.4 Finite size fluctuations produce stochastic oscillations

The stability of the fixed point focus is a result for the MF map that represents an infinite system without fluctuations. However, fluctuations are present in any finite size system and these fluctuations perturb the almost indifferent focus, exciting and sustaining stochastic oscillations that hover around the fixed point.

Without loss of generality, we fix $W=1$ in the simulations for the simple model with one parameter dynamics defined by Eq. (8). In Fig. 2a, we show the SO for the firing density $\rho[t]$ and average gain $\Gamma[t]$ ($\tau=500$ and $N=100,000$). As observed in the original LHG model [4], the stochastic oscillations have a sawtooth profile with no constant amplitude or period. In Fig. 2b, we show the SO in the phase plane $\rho$ vs $\Gamma$ for $\tau=100,500$ and $1,000$.

For small amplitude (harmonic) oscillations, the frequency is given by (see Supplementary Information):

[TABLE]

The oscillation period, also for small amplitudes, is given by $T=2\pi/\omega$. The full oscillations are non-linear and their frequency and period depend on the amplitude.

Notice that, in the critical case, we have full critical slowing down: $\lim_{\tau\rightarrow\infty}\omega=\tau^{-1/2}\rightarrow 0,\lim_{\tau\rightarrow\infty}T\rightarrow\infty$. This means that, exactly at the critical point, the 2d map corresponds to a center without a period, and the SO would correspond to critical fluctuations similar to random walks in the $\rho$ vs $\Gamma$ plane. However, since for any physical/biological system the recovery time $\tau$ is finite, the critical point is not observable.

2.5 Stochastic oscillations and avalanches

For large $\tau$, part of the SO orbit occurs very close to the zero state $\rho^{0}=0$, see Figs. 2b-c. Due to finite-size fluctuations, the system frequently falls into this absorbing state. The orbits in phase space are interrupted. As in usual SOC simulations, we can define the size of an avalanche as the number $S$ of firing events between such zero states. After a zero state, we force a neuron to fire, to continue the dynamics, so that Figs. 2b-c are better understood as a series of patches (avalanches) terminating at $\rho[t]=0$, not as a single orbit.

The presence of the SO affects the distribution of avalanche sizes $P(S)$, see Fig. 3. For small $\tau$, $\rho^{*}$ is larger and it is more difficult for the SO to fall into the absorbing state. We observe a bump of very large avalanches (dragon king events). Increasing $\tau$, we move closer to the Neimark-Sacker-like critical point and observe power law avalanches with exponent $-3/2$ similar to those produced in the static model that suffers a transcritical bifurcation. So, by using a different mechanism (Neimar-Sacker-like versus transcritical bifurcation), our aSOC model can reproduce experimental data about power law neuronal avalanches [3, 12, 2, 13]. It also predicts that these avalanches can coexist with dragon kings events [44, 59]. We notice that Fig. 3 for $\tau=5000$ seems to be subcritical, but this is a finite-size effect [16]. The $P(S)$ distribution is not our main concern here and a more rigorous finite size analysis can be found in previous papers [7, 16].

2.6 Stochastic neurons with LHG dynamic gains

We now return to the case of dynamic gains with LHG dynamics, see equation (7). Without loss of generality, we use $W=1$, so that the static model has $\Gamma_{c}=1/W=1$. The map is given by:

[TABLE]

The trivial absorbing fixed point is $(\rho^{0},\Gamma^{0})=(0,A)$, stable for $A<1$, see Supplementary Material. For $A>1$, we have a non-trivial fixed point given by:

[TABLE]

There is a relation between $\rho^{*}$ and $\Gamma^{*}$ that resembles the relation between the order parameter and the control parameter in the static model:

[TABLE]

valid for $\Gamma^{*}>\Gamma_{c}=1$. Here, however, $\Gamma^{*}$ is a self-organized variable, not a control parameter to be finely tuned.

We notice that to set by hand $A=1$, pulling all gains toward $\Gamma_{i}=1$, produces the critical point $(\rho^{*},\Gamma^{*})=(0,1)$. Nevertheless, this is a fine tune that should not be allowed for SOC systems. We must use $A>1$, and reach the critical region only at the large $\tau$ limit: $\rho^{*}\approx(A-1)/\tau u$ and $\Gamma^{*}\approx\Gamma_{c}+2(A-1)/\tau u$.

Proceeding with the linear stability analysis, we obtain (see Supplementary Information):

[TABLE]

To first order in $\tau$, we have:

[TABLE]

which means that, for large $\tau$, the map is very close to the Neimark-Sacker-like critical point. As the stable focus approaches the critical value, the system exhibits oscillations as it approaches the fixed point due to demographic noise, leading to SO in finite-sized systems.

For example, with the typical values $A=1.05$ and $u=0.1$ [14, 9], we have $\left|\lambda^{\pm}\right|\approx 1-1/\tau$, which gives $\left|\lambda^{\pm}\right|\approx 0.990$ for $\tau=100$ and $\left|\lambda^{\pm}\right|\approx 0.999$ for $\tau=1,000$. In Fig. 4, we present the exact $\left|\lambda^{\pm}\right|$, the square root of equation (20), as a function of $\tau$ for several values of $A$ and $u$. For large $\tau,$ the frequency for harmonic oscillations is given by $\omega\simeq\sqrt{(A-1)/\tau}$ (see Supplementary Information).

The analysis of the stochastic neuron model with LHG synapses $W_{ij}[t]$, see equation (6), is very similar to the one above with LHG dynamic gains $\Gamma_{i}[t]$. The only difference is that we need to exchange $\Gamma$ by $W$. The eigenvalue modulus is the same, so that LHG dynamic synapses also produce SO. But, instead of presenting simulations in our complete graph system, which would involve $N(N-1)$ dynamic equations for the synapses, we prefer to discuss the LHG synapses for another system well known in the literature. We will also examine how SO depends on the system size $N$.

3 Excitable cellular automata with LHG synapses

We now consider an aSOC version of a probabilistic cellular automata well studied in the literature [31, 35, 46, 42, 57, 60]. It is a model for excitable media that yields a natural interpretation in terms of neuronal networks. Each site has $n$ states, $X=0$ (silent), $X=1$ (firing), $X=2,\ldots,n-1$ (refractory). Here we will use $n=2$.

Each neuron $i=1,\ldots,N$ has $K_{i}$ random neighbours (the average number in the network is $K=\left\langle K_{i}\right\rangle$) coupled by probabilistic synapses $P_{ij}\in[0,1]$ (here we use $K=10$, for implementation details see [31, 14, 9]). If the presynaptic neuron $j$ fires at time $t$, with probability $P_{ij}$ the postsynaptic neuron $i$ fires at $t+1$. The update is done in parallel, so the synapses are multiplicative, not additive like in the stochastic neuron model. All neurons that fire go to silence in the next step. Only neurons in the silent state can be induced to fire. Since $1-P_{ij}X_{j}$ is the probability that neighbour $j$ does not induces the firing of neuron $i$, we can write the update rule as:

[TABLE]

The control parameter of the static model is the branching ratio $\sigma=K\left\langle P_{ij}\right\rangle$ and the critical value is $\sigma_{c}=1$.

In the cellular automata (CA) model with LHG synapses, we have [14, 9]:

[TABLE]

3.1 Mean field stability analysis

The 2d MF map close to the stationary state is:

[TABLE]

where in the second line we multiplied equation (23) by $K$ and averaged over the synapses.

There is a trivial absorbing fixed point $(\rho^{0},\sigma^{0})=(0,A)$, stable up to $A=1$, see Supplementary Material. For $A>1$, there exists a single stable fixed point given implicitly by:

[TABLE]

which we can find numerically.

For this model we have numerical but no analytic results. However, we can find an approximate solution close to criticality, where $\rho^{*}$ is small (see Supplementary Information). Expanding in powers of $1/\tau$, we get:

[TABLE]

This confirms that the same scenario of a weakly stable focus appears in the CA model.

This last result is particularly interesting. If we put $K=N-1$, the CA network becomes a complete graph. Performing the limit $N\rightarrow\infty$, equation (28) gives:

[TABLE]

which is identical to the $\left|\lambda^{\pm}\right|$ value for the LHG stochastic neuron model, see equation (21).

Concerning the simulations, we must make a technical observation. In contrast to the static model [31], the relevant indicator of criticality here is no longer the branching ratio $\sigma$, but the principal eigenvalue $\Lambda\neq\sigma$ of the synaptic matrix $P_{ij}$, with $\Lambda_{c}=1$ [35, 57]. This occurs because the synaptic dynamics creates correlations in the random neighbour network [9]. So, the MF analysis, where correlations are disregarded, does not furnish the exact $\sigma$ or $\Lambda$ of the CA model. But there exists an annealed version of the model in which $\sigma^{*}=\Lambda^{*}$ and the MF analysis fully holds [14, 9].

In this annealed version, when some neuron fires, the depressing term $-uP_{ij}$ is applied to $K$ synapses randomly chosen in the network. This is not biologically realistic, but destroys correlations and restores the MF character of the model. Since all our analyses are done at the MF level, we prefer to present the simulation results using the annealed model. Concerning the SO phenomenology, there is no qualitative difference between the annealed and the original model defined by Eq. (23).

The full stochastic oscillations in a system with $N=128,000$ neurons and $\tau=500$ can be seen in Fig. 5a and b. The angular frequency for small amplitude oscillations is given by equation (92) (see Supplementary Information). We have $\omega\simeq\sqrt{(A-1)/\tau}$ for large $\tau$, which is the same behavior found for stochastic neurons. The power spectrum of the time evolution of $\rho$ and $\sigma$ is shown in Fig. 6. The peak frequency gets closer to the theoretical $\omega$ (vertical line) for larger system sizes because the oscillations have smaller amplitudes, going to the small oscillations limit.

3.2 Dependence of the stochastic oscillations on system size

Now we ask if the SO survive in the thermodynamic limit. In principle, given that fluctuations vanish in this limit, and we have always damped focus for finite $\tau$, the SO should also disappear when $1/N\rightarrow 0$. In contrast, Bonachela et al. claimed that, in the LHG model, the amplitude of the oscillations basically does not change with $N$ and is non-zero in the thermodynamic limit [4]. Indeed, they proposed that this feature is a core ingredient of SOqC.

In Fig. 7, we measured the average $\left\langle\sigma[t]\right\rangle$ and the standard deviation $\Delta\sigma$ of the $\sigma[t]$ time series of the annealed model. We used $\tau=320,500,1,000$ and $2,000$ and system sizes from $N=4,000$ to $1,024,000$. We interpret our findings as a $\tau$ dependent crossover phenomenon due to a trade-off between the level of fluctuations (which depends on $N$) and the level of dampening (which depends on $\tau$): for a given $\tau$, a small $N$ can produce sufficient fluctuations so that the SO are sustained without change of $\Delta\sigma$. Nonetheless, starting from some $N(\tau)$, the fluctuations are not sufficient to compensate the dampening given by $|\lambda^{\pm}|<1$ and $\Delta\sigma$ starts to decrease for increasing $N$. The larger the $\tau$, the less damped is the focus and the SO survive without change to a larger $N$ (the plateau in Fig. 7b is more extended to the left).

So, the conclusion is that, for any finite $\tau$, the SO do not survive up to the thermodynamic limit, although they can be observed in very large systems. This is compatible with Fig. 6 where the power spectrum decreases with $N$. We also see in Fig. 7a that $\left\langle\sigma[t]\right\rangle\rightarrow\sigma^{*}$ for increasing $N$, so the system settles without variance at the MF fixed point $(\rho^{*},\sigma^{*})$ in the infinite size limit.

We can reconcile our findings with those of the LHG model [36, 4] remembering an important technical detail: these authors used a synaptic dynamics with $\tau=\tau_{0}N$, in an attempt to have the fixed point converging to the critical one in the thermodynamic limit (the same scaling is used in other models [23, 14]). In systems with this scaling, the fluctuations decay with $N$ but, at the same time, the dampening controlled by $\lambda(\tau(N))$ also decays with $N$. For this scaling, we can accord that the SO survive in the thermodynamic limit. However, we already emphasized that this is a non-biological and non-local scaling choice [9, 7, 16] because the recovery time $\tau$ must be finite and the knowledge of the network size $N$ is non-local.

However, from a biological point of view, this discussion is not so relevant: although the finite-size fluctuations (called demographic noise in the literature [38, 51, 11, 45]) disappear for large $N$, external (environmental) noise of biological origin, not included in the models, never vanishes in the thermodynamic limit. So, for practical purposes, the SO and the associated dragon kings would always be present in more realistic noisy networks and experiments.

4 Discussion and Conclusion

We now must stress what is new in our findings. Standard SOC models are related to static systems presenting an absorbing state phase transition [22, 30, 21, 50]. At the MF level, these static systems are described by an 1d map $\rho[t+1]=F(\rho[t])$ for the density of active sites, and the phase transition corresponds to a transcritical bifurcation where the critical point is an indifferent node. This indifferent equilibrium enables the occurrence of scale-invariant fluctuations in the $\rho[t]$ variable, that is, scale-invariant avalanches but not stochastic oscillations.

In aSOC networks, the original control parameter of the static system turns out an activity dependent variable leading to 2d MF maps. We performed a MF fixed point stability analysis for three systems: 1) a stochastic neuronal model with one-parameter neuronal gain dynamics, 2) the same model with LHG gain dynamics, and 3) the CA model with LHG synapses. For the first one we obtained very simple and transparent analytic results; for the LHG dynamics we also got analytic, although more complex, results; finally, for the CA model, we were able to obtain a first order approximation for large recovery time $\tau$. Curiously, in the limit of $K=(N-1)\rightarrow\infty$ neighbours, the stochastic neuron model and the CA model have exactly the same first order leading term. The complex eigenvalues have modulus $|\lambda^{\pm}|\approx 1-O(1/\tau)$ and, for large $\tau$, the fixed point is a focus at the border of indifference. This means that, in finite size systems, the dampening is very low and fluctuations can excite and sustain stochastic oscillations.

About the generality of our results, we conjecture they are generic and valid for the whole class of aSOC models, from networks with discrete deletion-recovery of links [39, 40, 23] to continuous depressing-recovering synapses [36, 4, 14, 9] and neurons with firing rate adaptation [7, 16, 15]). At the mean-field level, all of them are described by similar two-dimensional dynamical systems: one variable for the order parameter, another for the adaptive mechanism. For example, the prototypical LHG model [36, 4] uses continuous time LIF neurons (which is equivalent to setting $\Gamma\rightarrow\infty$ and $\theta,\mu>0$ in our model). Although with more involved calculations, in principle one could do a similar mean-field calculation and obtain a 2d dynamical system with an almost unstable focus close to criticality, explaining the SO observed in that model. In another aspect of generality, stochastic oscillations have been observed previously in other SOqC systems (forest-fire models [29, 55] and some special sandpile models [20, 52]) that share a similar MF description with aSOC systems. Indeed, it seems that stochastic oscillations is a distinctive feature of self-organized quasi-criticality, as defined by Bonachela et al. [5, 4].

We also emphasize that the fact that our network has only excitatory neurons is not a limitation of this study. In a future work [32], we will show that our model can be fully generalized to an excitatory-inhibitory network very similar to the Brunel model [8], with the same results.

So, it must be clear that what is new here is not the stochastic oscillations (quasicycles) in biosystems, since there is a whole literature about that [43, 38, 51, 56, 1, 11, 45]. What is new here is the interaction of the SO with a critical point with an absorbing state. This interaction interrupts the oscillations, producing the phenomenology of avalanches and dragon kings. This is our novel proposal for a mechanism that produces dragon kings coexisting with limited power law avalanches.

In conclusion, contrasting to standard SOC, aSOC systems present stochastic oscillations that will not vanish in the thermodynamic limit if external (environmental) noise is present. But what seems to be a shortcoming for neural aSOC models could turn out to be an advantage. Since the adaptive dynamics with large $\tau$ has good biological motivation, it is possible that SO are experimentally observable, providing new physics beyond the standard model for SOC. The presence of the Neimark-Sacker-like bifurcation affects the distribution of avalanche sizes $P(S)$, creates dragon king events, and all this phenomenology can be measured [17, 48, 47]. In particular, our raster plots results are compatible with system-sized events (synchronized fires) recently observed in neuronal cultures [44, 59]. So, we have experimental predictions that differ from standard SOC based on a transcritical bifurcation (and also from criticality models with Griffiths phases [41, 28]). We propose that the experimental detection of stochastic oscillations and dragon-kings should be the next experimental challenge in the field of neuronal avalanches.

Finally, we speculate that non-mean field models (square or cubic lattices) with dynamic gains and SO could be applied to the modeling of dragon-kings (large quasiperiodic ”characteristic”) earthquakes that coexist with power-law Gutemberg-Richter distributions for the small events [10, 58]. This will be pursued in another work.

Methods

Numerical Calculations: Numerical calculations were done by using MATLAB softwares. Simulation procedures: Simulation codes were made in Fortran90 and C++11.

Acknowledgements

This article was produced as part of the activities of FAPESP Research, Innovation and Dissemination Center for Neuromathematics (Grant No. 2013/07699-0, S. Paulo Research Foundation). We acknowledge financial support from CAPES, CNPq, FACEPE, and Center for Natural and Artificial Information Processing Systems (CNAIPS)-USP. LB thanks FAPESP (Grant No. 2016/24676-1). AAC thanks FAPESP (Grants No. 2016/00430-3 and No. 2016/20945-8).

Author contributions

AAC and JGFC performed the simulations and prepared all the figures. JGFC, LB, MC and OK made the analytic calculations. All authors contributed to the writing of the manuscript.

Competing interests

The authors declare that they have no competing financial interests.

Supplementary Information

Stability analysis of the stochastic

neuron model with one parameter ($\tau$) gain dynamics

The MF map is:

[TABLE]

Notice that $(\rho^{0}=0,\Gamma)$ is not a fixed point. The single fixed point is:

[TABLE]

To perform the fixed point stability analysis we calculate partial derivatives evaluated at the fixed point (denoted by $\left.\right]_{*}$):

[TABLE]

The eigenvalues of the system Jacobian matrix at the fixed point $(\rho^{*},\Gamma^{*})$ are determined by

[TABLE]

So, the characteristic equation is:

[TABLE]

That can be written as:

[TABLE]

with solutions:

[TABLE]

This equation has no meaning for $\tau<2$ since this would give $\rho^{*}=1/\tau>1/2$. When $2<\tau<2+\sqrt{2}$ we have $\Delta>0$, $\lambda$ real and $|\lambda|<1$, that is, the fixed point $(\rho^{*},\Gamma^{*})$ is a stable node. However, when $\tau>2+\sqrt{2}=3.4142\ldots$, we have $\Delta<1$ and two complex roots $\lambda^{\pm}$ with $|\lambda|=\sqrt{\lambda^{+}\lambda^{-}}<1$, that is, we have a stable focus with:

[TABLE]

which leads to equation (13).

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

[TABLE]

When $\tau\rightarrow\infty$,

[TABLE]

Stability analysis of the stochastic neuron model

with LHG gain dynamics

The MF map is given by:

[TABLE]

We have two fixed points. The absorbing state fixed point:

[TABLE]

and the non trivial fixed point $(\rho^{*},\Gamma^{*})$:

[TABLE]

Partial derivatives of $F$ and $G$ with respect to $\rho$ and $\Gamma$ are:

[TABLE]

Evaluating at the trivial fixed point, we get:

[TABLE]

The eigenvalues are found from $(a-\lambda)(d-\lambda)-bc=0$, or:

[TABLE]

with solution:

[TABLE]

Notice that $\Delta\geq 0$ so the eigenvalues are always real: the fixed point is an attractive or repulsive node. The final solution of Eq. 60 is:

[TABLE]

This means that the fixed point $(\rho^{0}=0,\Gamma^{0}=A)$ is stable for $\lambda^{+}=A<1$. The critical point is $A=1$.

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 $(\rho^{*},\Gamma^{*})$ is the unique stable solution. For $A=1$, the point $(\rho^{0}=0,\Gamma^{0}=1)$ is an indifferent node.

Evaluating the derivatives at the non trivial fixed point, we get:

[TABLE]

As before, the eigenvalues of the Jacobian matrix at the fixed point $(\rho^{*},\Gamma^{*})$ are determined by equation (37). After some algebra, we get:

[TABLE]

For large $\tau$, we have:

[TABLE]

which means that, for large $\tau$, 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 $\omega\equiv|\omega^{\pm}|=\arctan{\sqrt{\frac{4D}{T^{2}}-1}}$, where $D$ and $T$ are respectively the determinant and trace of the Jacobian matrix. When $\tau\rightarrow\infty$,

[TABLE]

which is similar to the $\tau^{-1/2}$ behavior of the simple gain dynamics, Eq. (45).

Stability analysis of the

Cellular Automata model with LHG synapses

The 2d MF map near the stationary state is:

[TABLE]

where in the second line we multiplied equation (23) by $K$ and averaged over the synapses. As before, the stability of the map is given by the eigenvalues of the Jacobian matrix $J(\rho,\sigma)$ at the fixed points. The matrix elements are:

[TABLE]

Thus, we can plug in the solutions of equations (26) and (26) in equations (73-76) and find the eigenvalues of $J(\rho,\sigma)$.

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

[TABLE]

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

[TABLE]

which is the same result of Eqs. (54–57). This matrix has two eigenvalues: $1-1/\tau$ and $A$. The eigenvalue $1-\tau$ corresponds to the $\sigma$ 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 $\rho^{*}$ is small, we can find an approximate solution for $J(\rho,\sigma)$ at the nontrivial fixed point. Expanding the binomial in equation (24) to first order and solving for $\rho^{*}$, we find:

[TABLE]

where we used equation (26) to compute $\sigma^{*}$. Note that, as before, this solution exists only if $A-1>0$. This is the exact condition to the absorbing state solution be unstable.

Inserting in equations (73-76), we have:

[TABLE]

where

[TABLE]

Note that they are all positive quantities. The eigenvalues of $J$ are:

[TABLE]

For systems to be close to criticality, we must have $\tau\gg 1$. Thus, $a$, $b$ and $d$ are of $O(1/\tau)$, while $c=O(1)$. Therefore, the $4bc$ term dominates the square root argument in equation (89). Hence, we have two complex conjugate eigenvalues with modulus and argument given by:

[TABLE]

We can expand $|\lambda^{\pm}|$ in powers of $1/\tau$:

[TABLE]

This confirms that the same scenario of weakly stable focus appears in the CA model. The harmonic frequency is:

[TABLE]

To first order in $\tau$, we get $\omega\simeq\sqrt{(A-1)/\tau}$, which shows that the frequency of oscillations is vanishingly small for large $\tau$, being identical to equation (70).