Self-organized bistability: is it relevant for brain dynamics?

Victor Buendía,1, 2, 3 Serena di Santo,4, 5 Pablo Villegas,6 Raffaella Burioni,2, 3 and Miguel A. Muñoz1, 2 Departamento de Electromagnetismo y Física de la Materia e Instituto Carlos I de Física Teórica y Computacional. Universidad de Granada, E-18071 Granada, Spain Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, via G.P. Usberti, 7/A 43124, Parma, Italy INFN, Gruppo Collegato di Parma, via G.P. Usberti, 7/A 43124, Parma, Italy Cognitive Neuroscience, SISSA International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy Departamento de Electromagnetismo y Física de la Materia e Instituto Carlos I de Física Teórica y Computacional. Universidad de Granada. E-18071, Granada, Spain Istituto dei Sistemi Complessi, CNR, via dei Taurini 19, 00185 Rome, Italy


I. INTRODUCTION
The theory of self-organized criticality (SOC) explains how systems can become self-organized to the edge of a continuous phase transition, i.e. to the vicinity of a critical point without the apparent need of parameter fine tuning [1][2][3][4][5]. This theory (or mechanism) is often invoked to explain the emergence of scale-free distributions in the sizes and durations of bursts of activity interspersed by periods of quiescence -often called "avalanches"-in very diverse settings from earthquakes to vortices in superconductors [1,3,5], as well as in cortical brain activity [6][7][8][9].
The basic mechanism for self-organization to the edge of a phase transition -as exemplified by its most paradigmatic representatives, sandpile models [2,10,11]-relies on two essential and intertwined features: The first one is the presence of two infinitely separated timescales: a fast dynamics characterizes the system intrinsic activity, while a slow dynamics modulates the control parameter. This dynamics of the control parameter acts differentially in each phase of the system state, thus creating a feedback loop that self-organizes the system to the edge of the transition separating both phases [5,[12][13][14]). The second important feature is that the bulk dynamics is conserved, i.e. dissipation occurs only at the system boundaries; note, however, that SOC models lacking conservation (as exemplified by forest-fire and earthquake models [15,16]) have also been long considered. The main difference between conserved and non-conserved dynamics in SOC is that while the first one drives the dynamics exactly to the critical point, non-conserved dynamics leads to a wandering around the critical point (i.e. with excursions to both sides of the critical point but not sitting exactly on it; we refer the reader to Bonachela and Muñoz [14] for a detailed and pedagogical discussion of these issues).
The paradigm of self-organized criticality has been very influential in the context of neuroscience to rationalize the empirical observation of scale-free avalanches of neuronal activity, whose sizes and durations are distributed as power laws [17], as robustly observed across species, scales, and experimental techniques [18,19]. Mathematical models inspired in SOC have been proposed to account for such scale-invariant avalanches [7,8,20,21]. Actually, the idea that the cerebral cortex -as well as other biological systems-might extract important functional benefits from operating at criticality, has attracted a lot of interest and excitement as well as some controversy (see e.g. [9,18,19,22,23]).
Recently, a mechanism similar in spirit to SOC -relying also on an infinite separation of timescales and on a conserved dynamic for the control parameter-has been discovered to be able to self-organize systems exhibiting a discontinuous phase transition to the edge of such a transition. As a matter of fact, self-organized bistability (SOB) -as the theory has been named-has been proposed as a new and very general paradigm for the selforganization to points of bistability or phase coexistence between two alternative phases.
SOB dynamics turns out to be also characterized by avalanches of activity, much as in SOC. In other words, localized perturbations to the quiescent state propagate in the form of scale-free avalanches with sizes S and durations T distributed as P (S) ∼ S −τ and P (T ) ∼ T −α ; and, also, the averaged avalanche size scales as S ∼ T γ , obeying the scaling relation γ = α−1 τ −1 . However, at odds with SOC, SOB avalanches are distributed in a bimodal fashion, i.e. for any finite system size their sizes and times distributions consist of a power-law complemented with a "bump", corresponding to anomalously large systemspanning events [24,25]. Moreover, the exponents τ, α and γ differ from their SOC counterparts and coincide with those of the usual mean-field branching process exponents [24,[26][27][28].
Similarly to SOC, this new paradigm has been argued to be of potential relevance in biological problems, given that living systems could have evolve to exploit the complementary benefits of two alternative coexisting phases [9].
In the original paper in which SOB was introduced it was also speculated that SOB might be relevant to explain the emergence of scale-free avalanches in the neuronal activity of the cerebral cortex [24] . Actually, during deep sleep or under anesthesia, the state of cortical activity is well-known to exhibit a form of bistability, with an alternation between two possible states of high and low levels of neural activity -called "up and down" states-, respectively [29][30][31]. This underlying bistability, together with the empirical observation of scale-free avalanches -appearing sometimes in concomitance with anomalously large outbursts-in the awake resting brain, made the authors of Ref [24] suggest (as a possibility to explore) that there could be some type of self-organization to the edge of bistability, rather than the usually postulated self-organization to criticality [7,8,20,21].
In a parallel endeavor, our research group has recently proposed a physiologically-motivated mesoscopic (Landau-Ginzburg) theory, specifically designed to shed light on the large-scale dynamical features of cortical activity [32]. The outcome of such an approach is that the relevant phase transition for cortical dynamics is a synchronization phase transition (which occurs in concomitance with scale-invariant avalanches). Such a synchronous-asynchronous phase transition is of a different nature of the quiescent-active phase transition assumed to describe avalanches in cortical dynamics. This observation has motivated a change of perspective in the field with important consequences.
The most remarkable fact for our purposes here is that the Landau-Ginzburg theory for cortical dynamics bears profound similarities with that of SOB. The main goal of the present paper is to scrutinize the analogies and differences between the above Landau-Ginzburg theory and the theory of SOB and to determine whether, actually, self-organized bistability might play any relevant role in neuroscience as initially conjectured.

II. THE THEORY OF SELF-ORGANIZED BISTABILITY
Let us overview here the main aspects of the theory of SOB (for the sake of clarity and self-containment, we present in Appendix A a simple one-site or mean-field formulation that might be helpful to gain insight for readers not familiar with SOB or SOC theories; for a complete description of SOB and its mathematical formulation we refer to di Santo et al. [24]).
Given a spatially-extended system, the theory of SOB can be written in terms of two equations, one for an activity field ρ(x, t) (which in sandpiles represents sites over threshold) and one for a "background" or "energy" field [33].
More precisely, the set of Langevin equations proposed to describe the evolution of the activity ρ( x, t) and "energy" field E( x, t) in the SOB theory [24], reaḋ (note that some dependencies on ( x, t) have been omitted in Eq.(1) for simplicity), a, b > 0 are constants, ∇ 2 ρ stands for a diffusive coupling with diffusion constant D, η( x, t) is a zero-mean multiplicative Gaussian noise with η( x, t)η( x , t ) = ρ( x, t)δ( x − x )δ(t − t ) describing particle-number ("demographic") fluctuations. This is the simplest set of equations extending the Langevin set of equations for self-organized criticality [5,14,34,35] to the case in which there is a discontinuous phase transition [36]. More in particular, the rationale of these equations is as follows. The first one for the field activity exhibits a discontinuous phase transition for the overall level of activity as the control parameter (e.g. value of E) is changed. The second one is responsible for creating a feedback loop between the control and the order parameters [12]. When there is no activity, the local energy field grows with the constant driving +h (shifting the system towards the active phase) while, in the presence of activity, dissipation at the boundaries dominates (shifting the system towards the absorbing phase). Importantly, h is an arbitrarily small driving rate, which, in analogy with the driving in sandpile models -where just one sandgrain is added at the time [1]-is assumed to scale with the inverse of the system size h ∝ 1/L 2 . Equivalently ε is the rate of the activity-dependent energy dissipation, which in sandpile models is a boundary effect, scaling as ε ∝ 1/L; thus, the ratio h/ε → 0 in the large system-size limit. In other words, the SOB limit is obtained when driving and dissipation goe to zero, i.e. there is conservation and the time scales of both equations are infinitely separated (in the jargon of non-linear systems this corresponds to a "slow-fast" dynamics).
As mentioned in the introduction, the outcome of this set of equations is that the system self-organizes to the point of bistability between the quiescent and the active phase with scale-free avalanches [24]. In particular, the associated exponents τ and α and γ take wellknown values, i.e. those of an unbiased branching process [26,37,38] (actually, to be more precise, they coincide with those of compact directed percolation or voter model universality class in all dimensions [39][40][41], and with those of the branching process above the upper critical dimension d = 2 [24]). Also, as said above, such scalefree distributed avalanches coexist (for any finite system size) with anomalously large events or waves, in which activity spreads ballistically all through the system [24].

III. THE LANDAU-GINZBURG THEORY OF CORTEX DYNAMICS
The set of Eqs.(1) can be directly compared with the physiologically-inspired Landau-Ginzburg model that has been recently proposed to shed light on the largescale features of brain activity, relying on synaptic plasticity as a chief regulatory mechanism [32]. In this modelling approach, neural activity is described at a coarsegrained level in the spirit of the approach of Wilson and Cowan to large-scale neural dynamics [42,43]. On the other hand, short-term synaptic plasticity is implemented as a main regulatory mechanism -using the celebrated Tsodyks-Markram model [29]-in line with previously proposed models of self-organization in brain dynamics [7,8,44,45]. In particular, the level of neural activity at each coarse-grained or mesoscopic region of the cortical tissue (representing a local subpopulation of neurons) is encoded in an activity variable ρ( x, t), while the local amount of available synaptic resources is called R( x, t). The following set of equations defines the dynamics of these two variables [32]: (where, again, some dependences on ( x, t) have been omitted for the sake of simplicity). In the equation for the activity field -much as abovea, b > 0 are constant parameters, I stands for a small external incoming input, ∇ 2 ρ stands for the diffusive coupling between local regions with diffusion constant D, and η( x, t) is a zeromean Gaussian white noise. The bilinear coupling between R and ρ in the first equation reflects the fact that the larger the amount of synaptic resources, the larger the rate at which further activity is generated. On the other hand, in the second the equation, τ R and τ D are the characteristic scales for the processes of recovery and depletion of synaptic resources, respectively, ξ is the maxi-mal level of available synaptic resources that can possibly be reached.
As discussed in detail in [32], let us emphasize that -importantly for what follows-this theory is not selforganized, i.e. the resulting phenomenology depends on on the value of the baseline level of synaptic resources, ξ: • For small values of ξ the system sets into the quiescent (or "down") state.
• For large values of ξ the system sets into an active (or "up") state with a sustained (high) level of activity, which is called "asynchronous irregular" phase.
• In between the above two regimes there is an intermediate one called "synchronous irregular" in which there are intermittent events which are synchronized; these correspond to fast waves traveling through the system, and generating co-activation of many units within a relatively small time window.
• Right at the critical point separating the synchronous from the asynchronous phase avalanches can be measured by using a protocol -relying on the definition of a discrete time binning-identical to that used by experimentalists to detect neuronal avalanches [9,17].
At the light of these facts, it was proposed that the actual scale-free avalanches of neuronal activity observed in the cortex stem from the dynamics operating in a regime close to the synchronous-asynchronous phase transition, rather than to the critical point of a standard quiescentto-active phase transition as had been assumed before [18]. However, the mechanism by which the cortex seemingly selects to operate near such a critical point remains undetermined [32].

IV. ANALOGIES AND DIFFERENCES BETWEEN BOTH THEORIES
Remarkably, the sets of equations Eq.(1) and Eq.(2) exhibit profound analogies (note that the notation in both of them has been fixed to make the similitudes self-evident). As a matter of fact, the first equation, which describes in the simplest possible (minimal) way a discontinuous/first-order phase transition into an absorbing state, is identical for both theories. On the other hand, formal differences appear in the second equation for the "energy" field in SOB theory, Eq.(1), which in Eq.(2) is replaced by a field describing the dynamics of synaptic resources.
A first important difference between both sets of equations is that in SOB the energy field is conserved in the bulk; only driving events and boundary dissipation make the total integral value fluctuate in time, but as specified earlier, they are both vanishingly small. On the other hand, the above equation forṘ( x, t) is not conserved in the bulk, as it includes a (positive) term for the charging/recovery of resources to their baseline level ξ, as well as a (negative) term for the activity-dependent consumption of resources.
The second (related) difference is that in Eq. (2) there is no perfect separation of time scales. Actually, the synaptic time constants τ R and τ D are biologically motivated parameters, e.g. τ R = 10 3 and τ D = 10 2 time units [32], that cannot be assumed to scale with sytem size (as would be required in SOB) [21].
In third place, in Eq. (2) there is a baseline level of synaptic resources ξ which plays the role of an overall control parameter: by tuning it one can shift the system between different phases from a quiescent regime with no activity, to one with oscillations, waves etc. (see [32]), and this has no counterpart in Eq.(1) for SOB. Indeed, the very existence of a control parameter allowing to shift between dynamical regimes is at odds with the concept of self-organization.
Thus our aim in what follows is to explicitly explore whether Eq.(2) can lead to the same behavior of SOB -i.e. to true self-organization-in some limit.

V. RESULTS
Let us start with a single unit (or mean-field) analysis of the local/individual components of Eq. (2) in the stationary state, for which one needs to detect where the two nullclines of the associated deterministic equations (for the one-site dynamics) intersect (see Appendix B for more details). The fixed points (ρ * , R * ) of the dynamical system depend on the maximal level of synaptic resources, ξ. For small values of ξ, the deterministic system settles in a "down" state where the activity (ρ * ) vanishes and R * = ξ. In contrast, an "up" state with sustained activity (ρ * > 0) with a depleted level of resources (R * < ξ) emerges for sufficiently large values of ξ.
Separating these two limiting (down and up) states, there is a range of values of ξ where a stable limit cycle emerges if the timescale separation is large enough (see leftmost inset of Fig.1A and purple-shaded region in Figure 1). As discussed in Appendix B, the limit cycle is created via a saddle-node-infinite-period bifurcation at ξ = a, and destroyed via a Hopf bifurcation at a value of ξ that depends on the timescale. In particular, for a plausible separation of the characteristic scales of synaptic resources, τ R and τ D as chosen in [32] (i.e. τ D /τ R = 0.1 and τ R = 10 3 ) the system exhibits an intermediate regime with oscillations as illustrated in Figures  1A and A1. These oscillations at the individual sites are at the origin of the propagating waves and emerging synchronous behavior in the extended system [32].
On the other hand, in the opposite limit in which the separation of timescales are both relatively fast and not very separated -a case that we do not discuss here in further detail-there is no intermediate regime of oscillations and the system just shifts from a quiescent (down) to an (2). Insets: the nullcline forρ = 0 is colored in black, while nullclines forṘ = 0 -for three different values of the baseline level of synaptic resources ξ-are plotted in orange, purple and green color, respectively. The small grey arrows represent the vector field for (ρ,Ṙ) and the light purple curve represents a limit-cycle trajectory. Inset (A) For τD/τR = 0.1 and τR = 10 3 the system may display a down state (orange nullcline, ξ = 0.2), a limit cycle (purple nullcline, ξ = 1.5) or an up state (green nullcline, ξ = 3.5). Inset (B) For τD/τR = 0.001 and τR = 10 7 the system approaches the SOB behavior and the dependence on the control parameter ξ becomes very weak, meaning that for a wide range of values of the control parameter above the down state (ξ = 1, 5, 10 for orange, purple and green nullclines respectively) the system sets in the oscillatory phase. Inset (C) When the separation of timescales is very low, the system only displays bistability between up and down states with no oscillations whatsoever, no matter the value of ξ. Main The global behaviour of the system depending on the timescale separation is coded in the main figure with different background colors: purple when the limit cycle appears, and green when only bistability is possible. Parameters are set to I = 10 −3 , a = 0.6, b = 1.3, in both cases.
active (up) state as ξ is increased (see the green-shaded region of Fig.1 and Fig.1c). This regime is likely to be useful to describe up-and-down transitions in the cortex.
Given that the mechanism of the recovery of synaptic resources (with a characteristic timescale τ R ) plays the role of driving, while the depletion of resources (with a characteristic scale τ D ) in presence of neural activity plays the role of dissipation, we conjecture that the SOB limit of infinite separation of timescales could be only possibly recovered taking formally the double limit 1/τ R → 0 and 1/τ D → 0, with τ D /τ R → 0. In particular, in what follows we consider different sets of values for  16 . The dashed lines are plotted as a guide to the eye, and have slopes corresponding to the expectations for an unbiased branching process (α = −2, τ = −3/2 and γ = 2, respectively). The "bumps" in the blue curves, correspond to anomalously large events, i.e. synchronized spiking events. The cut-offs/bumps change with N obeying finite-size scaling as in the theory of self-organized bistability [24]. Parameters: b = 0.5 a = 1, I = 10 −7 , τD = 10 4 , τR = 10 6 , D = 1. First of all, we observe that the larger the separation of timescales, the weaker the dependence on the control parameter ξ (see in particular in Fig. 1B how the three differentṘ = 0 nullclines -for three different values of ξbecome very close to each other). This allows us to reproduce one of the main features of the self-organization mechanism: the independence on any tuning parameter as the limit of infinite separation of timescale (see Appendix B for a more detailed description).
Second, when the separation of timescales tends to infinity, the slope of the nullclineṘ = 0 converges to zero at the intersection point (see Figure A1). Then, for sufficiently large values of ξ both nullclines intersect in a tangential way, giving rise to a limit cycle, much as in the mean-field SOB (see Appendix A and Fig. A1).
To go beyond the single-site or mean-field analysis we now study computationally a full (spatially explicit) system. Under the conditions of a large separation of timescales we have performed computer simulations of the full set of Eqs. (2). In particular, we considered square lattices with sizes from N = 2 12 to N = 2 16 to analyze finite size effects and performed up to 10 8 runs for different parameter sets. It is possible to measure avalanches of activity by defining a small threshold (e.g. θ = 10 −6 ) for the overall (integrated) activity and analyzing the statistics of excursions (sizes (S) and times (T)) above such a threshold [24,38]. Results of our extensive computer simulations are shown Fig. 2 and Fig. 3.
In particular, Fig.2 indeed reveals the existence of scale-invariant episodes of activity -whose distributions are very well fitted by the exponents of the unbiased branching process-appearing together with anomalous "king" avalanches in full analogy with the theory of SOB [24]. These computational results are very robust to changes in the control parameter ξ as illustrated in Fig.  3. Actually, hardly any difference is observed in the probability distributions when ξ is increased from a value 1 to 10. This reveals the emergence of true self-organization regardless of the specific value of the parameters in the limit in which an infinite separation of timescales is imposed.
Thus, at the steady state, the system with infinitelyseparated timescales actually self-organizes to the edge of a discontinuous/first-order phase transition, where active and absorbing phases coexist as in SOB.
Finally, observe that the lack of a diffusion term in the equation for synaptic resources does not seem to be a problem, indicating that diffusion in the second equation is not essential ("relevant" in the jargon of the renormalization group) feature and that possibly it could also be removed from the minimal set of Eqs.(1) describing SOB (this claim is also supported by our own computational analyses).

VI. DISCUSSION
We have shown, by using the recently proposed Landau-Ginzburg model of cortex dynamics [32], that self-organizing behavior to the edge of a discontinuous phase transition with bistability can be recovered by taking extremely slow synaptic timescales.
First, from a theoretical point of view, our analyses reveal that different regulatory mechanisms for the control parameter, -i.e. different equations for the "energy" field, with different meanings and possibly with diverse features such as conserved dynamics or not, or with or without a diffusion term-can be considered. In all cases, the SOB phenomenology -with scale-free avalanches controlled by branching process exponents, appearing together with anomalously large avalanches or waves-emerges once the limit of infinite separation of timescales is taken.
Second, this approach allows us to illustrate that, even if we can use diverse models to recover scale-free distributions of neuronal avalanches, the synaptic time scales needed to recover them in the limit of SOB are not compatible with those of a realistic short term plasticity dynamics. As reported by neurophysiological measurements, such scales are comprised between few hundreds of milliseconds up to few seconds [29,46,47]. Such timescales are outside the proposed SOB limit, which requires much slower dynamics. For instance, the unit of time, in the Landau-Ginzburg model, should be understood in terms of the spontaneous decay of neural ac-tivity: one time unit in Eq (2) is of the order of the millisecond, which means that in order to recover the SOB limits, the recovery timescale for synaptic resources should be of the order of several minutes! On the other hand, the Landau-Ginzburg model -leading to verifiable predictions-is based on the consideration of realistic timescales.
Second, as it has been pointed out [26,27,38], the exponents of the unbiased branching process appear in a number of diverse circumstances; but, reproducing them does not constitute a conclusive warranty of success for a given modelling approach or theory. Thus, further work is still needed to develop a clear framework of selforganization in neural dynamics beyond the paradigm of absorbing/active phase transitions and understand in detail why such values can emerge also in the context of synchronization phase transitions [32].
All of this leads to the conclusion that, in this approach, self-organization to the edge of bistability can be achieved only by considering an infinite separation of timescales that does not seem to be plausible in the cortex, at least not by considering short-term plasticity as a chief regulatory mechanism. The same criticism can be made -and has been explicitly made [21]-to models for neural dynamics based on SOC that also require an infinite separation of timescales to work.
Observe that the original Tsodyks-Markram model for short-time synaptic plasticity implements an additional variable describing a "facilitation" mechanism whose timescale could be much longer than the depression one considered here [29,48]. This would make the dynamics more complex and eventually regulate the baseline of synaptic resources, allowing self-tuning to the transition point. Also, one could explore the role played by other types of plasticities, including facilitation, such as spike-timing dependent plasticity. In particular, plasticity operates in the cortex through a wide spectrum of timescales. For example, long term potentiation is known to involve different mechanisms, such as the production of new ion channels, modifying the excitability of the units for a long-lasting period. This research line seems to be a very relevant in neuroscience and we leave a careful exploration of this possibility for a future work. In any case, adding dynamical variables to the system improves the realism of the modelling, but makes it less straightforward to establish a direct comparison with the theory of SOB.
Summing up, the SOB model can be seen as a limiting case of the Landau-Ginzburg model for cortex dynamics. It is still able to reproduce the presence of avalanches -using a definition of avalanche which is more standard for computational models but less related to experimental procedures-with sizes and durations distributed as power laws with the same exponents as observed empirically in neuroscience; however this only happens once the unrealistic limit of an extremely slow dynamics for synaptic resources is considered.
Thus, which theory is the most appropriate one can not be solely decided on the basis on the existence of avalanches and their critical exponent values. In future work, we will extend the Landau-Ginzburg theory of the cortex dynamics in a number of ways, exploring realistic ways of self-organization in the absence infinitely separated timescales, as well as analyzing the theoretical causes and origins of the resulting exponents, to pave the road for a more direct connection with biological systems.

ACKNOWLEDGMENTS
We acknowledge the Spanish Ministry and Agencia Estatal de investigación (AEI) through grant FIS2017-84256-P (European Regional Development Fund), as well as the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOM M 17/6105/U GR for financial support. We also thank Cariparma for their support through the TEACH IN PARMA project.

Appendix A: Mean field SOB
Here we perform a detailed single-site or mean-field analysis of the dynamical system describing SOB, i.e. Eqs.(1) but neglecting spatial dependence and noise. The first equation in Eqs.(1) can then be simply written aṡ where now ρ is the overall density of activity. Let us analyze this equation forĖ = 0, i.e. for a fixed value of E, with vanishing driving h and dissipation ε. Then Eq.(A1) has the typical form of a discontinuous transition for ρ. Actually, imposingρ = 0, one obtains three possible solutions: ρ 0 = 0 (stable for E < a) and a positive and negative pair, that emerge at saddle-node bifurcation. Indeed, as illustrated in Fig.A.1, these last ones exist (as real solutions) only for E > a−b 2 /4 and ρ + is stable for a−b 2 /4 < E < a (ρ − is unstable in this region, and becomes non-physical, i.e. negative for E > a). This gives a bistability region (shaded in magenta color in Fig.A.1) where two stable fixed points ρ 0 and ρ + coexist, with a branch of unstable fixed points ρ − in between. These are the solutions of the problem for constant "energy" E. For the full system, including the dynamics for E (i.e. h = 0, ε = 0) the mean-field dynamics of the control parameter E is given by: Observe that this dynamics acts differentially depending on the phase, i.e. the value of ρ: in the absence of activity it increases owing to the driving h, while in the presence of activity there is an additional dissipation (negative term). Observe also that the only possible fixed point -where these two contributions balance-is ρ * = h/ε. Importantly, this equation defines a nullcline which is a horizontal line in Fig.A.1. Note also that in the double limit of conservation h, ε → 0 with an infinite separation of timescales h/ε → 0, this is necessarily a point of vanishingly small activity and, as seen above (see Fig.A.1), the branch of points with low activity is unstable. Indeed, from eq.(A2) at stationarity one obtains E * = a − bh/ε + (h/ε) 2 , which corresponds to an unstable spiral focus. Inspection of the associated velocity field, (see Fig.A.1) it becomes clear that the dynamics exhibits a stable limit cycle around such an unstable fixed point, shifting cyclically between the upper and the lower branches, with high and low activity, respectively. Note that in the limit h/ε → 0, the timescale between succesive cycles increases (in particular, the imaginary part of the eigenvalues associated with the unstable spiral (ρ * , E * ), scales as ∼ b/2h/ε, meaning that the frequency of oscillations goes to 0 as happens in an infinite-period bifurcation [49]). Thus, this simple mean-field approach shows that the system behaves as an excitable system in which upon Behaviour of the Ginzburg-Landau system as a function of the control parameter ξ, as the limit ∆ → 0 is approached. Each color represents a different phase (see legend). Note that as the separation of timescales increases (i.e. as ∆ → 0) the oscillatory phase takes more space until it comprises the whole region ξ > a (SOB limit). Parameters are a = 0.6, b = 1.3.
driving the system it exhibits cycles of periodic activity and large amplitude. As shown in detail in [24], going beyond mean-field -i.e. turning back to a spatially explicit system such as a lattice-this mechanism generates scale-free avalanches alternating with system-wide spanning waves of activity. The reason behind such a change is that in the above situation driving may induce local jumps of the activity from the lower branch to the upper one, and from this it can propagate to nearest neighbors giving rise to an avalanche of activity, before falling down again (owing to dissipation) to the quiescent state [24].

Appendix B: Mean field Ginzburg-Landau
Here we discuss the single site or mean-field description of the Landau-Ginzburg theory for cortex dynamics. The first single-unit equation is identical to Eq.A1 while the second one readsĖ where E is now the overall "energy" or density of synaptic resources. We call it here E rather than R as in Eq.2) to make a more direct comparison with Appendix A; also to make the parallelism with the SOB theory even more explicit, we define h = 1/τ R and ε = 1/τ D , and study the limit ∆ ≡ h/ε → 0 in which the two timescales are infinitely separated.
The main difference between this mean-field theory and its SOB counterpart (in Appendix A) is the existence in Eq.(B1) of a key parameter ξ, which bounds the largest possible value of E.
Three cases can be distinguished depending on the value of ξ and the separation of timescales.
(i) If ξ < a obviously E < a, and then equation (A2) has a negative linear term, implying that ρ = 0 is the only possible stable state. In other words the two nullclines intersect in the absorbing phase, i.e. for ρ = 0 (see yellow line in the leftmost inset of Fig.1). Actually, this fixed point that can be easily obtained analytically, ρ = 0, E = ξ, which loses its stability for E ≥ a (the point E = a marks a saddle-node infinite-period transition that gives birth to a limit cycle).
(ii) More in general, imposingĖ = 0, one finds the fixed point value E * = hξ (h+ε)ρ (defined only for ρ > 0) which in the limit of large time-scale separation (h ε) becomes E * ∆ξ ρ . This implies that the nullclines intersect at an unstable fixed point (see red line in the leftmost inset of Fig.1). This leads to oscillations in the following way: within the absorbing state, energy slowly increases as h(ξ − E) and -if ξ > a-at some point the system reaches the value E = a, where the absorbing state ρ = 0 becomes unstable. When this happens, the energy fixed point is given by E * . If the energy is larger than this value, it decreases in a fast way, until E = E * . However, if ∆ is very small, the fixed point is E * 0, meaning that E decreases until the up-state branch becomes unstable again -falling to the absorbing state. When the absorbing state is reached, E * is unstable and the cycle starts again. This simple back-and-forth mechanism is responsible for the emergence of limit cycles in the dynamics for ξ > a (see also [50] where this is termed "Sisyphys effect".).
(iii) For larger values of ξ, such that the nullclines intersect in the (stable) up-state branch (see the brown line in the leftmost inset of Fig.1) there are no oscillations and the system sets in an actie state with E * = ∆ξ ρ (the oscillations disappear at a Hopf bifurcation). Observe that the condition for this case (iii) to emerge depends on the separation of timescales, ∆; in particular, in the limit ∆ → 0, a value ξ → +∞ is required to destabilize the cycle limit. This effect is illustrated in Figure B.2, where it can be seen that as the ∆ → 0 limit is approached the regime of oscillations broadens (up to infinite).
Finally, let us note that in case (ii) -which is the one for which avalanches are obtained once the spatially explicit version of the model is considered-the nullcline associated with Eq.(B1) becomes very similar to its counterpart for the SOB case: both of them intersect the unstable branch of ρ solutions in the same part and with very similar slope. The corresponding slopes are flat in both cases in the limit ∆ → 0. Also, observe in the upper inset of Fig.1 that the nullclines are quite insensitive to the specific value of ξ (as long as ξ > a) when such a limit is taken. Thus, the same phenomenology is expected to emerge in both cases.