{'fraction_non_alphanumeric': 0.043861807856128726, 'fraction_numerical': 0.010922858495030762, 'mean_word_length': 4.668025751072961, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-normality, reactivity, and intrinsic stochasticity in neural dynamics: a non-equilibrium potential approach

Intrinsic stochasticity can induce highly non-trivial effects on dynamical systems, including stochastic and coherence resonance, noise induced bistability, noise-induced oscillations, to name but a few. In this paper we revisit a mechanism first investigated in the context of neuroscience by which relatively small demographic (intrinsic) fluctuations can lead to the emergence of avalanching behavior in systems that are deterministically characterized by a single stable fixed point (up state). The anomalously large response of such systems to stochasticity stems (or is strongly associated with) the existence of a"non-normal"stability matrix at the deterministic fixed point, which may induce the system to be"reactive". Here, we further investigate this mechanism by exploring the interplay between non-normality and intrinsic (demographic) stochasticity, by employing a number of analytical and computational approaches. We establish, in particular, that the resulting dynamics in this type of systems cannot be simply derived from a scalar potential but, additionally, one needs to consider a curl flux which describes the essential non-equilibrium nature of this type of noisy non-normal systems. Moreover, we shed further light on the origin of the phenomenon, introduce the novel concept of"non-linear reactivity", and rationalize of the observed the value of the emerging avalanche exponents.

Appendix A Reactivity in non-normal systems 20 Appendix B Effective one-variable description: stochastic tunneling and the origin of the avalanche exponents 22

Introduction.
Noise is well-known to have a number of remarkable effects in many different physical, chemical, and biological systems, such that it can dramatically alter the predictions of deterministic approaches; noise-induced transitions, stochastic resonance, coherence resonance, and stochastic amplification of fluctuations are just a few well-acknowledged examples [1][2][3][4][5][6][7][8]. Here, we are ultimately interested in highly non-trivial stochastic effects emerging in concomitance with non-standard deterministic dynamics, i.e. with nonnormal (linearized) dynamics as described below. For the sake of specificity and clarity, we focus on a problem in the field of neural dynamics, but the ideas and mechanisms discussed here go beyond such an example and might be relevant in other apparently unrelated contexts, such as e.g. theoretical ecology and population dynamics, but, for the sake of concision, let us discuss in detail only the case of neural systems. The human brain exhibits persistent intrinsic activity even in the absence of any stimulus or task [9]. Understanding the origin and nature of such an irregular dynamics is a key challenge in neuroscience [10][11][12][13]. Pioneering experiments by Beggs and Plenz, revealed the existence of scale-invariant episodes of spontaneous electrochemical activity in neural tissues in vitro, thereafter named neural avalanches. Subsequently, neural avalanches were detected in a wide range of experimental settings, tissues and species both in vitro [14][15][16][17][18] and in vivo [19][20][21][22]. Experimental results produce rather consistently avalanche exponents compatible with those of an unbiased branching process [14,23]. Such a scale-invariant organization has been taken as an indicator that cortical dynamics might operate close to a critical state [14,[23][24][25]. Such a putative criticality might endow the system with huge sensitivity to stimuli, large spatio-temporal correlations, optimal transmission of information, and a number of other functional advantages [13,23]. Theoretical models have actually proposed a link between cortical dynamics and dynamical criticality [26][27][28][29][30][31].
Nevertheless, the so-called criticality hypothesis in cortical networks is still controversial [32,33] and some authors have highlighted that it is not clear whether the available empirical evidence actually calls for criticality or other alternative explanations -such as noise, bistability, neutral dynamics, etc-could be invoked [17,[33][34][35] (we refer to [23] for an extended overview and discussion of these issues). Thus,a careful scrutiny of possible alternative scenarios for the emergence of scale-invariant avalanches is an important general task. One of the main goals of the present paper is to contribute to this broad open problem.
In mathematical and computational models, the emergence of irregular bursts of neural activity in cortical tissues has often been related to the large excitability of the underlying dynamical system, as a result of which the intrinsic dynamics can result in large excursions from stable equilibria [2,36]. In particular, in models of neural networks, large excitability is often been explained as resulting from a delicate interplay between excitation and inhibition [37][38][39] that imposes that parameters in mathematical models need to be fine tuned to set the system close to the transition (critical) point between two distinct phases (excitation-dominated and inhibition-dominated, respectively).
Within this context, Benayoun et al. made the very interesting observation that avalanches, i.e. scale-invariant bursts of activity, may emerge in finite populations of neurons in a rather robust way, even when parameters are chosen away from the perfectly balanced (critical) condition [40]. To understand this -and inspired by previous work by Murphy and Miller [41]-they relied on a mechanism called "balanced amplification of fluctuations" [40], which is relevant for noisy neural systems when the levels of excitation and inhibition are relatively balanced, i.e. they are similar but not necessarily identical (in other words, when the system does not necessarily lie at the very critical point (see below)).
Our first aim here is to shed further light on the various factors contributing to the appearance of non-critical -but still approximately scale-invariant-avalanches of activity emerging under balanced amplification of fluctuations [40].
The structure of the paper is as follows. In Section 2 we briefly review the basic aspects of the dynamical model for neural systems, and describe in some detail the mechanism uncovered by Benayoun et al. for the emergence of avalanching behavior in finite populations. In Section 3 we revisit the results in [40], putting the emphasis of the importance of the interplay between "non-normal forms","reactive dynamics" and "demographic fluctuations" in neuronal populations. Moreover, we construct a statistical-mechanics description of the problem based on non-equilibrium potentials which helps shedding light onto the overall phenomenology. Finally, we present a broader discussion of our results.

Wilson-Cowan model for excitatory and inhibitory networks
Following Benayoun et al. [40], we consider the Wilson-Cowan mean-field description of a large-scale homogeneous population of excitatory and inhibitory neurons [42]. As usual, in this type of mean-field approaches, the connections between neurons within a large population are assumed to be dense enough so that heterogeneity and fluctuations can be neglected. In spite of its simplicity, the Wilson-Cowan model encompasses -as parameter values are changed-a plethora of possible scenarios, which strikingly resemble experimentally observed dynamical regimes of neural dynamics, such as multiple coexisting stable states (Up-Down states), oscillatory behavior, simple and multiple hysteresis loops, etc. [42][43][44][45]. We refer to Eef. [46] for a detailed analysis of the Wilson Cowan model and its rich phenomenology.
For the particular version of the Wilson-Cowan model considered in ref. [40] (sketched in Fig. 1A), the mean-field equations describing the overall dynamics of the activity (density of active neurons) for the two subpopulations of excitatory (E) and The excitatory population interacts with coupling parameter, ω E and the inhibitory population with coupling ω I (a more general case is that in which the effect of excitation on excitatory and inhibitory populations, is different from the effect of inhibition). Right: Phase diagram; the phase of the system depends on the control parameter, ie. with the difference ω 0 = ω E − ω I . If ω 0 > 0 the system is in the active or "up" phase, while for negative values it is in the low-activity "down" state. The blue and red arrows signal two specific choices of values that will be analysed in detail in what follows. inhibitory (I) neurons, respectively, read [42]: where α is the rate of spontaneous activity decay, s is the averaged incoming current which is simply the sum of all synaptic inputs, both excitatory and inhibitory, weighted by their synaptic efficacies (ω E and ω I , respectively), plus an external small constant input current h, and f (s) is a sigmoid response function: This set of equations exhibits a regime shift at ω 0 ≡ ω E − ω I = 0, from an inactive (inhibition dominated) to an active (excitation dominated) phase. This is illustrated in Fig.1 where we plot the averaged global stationary activity, Σ = (E + I)/2 as a function of ω 0 , while keeping ω s ≡ ω E + ω I constant. In particular, fixing h = 0 the regime shift Timeseries of the total excitatory and inhibitory activity Σ = (E +I)/2 in a balanced condition (ω E = 7, ω I = 34 5 ; blue arrow in Fig.1) for increasing levels of noise (from left to right, as marked in the figure legends). Above a certain noise-amplitude threshold the system reaches the inhibition-dominated state and exhibits avalanching behavior. Parameter values: α = 0.1, h = 10 −3 is a true bifurcation occurring at the "critical" transition point ω 0 = 0 ‡.
To go beyond this simple mean-field picture -describing infinitely-large neural populations-Benayoun et al. considered a large but finite population of binary neurons obeying some standard dynamical rules §. Starting from such microscopic (rate) model and employing a large-system expansion [48,49], Benayoun et al. were able to recover the above deterministic Wilson-Cowan dynamics, Eq. (1), in the infinite-size limit. But, additionally, considering next-to-leading-order corrections, they were also able to explicitly determine the stochastic term to be added to Eq. (1) accounting for the finite size of the population, leading to the set of stochastic equations where η E,I are uncorrelated Gaussian white noises, with network-size dependent amplitude σ ∝ 1/ √ N , interpreted in the Ito sense [49].

Model phenomenology: computational results
By performing computational simulations of their individual-neuron-based model, Benayoun et al. [40] observed that -when the difference between excitatory and inhibitory synaptic weights, ω 0 , is small with respect to their sum ω 0 ω s = ω E + ω I ; i.e. in a situation that is termed "balanced" [40,41], the actual dynamics departs quite dramatically from the mean-field expectations, even for relatively low noise amplitudes (i.e. for large but not infinite network sizes). In particular, and quite remarkably, for ‡ While, if h = 0, there is a rapid regime shift, but without a true discontinuity in the derivative at the transition point, i.e. without a true phase transition. § In particular, they consider a fully-connected network model of individual spiking neurons such that (i) each neuron is either active or quiescent, (ii) the probability that each quiescent neuron becomes active depends on (a sigmoid function of) its total synaptic input and (iii) each active neuron relaxes to the quiescent state at a constant rate [40] (see also [47] for a similar approach). Reasonably good fits to power-laws are obtained employing the well-known values for the return times of a random walk (see e.g. [50]); the quality of the fits improves with increasing noise, but in any case they are not perfect straight lines, revealing non-perfect scale invariance. (Right) The averaged avalanche size as a function of the duration T , also scales as predicted by the statistics of random walks, S ∼ T γ , with γ = 3/2.
parameter values (ω 0 0) such that the mean-field Wilson-Cowan equations predict a stable up state, the dynamics -in the presence of noise-turns out to spend most of the time close to the down state (with very low activity) while showing frequent bursts of activity (such as those illustrated in Fig.2).
The sizes S and durations T of such bursts or avalanches of activity turn out to be distributed as (approximated) power-laws, similar to those observed experimentally in actual neuronal networks, even if the exponent values are found to be detail-dependent rather than universal; i.e. do not necessarily coincide with the experimentally-observed branching process exponents [40]. Thus, the problem of determining the precise origin and values of these exponents remains open.
To illustrate and extend the above numerical analyses, we performed our own computer simulations but, opposite to Benayoun et al., not of the individual-neuronbased model, but directly of the set of stochastic equations, Eq.(3). For the numerical integration we used an Euler-Maruyama scheme with step size ∆t = 10 −4 (and checked that results are not very sensitive to this choice). In particular, we considered ω E = 7 and ω I = 34/5, which is a balanced case in the sense that ω 0 (1/5) is relatively small when compared with ω s (69/5). Our results are summarized in Fig.2 where we confirm the existence of a stable fixed point with small fluctuations around it for small noise amplitudes (e.g. σ = 10 −4 ). On the other hand, slightly larger values of noise amplitude (i.e. smaller system sizes) produce large fluctuations around such a stable fixed point, and, more remarkably, above some value of the noise amplitude, σ ≈ 5.10 −3 , the stability of the up state is severely compromised and the system ends up hovering around a lowactivity value with large excursions to the up state. Such excursions closely resemble avalanches of activity, growing, spreading and coming back to the inhibition-dominated down state.
The resulting avalanche-like dynamics suggests a definition of the duration T and size S of avalanches as the activity over a small arbitrary threshold, allowing us to confirm the results of ref. [40]; i.e. both avalanche duration and size follow an approximate power law of the form P (T ) ∼ T −α and P (S) ∼ S −τ with S ∼ T γ (see Fig.3). Contrarily to the case of [40] we can obtain relatively clean values for the exponents. In particular, although the size and time distributions are not perfect straight lines in the double-logarithmic plot, revealing a lack of strict scale invariance, the slopes of the best fits are compatible with the well-known exponents for the return times of a random walk, i.e. α = 3/2, τ = 4/3 and γ = 3/2, rather than those of the branching process (α = 2, τ = 3/2 and γ = 2: see e.g. [50] for a pedagogical discussion of these two classes). Let us finally remark that scale invariance is only approximate in all cases; there is no simple limit in which perfect power laws and scaling emerge.

Deterministic dynamics
To shed light onto this remarkable noise-induced effect, one can scrutinize the phase portrait of the deterministic dynamics of Eq.(1). Figure 4 shows the phase portrait plane for four different sets of parameter values. In the upper panels (corresponding to parameters signaled with a blue arrow in Fig.1) the stability is weak (ω 0 = 0.2) while in the lower ones (red arrow in Fig.1) it is stronger (ω 0 = 3). The most relevant case for our discussion here is that in Fig.4b, which corresponds to the balanced condition. One can readily observe that, in this last case, the force (or velocity) field is rather anomalous: it exhibits a sort of "scar" close to the diagonal, with shear flows in opposite directions above and below it, allowing for wild excursions far away from the fixed point if it is perturbed ¶ [40]. Similarly, in Fig.4d, even if the system is more stable, there also exists a scar allowing for such type of trajectories. On the contrary, in Figs. 4a and 4c, the stability of the system is quite strong and phase portraits are completely different. We now discuss these features step by step, starting from a linear stability analysis.

Linearized deterministic dynamics: non-normality and reactivity
To illustrate the phenomenon of non-normal dynamics (still following [40]), it is convenient to change variables to Σ = (E + I)/2, ∆ = (E − I)/2, so that Eq.1 becomes These exponents obey the usual scaling relationship γ = (α − 1)/(τ − 1) [50,51]. ¶ In other words, the vector field ( Ė , İ ) shows an abrupt jump around the diagonal, where the two nullclines are very close to each other, even if they intersect just at the only fixed point. (ω E = 3000, ω I = 2997). The equilibrium is reactive R = 0.106. The main observation to be made is that the deterministic dynamics is radically different in the cases with and without a shear flow, i.e. with or without a "scar" in the diagonal.
with θ = ω 0 Σ + ω s ∆ + h. Clearly, these equations can only have fixed points of the form (Σ * , 0), i.e. lying necessarily in the diagonal of the (E, I) phase portrait. In particular, in all the considered cases (see Fig.4) there is a "up-state" fixed point with E = I in which excited and inhibited neural popualtions coexist, and equilibrate each other.
A standard linear stability analysis around the (only) up-state fixed point leads to the Jacobian matrix where the eigenvalues are and In what follows, the spontaneous decay α and the spontaneous activation rate h are fixed to relatively small values (e.g. α = 0.1, h = 10 −6 ). In particular, this imposes that if ω 0 is small and positive, then, λ 1 and λ 2 are small [40], as corresponds to a weakly stable fixed point. Therefore, as usual, the (weak) stability of the fixed point is controlled by the distance to the transition point (at ω 0 = 0). The very structure of the Jacobian makes it clear that diagonal terms -the set of eigenvalues-do not enclose all the information about the linearized dynamics. The nonvanishing off-diagonal term (the so-called feed-forward term, ω f f , cannot be eliminated by changing variables, and establishes a clear asymmetry in the dynamics: ∆ boosts Σ but not the other way around).
Triangular matrices such as J are a particular case of non-normal forms/matrices (also called non-self-adjoint matrices), meaning that J * J = JJ * , where J * is the conjugate transpose of J. As a consequence, such matrices are not diagonalizable through a unitary transformation, or in other words, the associated basis of eigenvectors is not orthogonal [52,53] + . Non-normal matrices might be unfamiliar to many physicists, grown up with quantum mechanics, where observables are represented by Hermitian (self-adjoint) operators. However, non-normal forms have a long tradition in other realms of physics such as turbulence, where they play an important role in e.g. shear and pipe flows [54][55][56][57], in control theory [58], in ecology [59], lasers [60], and even, recently, in quantum mechanics [61,62].
In the case under scrutiny, the basis of eigenvectors in the variables (Σ, ∆) is /ω f f , leading to almost identical (degenerate) eigenvectors in the limit ξ → 0, i.e. ω f f λ 1,2 , i.e. the balance condition. Also, importantly, the scalar product of both eigenvectors does not vanish, i.e. they do not form an orthogonal basis (as illustrated in Fig.4, especially in cases b and d). + Let us remark that the change of variables employed above, (E, I) → (Σ, ∆), is what is usually called a "Schur transformation", generating a triangular matrix; indeed, the Schur decomposition is a decomposition of a given matrix A, such that A = M T M −1 where M is a unitary matrix and T is an upper triangular matrix: the associated Schur form [52].
A large value of the feed-forward term ω f f in Eq.(4) may induce a huge impact on the dynamics around a weakly stable fixed point. In particular, when the matrix J operates on a small perturbation vector along the ∆ direction, (0, ), applying the linearized dynamics once, gives (+ω f f , −λ 2 ) which results in a ("contracting") response in such a direction, but a much larger ("expanding") outcome along the Σ direction. This illustrates that a large value of the feed-forward term ω f f is able to strongly affect the linearized dynamics. Observe that this can occur even when the eigenvalues are not close to zero (i.e. the system is not necessarily close to the transition point) if the feedforward term is large with respect to their absolute value (see below). As a matter of fact, since the eigenvalues depend on the control parameters ω E and ω I only through the combination ω 0 = ω E − ω I -while ω f f only depends on ω s = ω E + ω Ithe condition for the appearance of the mechanism discussed above is the "balance condition" ω 0 ω s [41]. It is noteworthy that it is possible to quantify the degree of non-normality as the weight of the feed-forward interaction (off-diagonal elements) with respect to that of the eigenvalues. In particular, one can quantify the strength of the feedforward elements of a given triangular matrix as: where λ i are the diagonal elements (eigenvalues), and the term in the denominator is the sum of the squares of all matrix elements, m n . In our particular case, this implies that, fixing ω 0 , N N tends to grow with ω s . In other words, the balance condition ω 0 ω s implies that the strength of non-normality is large. In particular, in Fig.4 cases b and d are highly non normal (non-normality N N = 0.997 and N N = 0.706, respectively), while a and c are hardly non-normal (N N = 0.068 and N N = 10 −6 , respectively).
The anomalous behavior we have described in association with non-normal matrices is sometimes called "reactivity" in the (theoretical ecology) literature [59]. As carefully explained in Appendix A, reactivity describes the property of linear (or linearized) stable systems such that their dynamics -even if converging asymptotically to a stable fixed point-can exhibit unusually long-lasting transient behavior. In other words, the system can be strongly driven away from the fixed point (actually increasing the modulus of the perturbation vector) before converging to its steady state. Also, it is shown in Appendix A that all reactive matrices are non-normal, but the opposite is not true. In particular, in the examples of Fig.4 all stability matrices are non-normal, but only (b) and (d) are also reactive.

Linearized dynamics with noise
Thus far we have revisited results in [41] and [40] regarding the very peculiar features of the deterministic dynamics in systems with non-normal stability matrices. In order to proceed further in our analytical understanding of the avalanching phenomenon, we now study the stochastic dynamics (up to leading-order approximation) around the deterministic linearized equations, to explore the role of the so-introduced fluctuations around the deterministic fixed point. For this, we consider the set of linearized Langevin equations where the vector state q = (q 1 , q 2 ) describes deviations with respect to the deterministic fixed point (Σ * , ∆ * = 0), i.e. q 1 = Σ−Σ * and q 2 = ∆, d 1 and d 2 are non-trivial functions depending on the parameters and on (Σ * , 0). For the sake of generality, we prefer to use the following compact notationq where the vector force F , i.e. the vector with components F µ with ν = 1, 2, is F = Jq, G is the (diagonal) matrix (with constant noise amplitude evaluated at the fixed point) and η is a vector composed of two uncorrelated white Gaussian variables (as usual, repeated indices are implicitly summed over).
The above Langevin equation is equivalent to a Fokker-Plank equation [48]: where B µν = G µ γ G γ ν . Eq.8 as usual, can be written as a continuity equation: where the components of the current j are An overall stationary solution exists whenever the current is divergence free; ∂ µ j µ = 0. In general, if in the steady state the current is null, j µ = 0, detailed balance holds, and the system is symmetrical under time reversal transformations: the stationary state is at equilibrium. Otherwise, even in stationary conditions, there is a net probability current flowing through the system, the detailed balance is explicitly violated and the system is away from equilibrium.
Observe that the above Fokker-Planck equation corresponds to a multivariate Ornstein-Uhlenbeck process and, as such, can be solved exactly [48,63] even in the case of non-vanishing stationary probability current [64][65][66]. However, in order to illustrate the general procedure that shall be used afterwards (for more complex, non-linearized cases), in what follows we present an analytical approach based on path integrals and a weak-noise approximation to solve it.
Employing a simple (customarily employed) ansatz for the stationary solution in the weak-noise limit [67][68][69][70][71], one can write where C is a normalization constant and V (q) a scalar potential, one can rewrite Eq.(10) -up to leading order in σ-as: showing that the actual force in the Langevin dynamics can be decomposed as the gradient of a scalar potential, plus an additional force. Note that such an additional term, j ν /P st , can be seen as a curl flux force, as it is necessarily divergence free in the steady state (∂ µ j µ = 0); thus, it has neither sources nor sinks * . Figure 5. Non equilibrium potentials for the linearized dynamics, as calculated following [68]. White curves are some of the points calculated using V (q) = t * 0 pqdt along classical (or Hamiltonian) paths starting at the deterministic fixed-point, the surface represented in color code is a fit. The inset is the same plot, projected in the plane V = 0. Left panel: for the case ω f f = 0, the linearized dynamics is normal and the potential can be calculated analytically V = λ 1 x 2 /2 + λ 2 y 2 /2 and F = −σB∇V /2 (the force is conservative; i.e. it derives from the scalar potential V ). Right panel w = 10, the potential becomes much narrower in the y-direction. λ 1 = 0.1, λ 2 = 0.11. The projection in the inset shows that most-probable paths are much more curved and anisotropic than in the other case. In this case, the system dynamics do not simply follows the gradient of the scalar potential, but it has a curl component.
In order to determine explicitly the two components of the potential for Eq.(6), we first evaluate the scalar potential, V , by making use of standard methods of the nonequilibrium potentials literature relying on path integrals in the weak-noise limit [67][68][69][70][71] and subsequently we evaluate the curl flux force through Eq. (12).
In particular, plugging Eq.11 into the Fokker-Planck equation, the problem of finding the stationary solution up to leading order in σ implies solving where p µ ≡ ∂ µ V (14) * Observe that this decomposition into a component which derives from a scalar potential plus a curl flux vector (which eventually can be written as the rotor of a vector potential) is nothing but the usual Helmholtz decomposition of the deterministic force into a curl-free and divergence-free components [72].
which is formally identical to a Hamilton-Jacobi equation in classical mechanics, being H(q, p) a Hamiltonian. From a path integral viewpoint (see e.g. [68,69]), the probability to reach a given point from the deterministic attractor can be expressed as P (q) ∝ exp(−S min (q)), where S min (q) is the minimum action along a "classical" trajectory following the Hamilton's equations of motion,q µ = ∂H/∂p µ , such that starting at the deterministic fixed point q min (and with a value of p as close as possible to 0 ) reaches an arbitrary point q at time t * . Thus, the associated potential can be expressed (integrating Eq.14 along such a path) as and one just needs to compute this integral. In Fig.5 we plot the potential that we have obtained by numerically computing the above action for many classical paths starting at the deterministic fixed point with different values of p (with small modulus) pointing in different directions. Such different paths reach different points q † † spanning the whole space of q values, allowing us to reconstruct V (q) [68].
On the left of Fig.5 we represent the case with ω f f = 0, where the force is conservative, the resulting potential has a standard paraboloid shape centered at the origin. On the other hand, in the case ω f f = 0 (in particular, ω f f = 10), the classical trajectories are much more curved, and the resulting potential becomes asymmetric, and much narrower along the direction transversal to the diagonal, as illustrated in the right panel of Fig. 5. This reveals that fluctuations along the diagonal are strongly more likely than perpendicular ones.
Once the scalar potential has been determined, we compute the curl force, j µ /P st from Eq. (12). The result is illustrated in Fig.6: in the ω f f = 0 case, the system is conservative, there is no current nor a curl in the stationary state. In this case, the most probable path from any given point to the minimum of the potential is given by the steepest gradient descent (i.e. the deterministic path) and, vice versa, the most probable path leading to a fluctuation from the minimum to such a given point is the time-reversed of the deterministic path (as it is always the case in equilibrium problems [68]).
On the contrary, in the case ω f f = 10, the dominant contribution to the system dynamics is given by the curl flux force (represented by white arrows in Fig.5). In this case the paths clearly deviate from the steepest descent of the scalar potential and moreover they are no longer reversible, as clearly shown in Fig.6. Indeed, the figure illustrates that the most likely path to go from one point to another is quite different from the most likely path to complete the reverse trip. This reveals the breaking of the Taking strictly p = 0 trajectories do not leave the deterministic attractor, thus, vectors p with very small but not vanishing modulus -and varying directions-need to be considered. † † A linear interpolation algorithm has been employed to obtain a continuous curve, for visualization purposes detailed balance condition, illustrating the existence of internal probability currents, i.e. the non-equilibrium nature of the problem [73,74]. Figure 6. Representation of the scalar potential (color coded) and the curl flux force (vector field). Left panel: For the particular (normal) case w = 0, the residual curl flux force vanishes, since the force and the dynamics is fully specified as the gradient of the scalar potential. On the other hand (Right panel) in the (non-normal) case w = 10, the system is constrained close to the manifold y = 0 by the potential (the only fixed point lies at (x, y) = (0, 0) and is weakly stable) and is subject to a curl flux determined by the non-conservative part of the dynamics ((W x , W y ) (wy, 0) for σ → 0). Observe, in particular, that the most-probable fluctuating paths to go from the point marked with an asterisk to the one marked with a plus sign, is different from the one to make the reverse jump. This reveals the intrinsic non-equilibrium (i.e. broken detail balance) nature of the problem, while in the case w = 0 (Left) one trajectory is the reverse of the other, i.e. detailed balance is obeyed.

Stochastic description of the full dynamics
The mechanism that generates anomalously large fluctuations around the up-state fixed point has been already rationalized in the previous sections. However, in spite of this, there is an important aspect of the actual full dynamics that remains to be clarified. This is the fact that the full system becomes trapped into a state of extremely low activity (a down state) for relatively large times (see Fig.2). This phenomenon, which is essential for the system to exhibit avalanching behavior, is not predicted by the deterministic dynamics nor by the above weak-noise approach to the linearized stochastic dynamics. Thus, to shed further light onto it, one needs to go beyond such approximations, hopefully, constructing an effective potential for the full stochastic dynamics.
In principle, the effective potential for the full problem might be calculated using the same method as described above for the linearized dynamics [67][68][69][70][71]. However, for the full problem, owing to the presence of multiplicative/demographic noise, it is necessary to go beyond the Hamilton-Jacobi approach and consider the next-toleading order approximation [67]. We tried to perform such a calculation, but found a number of mathematical pitfalls; these technical difficulties make the problem highly non-straightforward and beyond the scope of the present paper (we leave a detailed mathematical analysis of this for a future work).
For this reason, we resorted to measure the steady state probability distribution P st in a purely computational way, and use it to determine the stationary scalar potential V (Σ, ∆) = − log P st (Σ, ∆). Fig.7 shows that such an effective potential has a deep minimum close to the origin, as the potentials of stochastic systems with multiplicative/demographic noise generally do [75] (see also [76]). Observe that the minimum is not a singularity just because h = 0, preventing a true absorbing state to exist. Observe the presence of a demographic-noise-induced minimum (absent in the previous linear-approximation approaches) very close to the origin. Interestingly, observe that although the potential tends to keep the system dynamical state confined very close to the diagonal there exists a strong curl flux force (as shown in the projection on the E-I plane) that pulls the system into the (multiplicative)noise-induced minimum. Reddish colors represent the different minima. Observe that we have also plotted the projection of the potential into the E variable, which is equivalent to the potential for the I variable. Parameters: h = 10 −6 , ω E = 7, ω I = 34/5, α = 0.1.
Moreover, Fig.7 shows that, on the one hand, the potential keeps the dynamics confined very close to the diagonal (E = I) and, on the other hand, as soon as the system moves slightly away from it, the dynamics is governed by the strong curl flux force -evaluated numerically form Eq.(12)-, which, either pulls the system directly towards the noise-induced minimum at the origin, or pulls it towards the deterministic minimum.
This allows us to rationalize the emerging avalanches: the system is trapped in the noise-induced minimum; small fluctuations drive it slightly away from the minimum, and then, the curl flux force drives the system towards the deterministically stable up state. However, the system does not dwell long around such an up state: given that the deterministic dynamics is non-normal and possibly reactive, and that any tiny stochastic fluctuations drive the system away from the fixed point while the curl flux force drags it back to the down state, closing the cycle of an individual avalanche.
Summing up, the emerging curl flux keeps the system jumping back and forth (circulating) from the down state to states of larger activity. This is the mechanism for the emergence of avalanching behavior in deterministically stable systems, once finitesize fluctuations are considered.
Finally, let us mention that in Appendix B, we analyze how all this phenomenology depends on system size (i.e. noise amplitude) and, moreover, through an effective onedimensional approach we explain why the measured avalanche exponents are compatible with those of a standard random walk.

Avalanching without reactivity
It should be noted that the linear stability analysis of the system at the deterministic fixed point (in particular, whether the linearized dynamics is non-normal and reactive, or not) gives a very good idea of the possibility of having the type of effects discussed in the main text, e.g. large transient departures from the fixed point in the presence of fluctuations.
However, it is important to keep in mind that the linearized dynamics does not include all the information of non-linear systems, and one could expect some deviations from linear-approximation-based knowledge, at least in some cases.
Thus we asked the question whether it could be possible to have avalanching behavior, as described in the main text, in the case in which the up-state fixed point is stable, weakly non-normal and non reactive at all. Surprisingly, the answer is yes! In order to illustrate this effect, we consider the set of stochastic Eqs.(3) with a particular choice of the parameter values (see caption of Fig.8) for which the linearized dynamics around the fixed point is strongly stable, quasi-normal and non-reactive but still -very importantly-the two nullclines are very close to each other.
As illustrated in Fig.8 (left) -even if the fixed point is strongly stable-the deterministic field of forces still has a scar near the diagonal, shear flows exist at both sides of it, and one can anticipate large excursions away from the up state once fluctuations are taken into account. As a matter of fact, Fig.8 (right) illustrates that noise amplitudes above some threshold (larger than in the cases discussed in the main text) lead the system to exhibit large excursions alternating up and down states, in a rather intermittent -avalanching like-way. Thus, in conclusion, under some circumstances it is possible to have very nearby nullclines, generating a "scar" in the velocity field -opening the door to avalanching phenomena when noise is included-even in the absence of non-normal/reactive deterministic dynamics at the fixed point. The key ingredient is the proximity of the two nullclines, regardless of whether when they intersect they do it in such region of proximity (inducing the two eigenvectors to be almost degenerate, and the matrix non-normal), or they intersect in a region whether they are separated, as in Fig.8 (Left). We propose to call this mechanism, non-linear reactivity. and ω I = 17. The proximity between the nullclines makes it easy that, above certain perturbation, the system exhibits large excursions from the up state to the down state, even though the non-normality is N N = 10 −4 and the reactivity R = −1.08. We have called this mechanism, non-linear reactivity.

Discussion and conclusions
The mechanism of amplification of fluctuations was proposed by Murphy and Miller [41] and later employed by Benayoun et al. [40] to explain the origin of avalanching phenomena in neural systems. It caught our attention as a possible alternative to criticality to justify the emergence of scale-invariance in neuronal networks, as first measured empirically by Beggs and Plenz [14].
In particular, Benayoun et al. observed that the model of Wilson-Cowanvery well-known and accepted as a mesoscopic description of neural dynamicsincluding excitatory and inhibitory neuron populations, is able to describe heavy-tailed distribution of bursts of spontaneous neuronal activity if it is tuned to a condition of an approximate balance between excitatory and inhibitory couplings. Importantly, this occurs without the need for a perfect fine tuning to the very critical point, at which excitation and inhibition exactly compensate to each other.
Such a mechanism of "balanced amplification" stems from the existence of a rather peculiar type of stable fixed point at the deterministic level (i.e. for infinitely large neural populations). In particular, as we have profusely illustrated here, if the Jacobian characterizing the linearized deterministic dynamics is a non-normal matrix -which happens if the conditions for balanced amplification hold-then, the corresponding eigenvectors do not form an orthogonal basis and, as a consequence, some directions can be poorly represented, implying that the components of vectors in such directions can become very large in the eigenvectors basis. Moreover, as discussed along the paper and more specifically in Appendix A, this implies that the system can be reactive, meaning that even if the fixed point is stable, perturbations from it can grow significantly in modulus before decaying back. We have shown that reactivity is a more stringent condition than non-normality: all reactive matrices are non-normal, but the opposite is not true.
This type of anomalous deterministic dynamics induces the possibility of transiently evolving away from the deterministic fixed point. When finite populations are considered, intrinsic or demographic noise needs to be taken into account. Such a noise induces the emergence of a new attractor of the dynamics for very low activities. Given that such an attractor -a down state-is deterministically unstable, fluctuations, can make the system escape from it, experiencing large excursions before falling back again into the down state, thus giving rise to avalanches of activity.
At the light of all this, one might be tempted to conclude that non-normality (and, possibly, reactivity) -together with intrinsic stochasticity-are the key ingredients for the described phenomena to emerge. However, as we have shown, such properties of the linearized dynamics are not actually necessary conditions. The main required ingredient in the deterministic dynamics is that the two nullclines need to be very close in a broad region of the phase portrait -even if they do not intersect in such a region-and this is a non-linear feature. Actually, we propose to call this non-linear reactivity. Still, the easiest way to have very close nullclines is when the two nullclines intersect when they are almost parallel, and in such a case the associated linear stability matrix is most likely non-normal, and possibly reactive.
Thus, the "balanced amplification" condition implies that the system is able to generate some type of "scar" or shear flow in the deterministic velocity field, allowing for very large excursions away from the up state when the system is perturbed. This together with demographic noise -which has two effects: (i) taking the system away from the deterministic fixed point, and (ii) creating a second noise-induced effective attractor of very low activity-explains the emergence of avalanches and large excursions between the down and up states.
In other to further rationalize all this, we have constructed effective potentials to describe the dynamics. We found, first that, in the case of balance amplification, the dynamics cannot be simply derived from the gradient of a scalar potential. There is a very important contribution from a curl flux (or, equivalently, a vector potential), that explains why the system circulates back and forth between states of low an high activity. Importantly, the global potential is a non-equilibrium one, and the most likely trajectory to jump from one state to another does not coincide with the most probable one to do the reverse jump.
Moreover, we have constructed an effective one-dimensional description of the full problem, that allowed us to understand that the avalanches around the down state, especially the small ones, can be well described by the excursions of a random walk around an almost flat potential. The resulting avalanche exponent values coincide with those of the unbiased random walks (α = 3/2 and τ = 4/3) and, thus, are not compatible with the empirically measured ones (in neuronal tissues) which are better fitted by the exponents of an unbiased branching process (α = 2 and τ = 3/2). In any case, these are not perfect power-laws, i.e. the system is not truly scale invariant, as above some level of activity, the system can get trapped for some time around the up state, breaking the freely-moving random walk picture.
On the other hand, preliminary results to be fully developed elsewhere show that the scenario described here -combining non-normal dynamics with stochasticity-can appear in other neural systems (e.g. including synaptic plasticity rather than inhibition as a chief regulatory mechanism), clearly illustrating its generality. Also, the present mechanism is very likely to play a role in population dynamics in ecology, a context in which reactivity has been profusely discussed at a deterministic level (see [59]).
Finally, we hope that this work fosters further analyses of non-standard deterministic dynamics -in the sense of showing either non-normal forms of the stability matrix and/or shear flows-in noisy dynamical systems of relevance in other fields and in the presence of more complex networked architectures.

Acknowledgments
We are grateful to the Spanish-MINECO for financial support (under grants FIS2013-43201-P and FIS2017-84256-P; FEDER funds). MAM also acknowledges the support from TeachinParma and the Cariparma foundation. Last but not least we are happy to acknowledge P. Moretti, P. Garrido, S. Suweis, A. Maritan, and especially V. Buendia for extremely useful comments.
dynamics, which -in some cases-can be affected by long transients. In Fig.A1 we illustrate a (non-normal) stable linear system can transiently evolve away from its stable fixed point or, in other words, a perturbation can transiently grow before decaying to zero [59]. Actually, the temporal evolution of a perturbation x(t) is described by: where H(J) ≡ (J + J T )/2 is called the Hermitian part (or symmetric part) of J.
The maximum amplification rate immediately following the perturbation, specifying the reactivity (or Rayleigh quotient of H) of the system can thus be written as: which is nothing but the maximum eigenvalue of H(J) [78]. This means that the strongest magnification that can occur must take place in the eigenspace belonging to the eigenvalue with the largest modulus -for all other vectors the magnification will be smaller-and the magnification is equal to the maximum eigenvalue of the Hermitian part itself. Thus, in summary, while the asymptotic behavior is determined by the maximum eigenvalue of matrix J, the transient behavior is determined by the maximum eigenvalue of the Hermitian part H(J).
If J is not symmetric, then J = H(J) and, thus, it is possible that the maximum eigenvalue of J is negative, while the maximum eigenvalue of H(J), λ H max , is positive, meaning that the equilibrium specified by J is stable but reactive [59].
Observe, in particular, that Rayleigh's principle does not hold for non-normal matrices † , i.e. if J is non-normal, then it is not true that the strongest amplification must occur in the direction of the eigenvector corresponding to the largest eigenvalue, thus, there can exist strongly-amplified directions. In particular, the Schur theorem [52,53] states that non-normal matrices can always be transformed in upper triangular matrices, having the eigenvalues on the diagonal and a "feedforward" dynamics outside the diagonal, which cannot be neglected and whose relative weight determines the strength of the non-normality.
Finally, let us remark that reactivity is a more stringent condition than nonnormality: reactive matrices are always non-normal, but the opposite is not true. Let us clarify the interplay between reactivity and non-normality through a simple example.
The condition b 2 > 4d specifies reactivity whereas b > 0 is the condition for non-normality. Thus reactivity is a more stringent condition than non-normality. † Note that H(J) is Hermitian and all Hermitian matrices are normal. Appendix B. Effective one-variable description: stochastic tunneling and the origin of the avalanche exponents Despite the fact that the main ingredients of the avalanching phenomenon have been carefully scrutinized in previous sections, we still do not have a transparent explanation for the computational observed avalanche exponent values (compatible with the statistics of standard unbiased random walks), nor a detailed analysis on how the phenomenon depends on noise amplitude (i.e. on system size).
Aimed at shedding further light on these issues, in this appendix we present a one-dimensional effective description to the full two-dimensional dynamics Eq.(3).
We have already computed the bivariate probability distribution P (Σ, ∆) = exp(−V (Σ, ∆)). Marginalizing over ∆, one obtains P (Σ) = ∞ 0 P (Σ, ∆)d∆, from which one can define an effective one-dimensional potential, V ef f (Σ) = − ln P (Σ), as illustrated in Fig.B1A for different values of the noise amplitude. Observe, in particular, that the effective potential becomes more and more peaked around the deterministic solution (Σ 0.5) as the noise amplitude is reduced (i.e. for large systems sizes), while for bigger noise amplitude a new noise-induced minimum emerges near the origin. This Figure B1. (a) Effective potential for the overall activity Σ for different values of the noise amplitude (σ, see legend). Its shape changes from a single potential well for low noise (dark red curve) to a bistable situation alternating between up and down states and, finally, to a new single potential peaked around the down state for larger noise (blue curve). explains the unusually large permanence times of the original two-variable system into low activity regimes, which is ultimately responsible for the avalanching behavior. Moreover, in Fig.B1B -where a double logarithmic scale is employed-the structure of the potential V ef f (Σ) can be further inspected, revealing the presence of a logarithmic decay. Also, in Fig.B1B we plot -with dashed lines-the complementary potential V ef f (∆) obtained from P (∆) = ∞ 0 P (Σ, ∆)dΣ, throughṼ ef f (∆) = − log(P (∆)). Observe that for all noise intensities the potential has a minimum around ∆ = 0, indicating that the dynamics spends most of the time close to the diagonal E = I.
To have some additional understanding of the origin of the noise-induced minimum of V ef f (Σ) near the origin one can imagine (just for argument's sake) that the dynamics of the system was strictly constrained to evolve along the E = I diagonal, and let us define a density variable ρ = E = I), which implies fixing ∆ = 0. Performing a Taylor expansion and keeping only leading terms in Eq.(3) gives rise tȯ ρ = h + aρ − bρ 2 +σ (ρ + h )η, (B.1) with a = (−α − h + ω 0 ), b = ω 0 ,σ = σ (ω 0 + α − h), h = h/(ω 0 + α − h), and η is, as above, a Gaussian white variable. The effective potential resulting from this simple one-dimensional dynamics can be easily calculated (by writing the associated Fokker-Planck equation and computing its steady state) A sketch of this potential is shown in Fig.B1C. Observe that it reproduces all the key features of the computationally determined one in Fig.B1B. Finally, let us remark that this approximation is not exact: actually, it is not true in the actual dynamics that the system remains confined to the diagonal (i.e. E = I, and thus ∆ = 0); indeed, one of the key ingredients of the complete dynamics is the shear flow occurring out of the diagonal, i.e. the non-conservative contribution from the curl flux. Thus, in order to mimic the effect of off-diagonal shear flow-we additionally introduce a "stochastic tunneling effect" in our one-variable (ρ) effective description, such that the system can spontaneously jump from the down state to the up one and the other way around with some fixed probabilities (p). In this way, we have constructed a new"effective" one-dimensional dynamics for the full (two-dimensional) problem. While the resulting potential keeps essentially the same shape as shown in Fig.B1B, we have computed the avalanche distribution using this effective dynamics; results are illustrated in Fig.B2. In particular, the distributions are consistent -at least for small avalanches-with the first-return statistics of unbiased random walkers, justifying the computationally measured exponents in the main text. For larger avalanches, there is always a bump in the probability distributions -much as in the neural modelcorresponding to avalanches remaining for a long time near the up fixed point. Figure B2. Avalanche size distribution P (S) (main) and duration distribution P (T ) (inset) for different values of p in the effective one-dimensional model. Observe that the exponent values are compatible with those of and unbiased random walk for small avalanche sizes. For larger sizes, deviations from perfect scaling occur (as they do in the original two-dimensional system), in particular "bumps" apper. These reveal the existence of a more complicated dynamics with jumps to the up state, where the system can get trapped for some time, before the avalanche ends. Parameter values: a = 0.5, b = 1, σ = 10 −3 , h = 10 −3 , T = 2.5 · 10 −4 .