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

TL;DR

This paper investigates networks of stochastic spiking neurons, demonstrating self-organized supercriticality with power-law avalanches and oscillations, relevant to understanding brain rhythms and phenomena like epilepsy.

Contribution

It introduces the concept of self-organized supercriticality in neuronal networks with dynamic gains, showing robustness and coexistence of avalanches and oscillations.

Findings

Networks reach a supercritical state with power-law avalanches.

Neuronal gains induce collective oscillations matching brain rhythms.

Self-organized supercriticality is robust to parameter changes.

Abstract

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.

Equations

$$V_{i}[t+1]=\left\{\begin{array}[]{lcl}\displaystyle 0&&\hbox{if $X_{i}[t]=1$,}\\ \displaystyle\mu V_{i}[t]+I_{i}[t]+\frac{1}{N}\sum_{j=1}^{N}W_{ij}X_{j}[t]&&\hbox{if $X_{i}[t]=0$.}\end{array}\right.$$

$$\Phi (V)=\frac{\Gamma (V-V_T)}{1+\Gamma (V-V_T)}\Theta (V-V_T).$$

$$\rho [t]=\frac{1}{N}\sum _{j=1}^N{X}_j[t].$$

$$\rho [t]={\int }_{-\infty }^{\infty }\Phi (V)p[t](V)dV,$$

$$U_0$$

$$U_k$$

$${\eta }_k$$

$$\rho$$

$$\rho \approx \left(\frac{1}{2+\mu +{\mu }^2/(1-\mu )}\right)\frac{\Gamma -{\Gamma }_C}{\Gamma }\propto {\Delta }_{\Gamma }^{\beta },$$

$$\rho =\rho \Phi (0)+(1-\rho )\Phi (I+W\rho ),$$

$$2\Gamma W{\rho }^2-(\Gamma W+2\Gamma (V_T-I)-1)\rho +\Gamma (V_T-I)=0,$$

$${\rho }^{\pm }$$

$$\Delta$$

$$\rho (W)$$

$${\Gamma }_C=1/W_C,$$

$${\left({\Gamma }_C{W}_C+2{\Gamma }_C(V_T-I)-1\right)}^2=8{\Gamma }_C^2{W}_C(V_T-I),$$

$${\Gamma }_C{W}_C$$

$${\Gamma }_C$$

$${\rho }_C=\frac{\sqrt{V_T-I}}{\sqrt{2W_C}}=\frac{1}{2}\frac{\sqrt{2{\Gamma }_C(V_T-I)}}{1+\sqrt{2{\Gamma }_C(V_T-I)}}.$$

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

$${\Gamma }_i[t+1]={\Gamma }_i[t]+\frac{1}{\tau }(A-{\Gamma }_i[t])-u{\Gamma }_i[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 [t].$$

$$\Gamma [t+1]=\Gamma [t]+\frac{1}{\tau }\Gamma [t]-\rho [t]\Gamma [t],$$

$$\frac{1}{\tau }{\Gamma }^{\ast }={\rho }^{\ast }{\Gamma }^{\ast }.$$

$${\Gamma }^{\ast }=\frac{{\Gamma }_C}{1-2/\tau }.$$

$${\Gamma }^{\ast }={\Gamma }_C\left(1+\frac{2}{\tau }\right).$$

$${\Gamma }^{\ast }\approx 1.002{\Gamma }_C.$$

$${\rho }^{\ast }=\frac{1}{\tau }$$

$$\rho$$

$${\eta }_k$$

$$U_0$$

$$U_k$$

$$=$$

$$\rho$$

$$\Phi (U_k)=\frac{\Gamma U_k}{1+\Gamma U_k}\approx \Gamma U_k-{\Gamma }^2{U}_k^2.$$

Full text

\articlenumber

x \doinum10.3390/—— \pubvolumexx

\externaleditorAcademic Editor: name \historyReceived: date; Accepted: date; Published: date

\TitleSelf-Organized Supercriticality and Oscillations in Networks of Stochastic Spiking Neurons \AuthorAriadne A. Costa 1,2, Ludmila Brochini 3 and Osame Kinouchi 4,* \AuthorNamesAriadne de Andrade Costa, Ludmila Brochini and Osame Kinouchi

\corresCorrespondence: osame@ffclrp.usp.br

1 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 Herz and Hopfield (1995); Beggs and Plenz (2003); Kinouchi and Copelli (2006). 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 Beggs and Plenz (2003), these behaviors have been reported also in biological networks, see reviews  Chialvo (2010); Marković and Gros (2014); Hesse and Gross (2015).

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 Kinouchi and Copelli (2006); Beggs (2008); Shew et al. (2009); Massobrio et al. (2015). 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 de Arcangelis et al. (2006); Pellegrini et al. (2007); de Arcangelis and Herrmann (2010); De Arcangelis (2012); De Arcangelis and Herrmann (2012); van Kessenich et al. (2016), and self-organization of dynamic synapses Levina et al. (2007, 2009); Bonachela et al. (2010); Costa et al. (2015); Campos et al. (2017) following Tsodyks and Markran Tsodyks and Markram (1997); Tsodyks et al. (1998). 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) Kole and Stuart (2012); Ermentrout et al. (2001); Benda and Herz (2003); Buonocore et al. (2016), can also lead to self-organized criticality Brochini et al. (2016). 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 Tang et al. (2001) – and thus $10^{4}N$ dynamic equations – but only $N$ equations for the neuronal gains.

It has also been observed in Brochini et al. Brochini et al. (2016) that, for achieving exact SOC, all papers in the literature used a time scale $\tau$ for synaptic recovery proportional to $N$, with $N\rightarrow\infty$. The use of this non-local information ($N$), and a diverging recovery time $\tau$, is not plausible biologically. A similar time scale $\tau$ is present in the neuronal gain recovery. When we use a biological range for $\tau$ 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 $\tau$, 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 $\tau$ 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 $\tau$. Surprisingly, we found collective oscillations produced by the gain dynamics that coexist with neuronal avalanches.

2 The model

We consider a fully connected network composed of $i=1,\ldots,N$ discrete-time stochastic neurons Brochini et al. (2016); Gerstner and van Hemmen (1992); Gerstner and Kistler (2002); Galves and Löcherbach (2013); Larremore et al. (2014). The synapses transmit signals from some presynaptic neuron $j$ to a postsynaptic neuron $i$ with synaptic strength $W_{ij}$. The Boolean variable $X_{i}[t]\in\left\{0,1\right\}$ denotes whether neuron $i$ fired between $t$ and $t+1$ and $V_{i}[t]$ corresponds to its membrane potential at time $t$. Firing $X_{i}[t+1]=1$ occurs with probability $\Phi(V_{i}[t])$, which is called the firing function Galves and Löcherbach (2013); Larremore et al. (2014); De Masi et al. (2015); Duarte and Ost (2016); Duarte et al. (2015); Galves and Löcherbach (2016).

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 $\mu\in[0,1]$, which models the effect of a current leakage.

The neuron membrane potentials evolve as:

[TABLE]

This is a special case of the general model from Galves and Löcherbach (2013) with the filter function $g(t-t_{s})=\mu^{t-t_{s}}$, where $t_{s}$ is the time of the last firing of neuron $i$ Brochini et al. (2016). 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\leq\Phi(V)\leq 1$ is sigmoidal, that is, monotonically increasing. We also assume that $\Phi(V)$ is zero up to some threshold potential $V_{T}$. If $\Phi$ is the shifted Heaviside step function $\Phi(V)=\Theta(V-V_{T})$, we have a deterministic discrete-time leaky integrate-and-fire (LIF) neuron. Any other choice for $\Phi(V)$ gives a stochastic neuron.

In Brochini et al. Brochini et al. (2016) we have studied a linear saturating function with neuronal gain $\Gamma$ similar to that used in Larremore et al. (2014). Here we study the so called rational function that does not have a saturating potential, see Fig. 1(a):

[TABLE]

Notice that we recover the deterministic LIF model $\Phi(V)=\Theta(V-V_{T})$ when $\Gamma\rightarrow\infty$. 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 Brochini et al. (2016).

3 Mean-field calculations

The network’s activity is measured by the fraction $\rho[t]$ of firing neurons (or density of active sites):

[TABLE]

The density of active neurons $\rho[t]$ can be computed from the probability density $p[t](V)$ of potentials at time $t$:

[TABLE]

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 $\rho[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\neq 0$.

In the mean-field limit, we assume that the synaptic weights $W_{ij}$ follow a distribution with average $W=\left\langle W_{ij}\right\rangle$ and finite variance. By disregarding correlations, the term in Eq. (1) corresponding to the sum of all presynaptic inputs simplifies to $W\rho[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 $\mu$ and the parameters that define the function $\Phi(V)$, that is, $\Gamma$ and $V_{T}$. In Brochini et al. Brochini et al. (2016) it is shown that the stationary $p(V)$ is composed of delta peaks with height $\eta_{k}$ situated at voltages $U_{k}$ given by:

[TABLE]

for all $k\geq 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 $\sum_{k=0}^{\infty}\eta_{k}\;=\;1$ must be included explicitly. Equations (6–8) can be solved numerically for any firing function $\Phi$, so this result is very general.

4 Results

4.1 Phase transitions for the rational $\Phi(V)$

In terms of non-equilibrium statistical physics, $\rho$ is the order parameter, $I$ is an uniform external field and $\Gamma$ and $W$ are the main control parameters. The activity $\rho$ also depends on $V_{T}$ and $\mu$.

4.1.1 The case with $\mu>0,I=0,V_{T}=0$

By using Eqs. (5)-(8), we obtain numerically $\rho(W,\Gamma)$ for several values of $\mu>0$, for the case with $I=0,V_{T}=0$ (Fig. 1(b)). Only the first $100$ peaks $(U_{k},\eta_{k})$ were considered, since, for the given $\mu$ and $\Phi$, there was no significant probability density beyond that point. The same numerical method can be used for studying the cases $I\neq 0,V_{T}\neq 0$.

We also obtained an analytic approximation (see Appendix) for small $\rho$:

[TABLE]

where $\Gamma_{C}=(1-\mu)/W$ defines the critical line and $\Delta_{\Gamma}=(\Gamma-\Gamma_{C})/\Gamma$ is the reduced control parameter. So, the critical exponent for the order parameter near criticality is $\beta=1$, characteristic of the mean-field directed percolation (DP) universality class Hinrichsen (2000). We also compare Eq. (9) with the numerical results for $\rho(\Gamma,\mu)$ in Fig. 1(c).

4.1.2 Analytic results for $\mu=0$

In the case $\mu=0$ it is possible to do a simple mean-field analysis valid for $N\rightarrow\infty$. This case is illustrative because it presents all phase transitions that occur with $\mu>0$.

When $\mu=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\rho$. Eq. (8) simplifies to:

[TABLE]

since $\eta_{0}=\rho$ and $\eta_{1}=1-\rho$.

By inserting the function Eq. (2) in Equation (10), and remembering that $\Phi(0)=0$, we get:

[TABLE]

with solutions:

[TABLE]

4.1.3 The case with $I=0,V_{T}=0$: continuous transition

For $V_{T}=I=0$, we have:

[TABLE]

where the phase transition line is

[TABLE]

and the critical exponent is $\beta=1$. This corresponds to a standard mean-field continuous (second order) absorbing state phase transition (Fig. 2(a) and 2(b) with $V_{T}=I=0$).

4.1.4 The case with $I<V_{T}$: discontinuous transition

For $I<V_{T}$, we have discontinuous (first order) phase transitions when $\Delta=0$, see Eq. (13):

[TABLE]

which, after some algebra, leads to the phase transition lines:

[TABLE]

which have the correct limit, Eq. (16), when $V_{T}\rightarrow 0,I\rightarrow 0$. The transition discontinuity is:

[TABLE]

In Fig. 2(a) we show examples of the phase transitions, which occurs when the unstable point $\rho^{-}$ collapses with the stable point $\rho^{+}$. It is important to notice that, for any $V_{T}>0$, the unstable point never touches the absorbing point $\rho_{0}=0$, so the zero solution is always stable. Only for the case $V_{T}=0$ the solution $\rho_{0}$ looses stability and $\rho^{+}$ is the unique solution above the critical line $\Gamma_{C}=1/W_{C}$. Fig. 2(b) gives the phase diagram $\Gamma\times W$ for some values of $V_{T}$ with $I=0$, see Eq. (19). Finally, we give the phase diagram for the variables $\Gamma W$ versus $\Gamma(V-I)$, see Fig. 3 and Eq. (18).

4.2 Self-organized supercriticality (SOSC)

through dynamic gains with $\mu=0$, $I=0$, $V_{T}=0$

If we fine tune the model to some point in the critical line $\Gamma_{C}=1/W$, we can observe perfect neuronal avalanches with size distribution $P_{S}(s)\propto s^{-3/2}$ and duration distribution $P_{D}(d)\propto d^{-2}$ Brochini et al. (2016). As expected, these are mean-field exponents fully compatible with the experimental results Beggs and Plenz (2003); Hesse and Gross (2015).

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]$ Levina et al. (2007, 2009); Bonachela et al. (2010). For example, in discrete time Costa et al. (2015); Campos et al. (2017):

[TABLE]

where $\tau$ is a synaptic recovery time, $A$ is an asymptotic value and $u\in[0,1]$ is the fraction of depletion of neurotransmitter vesicles when the presynaptic neuron fires.

In Brochini et al. Brochini et al. (2016), we proposed a new self-organization mechanism based in dynamic neuronal gains $\Gamma_{i}[t]$ while keeping the synapses $W_{ij}$ fixed Brochini et al. (2016). 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 Ermentrout et al. (2001); Benda and Herz (2003).

The dynamics for the neuronal gains studied in Brochini et al. (2016) has a form similar to that used in Levina et al. (2007); Bonachela et al. (2010); Costa et al. (2015); Campos et al. (2017) for synapses:

[TABLE]

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 ($\tau\in[1,\infty],A\in[1/W,\infty],u\in[0,1]$) that must be fully explored to characterize the stationary value $\Gamma^{*}(\tau,A,u,N)$. Here we propose a new simplified dynamics with a single free parameter, the gain recovery time $\tau$:

[TABLE]

The self-organization mechanism can be viewed in Fig. 4. So, we reduce our parametric study to determine the curves $\Gamma^{*}(1/\tau,1/N)$, see Figs. 5(a) and 5(b). The fluctuations measured by the standard deviation $SD$ of the $\Gamma[t]$ time series, after the transient, diminishes for increasing $\tau$ (Fig. 5(c)) and probably goes to zero for $\tau\rightarrow\infty$, in accord to Campos et al. Campos et al. (2017). However, in contrast to this idealized $\tau\rightarrow\infty$ limit, as discussed in the references Bonachela et al. (2010); Campos et al. (2017), the fluctuations do not converge to zero for finite $\tau$ in the thermodynamic limit $N\rightarrow\infty$ (see Fig. 5(d)). This occurs because, for low $\tau$, the adaptation mechanism produces oscillations of $\Gamma[t]$ around the value $\Gamma^{*}(\tau)$

We can do a mean-field analysis of Eq. (23) to find the value $\Gamma^{*}(\tau)$. Denote the average gain as $\Gamma[t]=\left\langle\Gamma_{i}[t]\right\rangle$. Averaging over the sites, we have :

[TABLE]

since $\rho[t]=\left\langle X_{i}[t]\right\rangle$. In the stationary state, we have $\Gamma[t+1]=\Gamma[t]=\Gamma^{*},\>\rho[t]=\rho^{*}$, so:

[TABLE]

A solution is $\Gamma^{*}=0$, but this is unstable, see Eq. (23). Another solution is obtained by inserting Eq. (15), $\rho^{*}=(\Gamma^{*}-\Gamma_{C})/(2\Gamma^{*})$, in Eq. (25):

[TABLE]

Notice that this is valid only when Eq. (15) is valid, that is, for $\Gamma^{*}\geq 1/W$. Also, Eq. (26) presumes that $\rho^{*}$ is a stable fixed point, which can not be true for some interval of values of $\tau$, see below.

A first order approximation leads to:

[TABLE]

This mean-field calculation shows that, if $\tau\rightarrow\infty$, we obtain an exact SOC state $\Gamma^{*}\rightarrow\Gamma_{C}$. Or, for finite networks, it would be required a scaling $\tau=O(N^{a})$ with an exponent $a>0$, as done previously for dynamic synapses Levina et al. (2007); Bonachela et al. (2010); Costa et al. (2015); Campos et al. (2017). However, this scaling for $\tau$ cannot be justified biologically.

So, biology requires a finite recovery time $\tau$ 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 $\tau=1000$ ms (assuming 1 time step equals 1 ms in the model):

[TABLE]

Even a more conservative value $\tau=100$ ms gives $\Gamma^{*}\approx 1.02\>\Gamma_{C}$. Although not perfect SOC Marković and Gros (2014), 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:

[TABLE]

showing that the network presents supercritical activity for any finite $\tau$. This result, however, is valid only in the infinite size limit. For finite size networks, fluctuations interrupt this constant $\rho^{*}=1/\tau$ activity, leading the system to the absorbing state and defining the end of the avalanches.

Simulations reveal that these fixed points $(\Gamma^{*},\>\rho^{*})$ correspond only to mean values around which both $\Gamma[t]$ and $\rho[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 $\Gamma[t]$ oscillations).

5 Discussion

We examined a network of stochastic spiking neurons with a rational firing function $\Phi$ 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 $\Phi$ function, but not the usual one where the synaptic strength $W$ and the neuronal gain $\Gamma$ 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) Kole and Stuart (2012). Our dynamic $\Gamma_{i}[t]$ mimics the well known phenomenon of spike frequency adaptation Ermentrout et al. (2001); Benda and Herz (2003) and is a one-parameter simplification of the three-parameter dynamics studied by us in Brochini et al. (2016).

We observe that this gain dynamics is equivalent to approach the critical line with fixed $W$ and variable $\Gamma[t]$, that is, performing vertical movements in Fig. 2(b). Previous literature approaches the critical point $W_{C}$ by dynamic Levina et al. (2007, 2009); Bonachela et al. (2010); Costa et al. (2015); Campos et al. (2017) or Hebbian de Arcangelis et al. (2006); Pellegrini et al. (2007); de Arcangelis and Herrmann (2010); De Arcangelis (2012); De Arcangelis and Herrmann (2012); van Kessenich et al. (2016) synapses. This corresponds to fix $\Gamma$ and allow changes in $W_{ij}[t]$ along the horizontal axis, see Fig. 2(b).

The two homeostatic strategies are similar, but we stress that we have only $N$ equations for the gains $\Gamma_{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 Levina et al. (2007); De Arcangelis (2012) 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 $\Gamma^{*}$ 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 $\tau$. Previous claims about achieving exact SOC by using dynamic synapses are based on the erroneous assumption that we can use a synaptic recovery time $\tau\propto N^{a}\rightarrow\infty$ Levina et al. (2007); Bonachela et al. (2010); Costa et al. (2015); Campos et al. (2017). But if we use a finite $\tau$, which is not only plausible but biologically necessary, we obtain SOSC, not SOC Campos et al. (2017); Brochini et al. (2016). However, we found that for large but plausible values of $\tau$, 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. Friedman et al. (2012) and Fig. 4 of Scott et al. Scott et al. (2014). 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 Priesemann et al. (2009); Girardi-Schappo et al. (2013); Levina and Priesemann (2017), so the analysis must be done with some care.

Supercriticality, in the form of the so called Dragon-king avalanches Sornette and Ouillon (2012); De Arcangelis (2012); Lin et al. (2017), has been conjectured to be at the basis of hyperexcitability in epilepsy Hesse and Gross (2015); Hobbs et al. (2010); Meisel et al. (2012). Also, networks can be put artificially in the hyperexcitable state and show bimodal distributions $P_{S}(s)$ with large supercritical bumps Haldeman and Beggs (2005). The SOSC phenomenon seems to be a natural explanation for such hyperexcitability. In Haldeman and Beggs (2005), 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 $\tau$ recovery time.

The unexpected oscillations in $\Gamma[t]$ around $\Gamma^{*}$ have amplitudes that depend on $\tau$ and vanish for large $\tau$ (Fig. 5(d) and Fig. 7(a)). These oscillations in $\Gamma[t]$ induce oscillations in the activity $\rho[t]$ (Fig. 7(b) and 7(c)). In our model, the discrete time interval $\Delta t$ is postulated as describing the width of a spike, that is, $\Delta t\approx 1-2$ ms. From our simulation data, with these values for $\Delta t$, we obtain frequencies $f\approx 0.5-16$ Hz, depending on the $\tau$ value.

Interestingly, this frequency range includes Delta, Theta and Alpha rhythms. Coexistence of Theta waves and neuronal avalanches has been observed experimentally Gireesh and Plenz (2008). Also, some theoretical work recently discussed coexistence of oscillations and avalanches Poil et al. (2012).

The presence of oscillations can mean that the fixed point $(\Gamma^{*},\rho^{*})$ is unstable below some biffurcation point $\tau_{b}$ (even for $N\rightarrow\infty$) or that it is stable but has a very small negative Lyapunov exponent, such that finite-size fluctuations drive $\Gamma[t],\rho[t]$ away from equilibrium, producing excursions (oscillations) in the $(\Gamma,\rho)$ 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 Levina et al. (2007); Bonachela et al. (2010); Costa et al. (2015), 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. 5(b) and 5(d)) 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 Bedard et al. (2006); Tetzlaff et al. (2010); Priesemann et al. (2014). How we can conciliate these findings? Here we offer an answer based on the findings of Priesemann et al. Priesemann et al. (2014). These authors found that, in order of explain in vivo experiments with awake animals, they need three ingredients: subsampling Priesemann et al. (2009), 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 $\Gamma_{i}[t]$ in Eq. (23), so that a new equilibrium is achieved with $\Gamma^{*}<\Gamma_{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.

6 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 $\Phi(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.

7 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 Campos et al. (2017). 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 Sornette and Ouillon (2012) 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 De Arcangelis (2012).

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 $\mu>0$.

Acknowledgements.

This article was produced as part of the activities of FAPESP Research, Innovation and Dissemination Center for Neuromathematics (grant #2013/07699-0, S.Paulo Research Foundation). The RIDC for Neuromathematics covered publication costs. AAC also thanks grants $\#$2016/00430-3 and $\#$2016/20945-8 São Paulo Research Foundation (FAPESP). LB also acknowledges grant $\#$2016/24676-1 São Paulo Research Foundation (FAPESP). OK also received support from Núcleo de Apoio à Pesquisa CNAIPS-USP.

\authorcontributionsAAC performed the network simulations, and prepared all the figures. OK and LB performed analytic and numerical calculations. All authors analyzed the results and wrote the manuscript.

Conflicts of interest:

The authors declare no conflict of interest.

Appendix A Phase transition for $\mu>0,V_{T}=0$

We want to derive the critical point for the leakage case $\mu>0$ and also to obtain approximate curves for the activity $\rho$ near the critical region. We start from the exact formulas (supposing $\Phi(0)=0$):

[TABLE]

then, we use the recurrence relations Eq. (31) and (33) into Eq. (30):

[TABLE]

We notice that, due to Eq. (33), all terms $U_{k}$ are small in the critical region where $\rho\rightarrow 0$. So, we approximate the rational $\Phi(U)$ function for small $U$:

[TABLE]

Inserting in Eq. (36), and using Eq. (33) we get:

[TABLE]

Since each $U_{k}$ is proportional to $\rho$, from now we conserve only terms proportional to $\rho$ and $\rho^{2}$. After recombining the terms up to order $V_{k-1}^{2}$ according to Eq. (37), we obtain:

[TABLE]

Notice that $\sum_{k=1}^{\infty}\eta_{k-1}=\sum_{k=0}^{\infty}\eta_{k}=1$ by normalization. Using this fact and also Eq. (30) (using $\Phi(U_{0})=0$), after some rearrangement we obtain:

[TABLE]

With respect to the last term, we use Eq. (34) and (30) to obtain

[TABLE]

which is valid for $\mu<1$. Here, the sum $\rho\sum_{k=1}^{\infty}\eta_{k-1}\phi(U_{k-1})\mu^{k}$ is composed of terms in $\rho^{3}$ than can be dismissed. Using this approximation in Eq. 39 we obtain two solutions. One is the absorbing state $\rho=0$. The other solution is

[TABLE]

where we considered $\Gamma_{c}=(1-\mu)/W$. Moreover, in the critical region we can approximate $\Gamma W\approx 1-\mu$, leading to:

[TABLE]

We compare this analytical approximation with numerical solutions for $\rho(\Gamma W,\mu)$ near the critical point, see Fig. 1(c).

A similar calculation for the monomial function $\Phi(V)=\Gamma V\Theta(V)\Theta(\Gamma V-1)+\Theta(1-\Gamma V)$ gives:

[TABLE]

with the same critical line $\Gamma_{C}=(1-\mu)/W$. The monomial function with $\mu>0$ was studied numerically in Brochini et al. Brochini et al. (2016) but this analytic proof for $\Gamma_{C}$ is new.