Unconventional criticality, scaling breakdown, and diverse universality classes in the Wilson-Cowan model of neural dynamics

PHYSICAL REVIEW E 108, 034110 (2023) Unconventional criticality, scaling breakdown, and diverse universality classes in the Wilson-Cowan model of neural dynamics Helena Christina Piuvezam ,1,*Bóris Marin ,2Mauro Copelli ,1and Miguel A. Muñoz 3,† 1Departamento de Física, Universidade Federal de Pernambuco, Recife PE 50670-901, Brazil 2Centro de Matemática, Computação e Cognição, Universidade Federal do ABC, São Bernardo do Campo, Brazil 3Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Granada, Spain (Received 18 January 2023; accepted 4 August 2023; published 11 September 2023) The Wilson-Cowan model constitutes a paradigmatic approach to understanding the collective dynamics of networks of excitatory and inhibitory units. It has been profusely used in the literature to analyze the possible phases of neural networks at a mean-field level, e.g., assuming large fully connected networks. Moreover, its stochastic counterpart allows one to study fluctuation-induced phenomena, such as avalanches. Here we revisit the stochastic Wilson-Cowan model paying special attention to the possible phase transitions between quiescent and active phases. We unveil eight possible types of such transitions, including continuous ones with scaling behavior belonging to known universality classes—such as directed percolation and tricritical directed percolation—as well as six distinct ones. In particular, we show that under some special circumstances, at a so-called “Hopf tricritical directed percolation” transition, rather unconventional behavior is observed, including the emergence of scaling breakdown. Other transitions are discontinuous and show different types of anomalies in scaling and/or exhibit mixed features of continuous and discontinuous transitions. These results broaden our knowledge of the possible types of critical behavior in networks of excitatory and inhibitory units and are, thus, of relevance to understanding avalanche dynamics in actual neuronal recordings. From a more general perspective, these results help extend the theory of nonequilibrium phase transitions into quiescent or absorbing states. DOI: 10.1103/PhysRevE.108.034110 I. INTRODUCTION A large variety of natural systems exhibit a continuous (second-order) phase transition between an active phase and a quiescent (or absorbing) one where all activity ceases [1–5]. This phase transition is often associated with scaling behavior and is typically described by the directed percolation universality class, as originally conjectured by Janssen and Grassberger [6,7]. Actually, directed percolation (DP) is one of the most robust classes of universal critical behavior away from thermal equilibrium [1–5,8], as it describes all possible phase transitions into an absorbing state in the absence of additional symmetries or conservation laws, even for multicomponent systems [1,3–5,9,10]. Moreover, some of the representative models of this class, such as the branching process and the contact process [11,12], have been broadly studied in a large variety of contexts, including countless applications in materials science, turbulence, epidemics, theoretical ecology, social sciences, and neuroscience. Conversely, under some circumstances, phase transitions into quiescent states occur in a discontinuous (or first-order) rather than continuous manner. This is often the case when higher-order reactions are considered, e.g., where more than one active unit is required to activate the third one [5,13,14]. This situation usually involves a bistable regime with phase coexistence and hysteresis. There are also well-studied systems that include both types of transitions, continuous and *[email protected] †[email protected] discontinuous, as well as a tricritical point with a scaling behavior that differs from DP and is described by the so-called tricritical directed percolation (TDP) universality class (see e.g., a modified contact process [13,15]). In the context of neuronal systems, the experimental work by Beggs and Plenz reported on the existence of neuronal avalanches (i.e., outbursts of neuronal activity between quiescent periods). These exhibited highly variable sizes and durations, which were both distributed as power laws. Moreover, the associated exponents were found to be consistent with those of critical systems in the mean-field DP universality class [16], suggesting that brain dynamics could be poised near the edge of a phase transition [17–22]. Further experimental works reported evidence of scaling exponents deviating from DP [23,24], so that the interpretation of the scaling behavior in terms of universality classes remains a current matter of debate [19]. In particular, the possible departure from the standard DP class (together with the possible existence of discontinuous transitions in brain dynamics [25–27]) raises a number of questions from a theoretical point of view. For example, the fact that neuronal networks include inhibitory units—which hinder activity propagation and are not usually included in simple models in the DP class, such as the standard branching process—triggered a renewed interest in the scrutiny of novel types of phase transitions in networks of excitatory and inhibitory units [28–36]. Do different types of quiescent to active phase transitions emerge in simple models of activity propagation once inhibitory effects are considered? Here, to further advance our knowledge along these lines, we revisit one of the most broadly studied parsimonious models in neuroscience: the Wilson-Cowan model [37]as 2470-0045/2023/108(3)/034110(18) 034110-1 ©2023 American Physical Society PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) well as its stochastic counterpart [28,36]. We systematically analyze the resulting phase diagram, carefully analyzing all the possible phases and phase transitions. In particular, we reveal that, depending on the relative strengths of excitatory and inhibitory couplings, there can be up to eight different types of phase transitions into quiescence. Some of them exhibit well-known scaling behavior (such as DP or TDP), while others are discontinuous or show different types of anomalies in scaling or even mixed features of continuous and discontinuous transitions. Finally, we elucidate a unique type of phase transition that is highly nontrivial, exhibiting unconventional behavior and breakdown of scaling. Our results help rationalize and categorize the possible types of criticality in networks of excitatory and inhibitory units, contributing to the advance of the brain-criticality hypothesis and of the general theory of nonequilibrium phase transitions [4]. II. THE WILSON-COWAN MODEL AND ITS STOCHASTIC COUNTERPART In its original formulation, the Wilson-Cowan model describes the collective deterministic (or “mean-field”) behavior of a local population of both excitatory and inhibitory neurons by means of two coupled differential equations [37]. Such equations reproduce—as a function of a set of couplingstrength parameters—a variety of possible dynamical regimes (all of which have counterparts in actual neuronal systems [33,38–40]) that are delimited by phase transitions (bifurcation lines) [41,42]. To go beyond this deterministic or mean-field picture, Benayoun et al. [28] proposed a microscopic version of the Wilson-Cowan model in the form of a Markovian process for a population of coupled excitatory and inhibitory individual binary neurons that can be either active or inactive [43]. In the so-called stochastic Wilson-Cowan (SWC) model, the state of each unit at a given time t—which can be either excitatory (E) or inhibitory (I)—is given by σE/I (t)= 1 for active neurons and σE/I (t)=0 for inactive ones. These discrete state variables change according to a master equation specified by transition rates, which are defined as follows. Each active neuron, regardless of its type, shifts from the active to the quiescent state, 1 →0, at a constant decay rate α. The reverse transition (activation), 0 →1, occurs at rate (s), defined as (s)=tanh(s),if s>0 0,otherwise,(1) where the input sto neuron, is s= m wmσm+h,(2) wmis the synaptic weight from neuron mto neuron , and his a constant external input. Observe that the form of the response function,(s), in Eq. (1) enforces the nonnega- FIG. 1. (a) Sketch of the Wilson-Cowan model, including both excitatory and inhibitory populations and synaptic couplings between them. (b) Owing to the piecewise definition of the response function in Eq. (1), there are three different regions in the state space depending on the total input son excitatory and inhibitory populations, respectively. In particular, in region I, sEand sIare both positive; in region II, sEis negative and sIis positive; and, finally, in region III, sEand sIare both negative (the figure illustrates these three regions for wEE =2.5, wEI =1.5, wIE =1.5, wII =0.5, and h=0). Observe that when wEE/wEI >wIE/wII, the lines sE=0 and sI=0 switch position, and region II then shows sE>0and sI<0 (a condition that is not explicitly explored here). The trajectories (red-dashed arrows) show how the system is attracted either to the origin (quiescent state) or to a nontrivial active state or fixed point (FP). As discussed in Sec. IV, initial conditions in region I that are close to the switching manifold sE=0 can cross over to region II. In this case, regions II and III are trapping ones: once trajectories cross the switching manifold, they are unable to return to region I (see Appendix Cfor a more detailed explanation). tivity of the transition rates. In what follows, the synaptic weights are chosen to depend only on the type (excitatory or inhibitory) of both the presynaptic and the postsynaptic neuron, leaving [as sketched in Fig. 1(a)] only four free parameters: wm∈{wEE,wIE,wIE,wII}∀, m, where, e.g., wIE is the excitatory coupling strength to inhibitory neurons, and so forth. Previous works on this model have often employed symmetric weights—common excitatory (wE≡wEE =wIE) and inhibitory (wI≡wII =wEI) inputs [28,33,44]— as a way to reduce the dimensionality of the phase diagram. In order to systematically explore the full set of possible phase transitions, here we do not impose such constraints. This stochastic process can be implemented on different types of networks, as specified by a connectivity matrix. As a first approach, one can assume a large fully connected network of size N. Indeed, performing a standard size expansion [45,46], one recovers (up to leading-order) the well-known Wilson-Cowan equations [28] (see Fig. 1) 034110-2 UNCONVENTIONAL CRITICALITY, SCALING … PHYSICAL REVIEW E 108, 034110 (2023) written as ˙ E=−αE+(1−E)(wEEE−wEII+h),(3) ˙ I=−αI+(1−I)(wIEE−wII I+h),(4) where Eand Iare the densities of active excitatory and inhibitory neurons, respectively [28] and the coupling constants are defined including a rescaling with the network size. Similarly, by adding next-to-leading corrections, one obtains a set of two Langevin equations which we do not write explicitly here (we just remark that they include square-root noise following the lead of DP and TDP classes [1,3–5]). The pair of coupled stochastic equations allows one to describe fluctuation effects in finite-size (fully connected) networks [28,45,47]. Nevertheless, let us stress that the forthcoming computational analyses refer to simulations of the full microscopic model. Observe that, owing to the discontinuous derivative of  at zero, Eq. (3) and Eq. (4) are piecewise smooth differential equations [48–50]—i.e., they are smooth everywhere except at switching manifolds, which are defined by the conditions of vanishing input in the response function in the previous equations: sE≡wEEE−wEII+h=0 and sI≡wIEE−wII I+ h=0. These two conditions divide the state space into three regions: I, II, and III, as illustrated in Fig. 1(b) (for h=0): (1) In region II, the input to excitatory populations is negative, so that Eq. (3) becomes ˙ E=−αEand trajectories in this region decay exponentially fast to the quiescent phase (either crossing to region III or not; see Appendix Cfor a more detailed explanation). (2) Similarly, in region III, the input to both excitatory and inhibitory populations is negative, so that ˙ E=−αEand ˙ I= −αI, leading to an even faster decay to quiescence. (3) Conversely, in region I, the total input does not vanish for either subpopulation and the dynamics can be more complex, possibly reaching nontrivial (active) fixed points. Inspection of Eq. (3) and Eq. (4) readily reveals that trajectories starting in regions II and III do not cross over to region I (as excitation always diminishes in these regimes), but the opposite can happen (see, e.g., the central trajectory shown in Fig. 1as well as Appendix C). In the next sections, we explore in detail, both analytically and numerically, the features of each of the possible phase diagrams as well as all the possible phase transitions between quiescent and active states. III. MEAN-FIELD PHASE DIAGRAMS: GENERAL AND SPECIFIC FEATURES To avoid confusion, let us first underline that in what follows we refer indistinctly to phase transitions or to bifurcations, as the present focus is on the description of fully connected networks (i.e., mean-field systems). Thereby, DP transitions correspond to transcritical bifurcations, discontinuous transitions to saddle-node bifurcations, tricritical points to saddle-node-transcritical (codimension-2) bifurcations [51,52], and so forth. In the absence of any external driving force (h=0), the steady-state conditions for Eq. (3) and Eq. (4) always admit a trivial solution E∗=I∗=0, which defines the quiescent phase as well as, possibly, some nontrivial solutions (E∗>0 and I∗>0) of the following equations: E∗=1 wIE wIII∗+−1αI∗ 1−I∗,(5) I∗=1 wEI wEEE∗−−1αE∗ 1−E∗,(6) which define the active phase. Observe that Eq. (5) and Eq. (6) are well-defined only as long as −1exists, i.e., in region I [Fig. 1(b)], so that nontrivial solutions exist only inside said region. In what follows, we analyze the overall phase diagram, describing the stable phases as a function of the model parameters. In particular, without loss of generality, we keep the activity-decay rate α= 0 and the self-inhibition weight wII ⩾0 fixed. Choosing wEE and wEI as control parameters, the system may display three qualitatively different types of phase diagrams in the (wEE,wEI) plane depending on the value of the remaining free parameter, wIE, as illustrated in Fig. 3. Other parameter choices are possible, but the system is always described by one of these three qualitatively different types of phase diagrams. A. Quiescent phase and its stability limits First of all, let us stress that the quiescent state is always stable (and is the only stable state) with respect to the introduction of inhibition-dominated perturbations (i.e., for initial conditions in regions II and III); in other words, activation of inhibitory neurons does not allow the system to escape from the quiescent state. Therefore, in what follows, we focus on its stability and the resulting phase diagram as a result of excitation-dominated perturbations. Importantly, there are two different types of quiescent phases: The first one is a standard quiescent phase, i.e., a regime in which the quiescent phase is locally stable to excitation-dominated perturbations [Fig. 2(a)]. This occurs if the eigenvalues of the associated stability matrix, as specified by λ±=wEE −2α−wII ±(wEE +wII)2−4wEIwIE 2,(7) have negative real parts (white zone in the diagrams of Fig. 3). Alternatively, if the eigenvalues have positive real parts and an imaginary component, then, in principle, one could expect oscillations away from quiescence to emerge. However, given the nonsmooth piecewise dynamics, the resulting “curvy” trajectories end up crossing over to region II, where the dynamics follow the equation ˙ E=−αEand the quiescent state is the only attractor. Therefore, in this regime, small excitatory perturbations to the quiescent state may give rise to large trajectories in state space before returning to quiescence (see Fig. 2(b) and [28,36]). This property is called “excitability” [53] (or “reactivity” [54]) and the corresponding quiescent state is called “excitable quiescent.” Observe that, inside the excitable quiescent phase, there is a region where the eigenvalues are both real and positive. The trajectories continue to fall into regions II and/or III, therefore, the phase remains the same. When in the standard quiescent phase, the 034110-3 PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) FIG. 2. Two different types of quiescent phases. Time series towards the absorbing state of trajectories (a) on the standard quiescent phase (wEE =1.2andwIE =0.2) and (b) on the excitable quiescent phase (wEE =2.2andwIE =1.0). Observe the nonmonotonicity in the second case, which is a manifestation of the excitability of the quiescent state: perturbations can be amplified before trajectories finally decay to quiescence. The insets show the phase space for these two cases, respectively, as well as the corresponding switching manifolds and some sample trajectories (arrows). In the second case, trajectories cross the switching manifold. Parameter values are α=1.0, wII =0.2, and wEI =2.0. quiescent state loses its local stability when the real part of the largest eigenvalue becomes positive. First, considering a real eigenvalue, it occurs at wT EE =α+wEIwIE α+wII (8) as can be easily seen from Eq. (7); note that this condition corresponds to the diagonal blue dashed line in all the plots in Fig. 3, representing a line of transcritical bifurcations (as well as its continuation to the right). Not surprisingly, separating the previous two phases (standard quiescent and excitable quiescent), there is a line of (supercritical) Hopf bifurcations (dot-dashed vertical lines in Fig. 3) occurring at wH EE =2α+wII,(9) which exist only when there is a nonvanishing imaginary part, i.e., from Eq. (7), at (wEE +wII)2−4wEIwIE <0,(10) so that the Hopf bifurcation is defined only above the Hopftranscritical point, i.e., above the diagonal blue line in all the phase diagrams in Fig. 3. Summing up, there are two types of quiescent phases, separated by a line of Hopf bifurcations: (1) A standard quiescent phase, which is represented by a locally stable fixed point (upper left regions in Fig. 3) (2) An excitable quiescent phase, in which local stability is lost, but the origin remains globally stable because the trajectories cross to trapping regions (upper right regions in Fig. 3). Finally, let us note that the excitable-quiescent state may lose its global stability (in favor of a nontrivial active state, as described below) at a transition line that needs to be numerically determined (black dashed line in Fig. 3), beyond which trajectories flow to the alternative active state. FIG. 3. Different phase diagrams of the model [Eq. (3)andEq.(4)], with excitation-dominated initial conditions and parameter values: α=1, wII =0, and wIE =3 (Case A), wIE =1 (Case B), and wIE =0.8 (Case C). In all cases, the (vertical dot-dashed) line of Hopf bifurcations separates the standard quiescent phase from the excitable quiescent one. (a) In case A, the Hopf line collides with the (blue dashed) transcritical line to the right of the tricritical point (black dot), within a region of bistability. (c) In case C, on the contrary, the intersection of the Hopf line with the diagonal line occurs to the left of the tricritical point. (b) Between the previous two cases, case B (for which fine-tuning a third parameter, wIE =1, is needed) the Hopf line collides with the transcritical one at a codimension 3 bifurcation point that we call Hopf-tricritical point [Hopf saddle-node-transcritical bifurcation that corresponds to, as we called it, the Hopf tricritical directed percolation (H+TDP) point]. In all cases, the standard quiescent phase loses its stability at a line of transcritical bifurcations (blue dashed lines), while the excitable phase loses its global stability from below (i.e., as wEI decreases) at the black dashed line—that may be very close to (or coincide with) the (orange) line of saddle-node bifurcations. Finally, the horizontal black segments T1,...,T8represent eight qualitatively different ways to transition from a quiescent to an active state as the control parameter, wEE, increases. 034110-4 UNCONVENTIONAL CRITICALITY, SCALING … PHYSICAL REVIEW E 108, 034110 (2023) B. Active phase and its stability limits The active state or phase (with E∗= 0 and I∗= 0) becomes a stable solution either at (1) a transcritical bifurcation (i.e., it emerges continuously once the quiescent phase loses its stability in a DP transition), which occurs at Eq. (8)as represented by the blue dashed lines in Fig. 3; (2) a saddlenode (SN) bifurcation, i.e, emerging discontinuously (orange solid line in Fig. 3)at wSN EE =min E∗wEII∗(E∗) E∗+1 E∗−1αE∗ 1−E∗,(11) where E∗and I∗are solutions of Eq. (5) and Eq. (6) that can be solved numerically; or (3) at a tricritical point, where (if) the previous two lines meet (black dot in Fig. 3), to which one can also refer as “saddle-node-transcritical” (SNT) point [its location (wSNT EE ,wSNT EI ) in the phase diagram is explicitly derived in Appendix B; see, in particular, Eq. (B3) and Eq. (B4)]. Let us emphasize that, even if for the naked eye, the saddlenode (orange) curve in all the plots in Fig. 3seemstobea continuation of the transcritical (blue dashed) line, parallel to it, that is not exactly the case; the mathematical conditions for both are different. C. Relative location of the line of Hopf bifurcations Observe that the line of Hopf bifurcations—which as shown in Fig. 3(a) is always vertical in the (wEE,wEI) plane—collides with the line of transcritical bifurcations at a special point [here named, in general, Hopf-transcritical (HT) bifurcation, which is marked with an empty circle in the different panels of Fig. 3]. From Eq. (8) and Eq. (9), one can easily derive the conditions for the HT point: wHT EI =(α+wII)2 wIE ,(12) wHT EE =2α+wII.(13) The key aspect distinguishing the three possible topologies of the phase diagram (A/B/CinFig.3) is whether this HT point lies to the right (case A), left (case C), or on top of the tricritical (SNT) point (case B) in phase space: (1) Case A: wSNT EE <wHT EE , (2) Case B: wSNT EE =wHT EE , (3) Case C: wSNT EE >wHT EE . As already mentioned, these three possibilities are illustrated in Fig. 3, in which the value of wIE is changed (from 3, to 1, and 0.8) to switch between regimes. Also, note that case B requires a higher level of fine-tuning than the other two cases, which appear in broad regions of parameter space. From here on, one needs to separately discuss the three aforementioned possible structures of the phase diagram and carefully classify the diverse types of phase transition emerging in each of them. 1. Case A: Left panel in Fig. 3 In this case, the HT point lies to the right of the tricritical point. Visual inspection of Fig. 3(a) reveals that there are four different ways to go from a quiescent phase (either standard or excitable) to the active one. These are labeled as: T1,for the transcritical bifurcation from the standard quiescent; T2, for a standard tricritical transition; T3, for a transition through a bistable regime (saddle-node bifurcation with coexistence between the standard quiescent and the active phase); and T4, also for a discontinuous transition with bistability, although in this case, between the excitable quiescent phase and the active one. 2. Case B: Central panel in Fig. 3 Here the HT point lies exactly on top of the tricritical point. This structure leads to only three possible types of transitions: T1and T4(as described above), and a transition labeled T5, which occurs through the tricritical (SNT) point that coincides with the special HT point in a codimension 3 bifurcation corresponding to a transition that we call Hopf tricritical directed percolation (H+TDP) transition. 3. Case C: Right panel in Fig. 3 In case C, the HT point lies to the left of the tricritical point and there are five types of transitions, including the standard transcritical (T1) and saddle-node (T4) ones, as well as three unique ones: T6, a transition through the special HT point; T7, a transcritical bifurcation but into the excitable quiescent phase; and, finally, T8, a tricritical (or SNT) transition into the excitable quiescent phase. In the next section, we analyze these eight types of phase transitions (or bifurcations)—from T1to T8—scrutinizing the corresponding peculiarities for each of them. IV. SCALING PROPERTIES AT THE DIFFERENT TYPES OF TRANSITIONS Standard linear stability analysis of the fixed points of the (mean-field) dynamics [Eq. (3) and Eq. (4)] allows one to study the nature of bifurcations and make analytical predictions for the scaling behavior [1,4,55]. In particular, a linear approximation of Eq. (5) around the quiescent solution yields a value of I∗proportional to the density of active excitatory neurons E∗, hence in what follows we employ indistinctly either the latter or the sum of both as an order parameter. In all cases, and for all possible types of transitions, we compute the usual quantities and scaling exponents (as long as they are well defined) customarily employed in the analysis of quiescent-active phase transitions. Even if, generally, three independent exponent values suffice to fully determine the universality class [55,56], here, for the sake of completeness, we determine a larger number of them, which also allows us to check for consistency. In particular, we compute the following: “Static exponents” such as: (1) β, the control parameter one (E∗∝β), where =wEE −w∗ EE is the distance to the transition point and w∗ EE stands, generically, for the value of wEE at the specified bifurcation and (2) δh, defined by E∗∝ h1/δh, representing the response to a constant external field h at criticality. “Correlation exponents”(ν) such as: the one for (3) the correlation length, ξ⊥,ξ⊥∝ν⊥, and for (4) the correlation time, ξ,ξ∝ν. 034110-5 PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) “Dynamic exponents” such as: (5) θ, that governs the time decay of the order parameter E(t)∝t−θ. “Spreading exponents” such as those describing: (6) the total number of active sites, N(t)∝tη; (7) the mean-squared radius in surviving runs R2(t)∝tz; and (8) the survival probability Ps(t)∝t−δ[1]. “Avalanche exponents” defined by: (9) P(S)∼S−τ,forthe distribution of avalanche sizes, S; (10) P(T)∼T−τt, for durations, T; and (11) S∼Tγlinking durations with averaged sizes, S(T). Note that these last exponents (the spreading and avalanche ones) are not independent of each other, but related through scaling relations as derived in, e.g., [55] τ=1+η+2δ 1+η+δ,(14) τt=1+δ, (15) γ=τt−1 τ−1=1+δ+η, (16) where the last one describes the “crackling noise” scaling relation [57]. Other scaling relations can be found in [4,5,55], in particular, θ=β/ν,(17) relates static and dynamic exponents. Importantly, for standard processes with absorbing states (e.g., DP and TDP), the averaged shape of avalanches with different durations and sizes (or “mean temporal profile of avalanches”) collapses onto a universal curve that typically has a symmetric inverted parabolic form (see Sec. V)[58,59]. It is noteworthy that there is a set of exponents that can be argued to remain unchanged across transition types (a fact that is also confirmed numerically [60]). This is due to the mean-field and diffusive character shared by all the transitions discussed here. In particular, owing to the diffusive nature of the system in all continuous transitions, correlations (ξ) should diverge at the critical point with the well-known meanfield exponents [1,4,5] as follows: ξ⊥∝ν⊥with ν⊥=1/2, for the correlation length; and ξ∝νwith ν=1, for the time correlation. From this, given that [55]z=2ν⊥/ν= 1, z=1 for all continuous transitions here. Similarly, the survival probability exponent (whose scaling behavior was determined in [61]) always takes a value δ=1 for all the continuous transitions studied here, implying that τt=2 [see Eq. (15)] is conserved across transitions. Finally, the exponent ηis expected to vanish for all mean-field transitions (for which there is no “anomalous dimension” [3]). However, remarkably, here we report on a possible exception to this general rule (η=2) for some of the “anomalous” transitions (see below). Finally, it is worth stressing that many aspects of the transitions, especially those related to discontinuous or hybrid-type transitions, are not expected to exhibit universal features. A. T1: Directed percolation (transcritical bifurcation) T1corresponds to a transcritical bifurcation, describing a continuous transition between the standard quiescent and active phases. As discussed in Sec. I, guided by universality principles, one expects it to lie in the usual (mean-field) directed percolation universality class (DP) [1,3,4,6,7]. Indeed, this is the case, as explicitly shown in what follows. Transcritical bifurcations occur when the quiescent steady state loses its local stability, i.e., Eq. (8). Expanding Eq. (3) and Eq. (4) in power series of E∗and I∗, one finds E∗(;h=0) =(α+wII)3 (α+wII)3−wEIw2 IE +O(2),(18) from which β=1 follows (see Fig. 4). The introduction of an external field hsmooths out the transition [as illustrated with dashed lines in Fig. 4(a)]. Hence, expanding the fixed point in powers of h,at=0, yields E∗(h;=0) =(α+wII)2(wEI −α)h αwIEw2 EI −(wII +α)3+O(h),(19) so that δh=2[seeFig.5(a)]. Similarly, one can derive the solution for I(t) and, by expanding it in a power series, obtain I(t)≈[wIE/(wII + α)]E(t). It is, thus, convenient to define two new variables: and , as the weighted linear combinations: 2=wIE E+(wII +α)I,(20) 2=wIE E−(wII +α)I,(21) in terms of which the mean-field dynamics [Eq. (3) and Eq. (4)] take a simpler form: 2˙ ≈ +(2wEE +2wII +)+O,(22) 2˙ ≈ −[2(wEE +wII −2α)+]+O,(23) where O=O(2, 2,,...) stands for higher-order terms. Observe that, right at the transition (=0), the stability matrix around the origin is A=0wT EE +wII 0wT EE −wII −2α,(24) which has a vanishing eigenvalue, while the second one is strictly negative at criticality for T1transitions. This means that decays exponentially fast; therefore, it is an irrelevant field for scaling. Only one “slow mode” or “relevant field” exists, , and—as theoretically predicted in Grinstein et al. [10] for these conditions—the scaling behavior should coincide with standard DP. In particular, at the transition point—where the linear term of Eq. (22) vanishes—the quadratic term dominates and therefore (t)∝t−1, so that θ=1 [as numerically confirmed in Fig. 6(a)]. Considering the previous three independent exponent values, one can already conclude that the T1transition actually belongs in the DP universality class (see Table I). Nevertheless, for the sake of completeness, the survival probability and the total number of particles, at =0, are numerically verified to scale with spreadingexponent values δ=1 [Fig. 7(a)] and η=0 [Fig. 8(a)], respectively, as expected for the DP class. We have also confirmed the consistency with DP by numerically analyzing the statistics of avalanches at the transition, revealing exponent values compatible with the DP predictions τ=3/2, 034110-6 UNCONVENTIONAL CRITICALITY, SCALING … PHYSICAL REVIEW E 108, 034110 (2023) FIG. 4. Order parameter as a function of the control parameter, =wEE −w∗ EE, across the eight possible types of transitions as represented in Fig 3. In all plots revealing continuous transitions, the dot-dashed blue curves correspond to asymptotic behavior E∗(;h=0) ∼β, with the corresponding values of the βexponent. The shaded gray areas correspond to bistability between active and standard quiescent phases and the pink shaded area between active and excitable quiescent phases. The green shaded areas represent the excitable quiescent phase (same colors as Fig. 3). (a) At T1(wIE =3), the system exhibits a second-order transition from the standard quiescent to the active phase, consistently with the directed-percolation (DP) universality class (E∗∼1for ⩾0). (b) T2(wIE =3) is also a continuous phase transition occurring through a tricritical point and is consistent with the tricritical directed percolation (TDP) universality class (E∗∼1/2). (c, d) In contrast, T3and T4are first-order or discontinuous phase transitions with coexistence between an active phase and one of two possible kinds of quiescence (wIE =3): first, a standard quiescent state (gray shaded area) and, second, an excitable quiescent state (pink shaded areas). (e) Case B(wIE =1) allows for a special tricritical transition (T5) occurring through a Hopf-tricritical point (E∗∼1/2). (f–h) In case C (wIE =0.8), both T6and T7are continuous phase transitions when h=0, with E∗∼1and T8,E∗∼1/2, respectively. However, once a nonvanishing external field h= 0 is introduced (dash-dotted and dashed lines), there is bistability driven by the external field [for more details see gray shaded areas in Figs. 5(f)–5(h)]. Parameter values are set as in Fig 3. FIG. 5. Order parameter as a function of the external field right at the transition (=0). Observe that three out of the eight types of transitions described here exhibit power-law scaling with the external field, i.e., E∗(h;=0) ∼h1/δh.(a)ForT1,δh=2, consistently with the DP universality class. (b) For T2,δh=3, consistently with TDP. (c, d) For T3and T4, the order parameter’s response to an external field shows bistability (shaded area). (e) Transition T5(H+TDP), differently from the usual tricritical transition (TDP), scales with δh=2. (f)–(h) Remarkably, for T6,T7,andT8, contrary to the behavior with h=0, the order parameter becomes bistable (shaded area) as the external field increases. Parameter values as in Fig. 3. 034110-7 PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) FIG. 6. Order parameter time series for the eight types of transitions. Observe that only three of them exhibit dynamical scaling E(t)∝t−θ. (a) T1exhibits an asymptotic power-law decay with the expected DP value θ=1. (b) T2shows a slower asymptotic time decay in the TDP class, θ=1/2. (c) For T3, the saddle-node bifurcation gives rise to bistability between an active and a quiescent phase. (d) T4behaves very similarly to T3but frustrated oscillations drive the system more easily to regions II and III, so the bistability is between active and excitable quiescent phases. (e) T5is a genuine second-order phase transition with θ=1. (f)–(h) T6,T7,andT8are not genuine continuous transitions and show no signatures of dynamic scaling, but rather an exponential decay to quiescence. Parameter sets as in Fig. 3. τt=2, and γ=2 [see the distributions of sizes S, durations T, and average sizes as a function of durations in Figs. 9(a),9(d), and 9(g), respectively]. Moreover, the averaged avalanche shape is approximately an inverted parabola throughout the T1line [Figs. 10(a) and 10(b)], collapsing for different durations with γ=2, even if with some asymmetry (see Sec. Vfor a more in-depth discussion on avalanche shapes). Thus, in summary, at the line of transcritical bifurcations (T1) separating a standard quiescent from the active phase, the Wilson-Cowan stochastic model exhibits a genuine critical point in the DP class, a result that is consistent with recent analyses of de Candia et al. [33] for their specific choice of parameter values. TABLE I. Summary of mean-field exponents for the discussed continuous phase transitions [13,55]. DP TDP H+TDP Codim. 1 2 3 β11/21/2 δh23 2 θ11/21 δ11 1 η00 2 ν11 1 τ3/23/25/4 τt22 2 γ22 4 B. T2: Tricritical directed percolation (saddle-node-transcritical bifurcation) The tricritical point in case A [see Fig. 3(a)] corresponds to a saddle-node transcritical (SNT) bifurcation—i.e., where the lines of transcritical and saddle-node bifurcations intersect without further degeneracies [62,63]. Thus, in order to tune to this transition point one needs to set two parameters in the phase diagram (wEE,wEI), as explicitly calculated in Appendix B. An analysis in terms of the fields and  (analogously to the previous case) shows that there is only one vanishing eigenvalue at the transition, and, thus, the second field is irrelevant for scaling. Therefore, T2is expected to be described by the mean-field tricritical directed percolation universality class (TDP) [13]. Indeed, considering the leadingorder terms in a power expansion in both and h, one has E∗(, h=0) ≈ wIE +O(),(25) E∗(h,=0) ≈3w2 IE −(α+wII)2h w2 IE[(α2−3)wIE −(α2+3)α]1 3 +O(h1 2),(26) from where β=1/2 [Fig. 4(b)] and δh=3 [Fig. 5(b)], as expected for the TDP universality class. At the transition, the lowest order correction of Eq. (22) in is O(3), so that asymptotically ∝t−1/2and, hence, θ=1/2, as numerically confirmed in Fig. 6(b). Once again, considering the linear relationship between E(t) and I(t), both densities share this scaling. 034110-8 UNCONVENTIONAL CRITICALITY, SCALING … PHYSICAL REVIEW E 108, 034110 (2023) FIG. 7. Survival probability as a function of time. The survival probability at second-order phase transitions scales as Ps(t)∝t−δ. Black dots stand for the numerical simulations (for the same parameters as Fig. 3and N=108) and dashed lines show the corresponding exponent value. (a) T1belongs to the DP universality class, i.e., δ=1. (b) T2belongs to TDP so that δ=1 (the same as DP). (c), (d) For the first-order phase transitions, the system’s survival probability converges to a nonvanishing value as t→∞due to the possibility of being attracted to the active phase. (e) For T5,δ=1 as in the previous continuous transitions, but with stronger finite-size effects. (f) For T6, the system shows a behavior similar to T5: a decay with δ=1 and strong finite-size effects. (g), (h) The survival probability shows a peculiar behavior of several sharp decays with some small plateaus, which stem from the excitability of the quiescent phase. Finally, the exponent for the survival probability remains δ=1[seeFig.7(b)], η=0 [see Fig. 8(a)], τ=3/2 [Fig. 9(b)], τt=2 [Fig. 9(e)], and γ=2 [Fig. 9(h)], all of which are consistent with the expected values in the TDP class (see Table I). C. T3: Standard discontinuous transition (saddle-node bifurcation) The line of saddle-node bifurcations [see Fig. 3;Eq.(11)] defines the third type of transition, T3, to go from a standard quiescent state to the active phase. This type of transition is FIG. 8. Mean number of particles N(t) in spreading experiments in bona fide continuous phase transitions (i.e., T1,T2,andT5), at which one expects N(t)∼tη. Simulations with the same parameters as Fig. 3with N=108(a) and N=108and N=1010 (b). (a) For T1and T2, we obtain results compatible with η=0, as expected for DP and TDP as well as, in general, for mean-field theories. (b) On the other hand, for T5we obtain the unusual result η=2, with strong finite-size effects. characterized by a discontinuous jump in the order parameter and includes an intermediate regime of bistability, where both the active and the standard quiescent state are stable [see Figs. 3(a) and 4(c)]. The regime of coexistence lasts until, at a second bifurcation, the quiescent phase loses its local stability. Given that the transition is discontinuous, the exponents β and δhare not properly defined [Fig. 5(c)]. Similarly, neither the activity nor the survival probability decay to 0 for initial conditions in the basin of attraction of the active phase [see Figs. 6(c) and 7(c)], so that the exponents θand δare not well defined either. Thus, in summary, the T3transition is a standard firstorder or discontinuous transition into a quiescent or absorbing state [5,13]. D. T4: Discontinuous transition from an excitable quiescent state (saddle-node bifurcation) A scenario very similar to T3occurs at T4, which appears in all three possible phase diagrams (A, B, and C; see Fig. 3). Transition T4is also discontinuous with phase coexistence, but it differs from T3in the fact that—as illustrated in Fig. 4(d)— the quiescent phase that coexists with the active one in the regime of bistability is of the excitable type, rather than the standard one. The bistability regime ends where the excitable quiescent phase loses its global stability in favor of the active one, at a discontinuity-induced transition (black-dashed line in Fig. 3). For the same reasons as in T3, none of the critical exponents is well defined [see Figs. 5(d),6(d), and 7(d)]. Thus, in summary, T4is a discontinuous transition with bistability, but with the peculiarity of having an excitable quiescent state coexisting with the active one. 034110-9 PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) and T8). In cases A and C, this bifurcation has codimension 2 and occurs at wSNT EI =(α+wII)3 w2 IE ,(B3) wSNT EE =α+(α+wII)2 wIE .(B4) The nontrivial solution emerges from the trivial solution with wEE, and it scales with the distance to the critical value, ,as E∗∝I∗∝1/2. Finally, in Fig. 3case B, a codimension 3 bifurcation emerges from an extra fine-tuning of the parameters when wIE =α+wII. For this choice of parameters, at wEI =wIE, the saddle-node transcritical collides with the Hopf right at the tricritical point, T5. Combining Eq. (12) and Eq. (B3), the values of the control parameters, for this bifurcation, are wEI =α+wII,(B5) wEE =2α+wII.(B6) APPENDIX C: DO TRAJECTORIES CROSS OR SLIDE ONTO THE SWITCHING MANIFOLDS? Piecewise continuous dynamics have two possible behaviors at the switching manifolds: sliding or crossing [48]. To determine the behavior of the Wilson-Cowan model system, we consider the Heaviside function (1)inEq.(3) and Eq. (4): ˙x=⎧ ⎪ ⎨ ⎪ ⎩ f+ x≡−αx+(1 −x) tanh(wix−wjy), if s≡wix−wjy>0 f− x≡−αx,ifs<0 ,(C1) where f+ x(f− x) is evaluated to the right (left) of the switching manifold, s=0. Let us consider the switching manifold sE=0, where E= (wEI/wEE )I[Fig. 1(a)]. One can then write  ∇sE=∂ ∂EsE ∂ ∂IsI=wEE −wEI,(C2)  f+=−αE+(1 −E) tanh (wEEE−wEII) −αI+(1 −I) tanh (wIEE−wIII)T ,(C3)  f−=−αE −αI+(1 −I) tanh (wIEE−wIII)T .(C4) When the system reaches the switching manifold, the trajectories will cross it or slide on it depending on the sign of ( f+· ∇sE)(  f−· ∇sE) at the switching manifold:  f+· ∇sE=−α(wEEE−wEII) +wEE(1 −E) tanh (wEEE−wEII) −wEI(1 −I) tanh (wIEE−wIII),(C5)  f−· ∇sE=−α(wEEE−wEII) −wEI(1 −I) tanh (wIEE−wIII).(C6) Given that at the switching manifold, wEEE=wEII:  f+· ∇sE=−wEI(1 −I) ×tanh wIE(wEI −wII) wEE I,(C7)  f−· ∇sE=−wEI(1 −I) ×tanh wIE(wEI −wII) WEE I,(C8) ( f+· ∇sE)(  f−· ∇sE)=[wEI(1 −I) ×tanh wIE(wEI −wII) wEE I2 . (C9) For (  f+· ∇sE)(  f−· ∇sE)>0, the trajectories cross the switching manifold and cannot cross back. Observe that, when wEE/wEI <wIE/wII, the flow points to region II, making it a trapping region. This condition is always satisfied for the parameters we explored in our simulations. When wEE/wEI >wIE/wII, region II amounts to sE>0 and sI<0 and the calculations are analogous. However, the flow can be reversed, so that a trajectory beginning in region II may cross to region I, and the latter becomes the trapping region instead. [1] J. Marro and R. Dickman, Nonequilibrium Phase Transition in Lattice Models (Cambridge University Press, Cambridge, 1999). [2] H. Hinrichsen, Adv. Phys. 49, 815 (2000). [3] G. Grinstein and M. A. Muñoz, The Statistical Mechanics of Absorbing States, Lecture Notes in Physics, Vol. 493 (Springer, Berlin, 1997), p. 223. [4] M. Henkel, H. Hinrichsen, and S. Lübeck, Non-Equilibrium Phase Transitions: Absorbing phase transitions, Theoretical and Mathematical Physics (Springer, Berlin, 2008). [5] G. Ódor, Universality in Nonequilibrium Lattice Systems: Theoretical Foundations (World Scientific, Singapore, 2008). [6] H.-K. Janssen, Z. Phys. B 42, 151 (1981). [7] P. Grassberger, in Nonlinear Phenomena in Chemical Dynamics, edited by C. Vidal and A. Pacault (Springer, Berlin, 1981), pp. 262–262. [8] J. Binney, N. Dowrick, A. Fisher, and M. Newman, The Theory of Critical Phenomena (Oxford University Press, Oxford, 1993). [9] M. A. Muñoz, G. Grinstein, R. Dickman, and R. Livi, Phys. Rev. Lett. 76, 451 (1996). [10] G. Grinstein, Z.-W. Lai, and D. A. Browne, Phys.Rev.A40, 4820 (1989). [11] T. E. Harris, The Theory of Branching Processes (Courier Corporation, Berlin, 2002). [12] T. Liggett, Interacting Particle Systems, Classics in Mathematics (Springer, Berlin, 2004). 034110-16 UNCONVENTIONAL CRITICALITY, SCALING … PHYSICAL REVIEW E 108, 034110 (2023) [13] S. Lübeck, J. Stat. Phys. 123, 193 (2006). [14] P. Villa Martín, J. A. Bonachela, and M. A. Muñoz, Phys. Rev. E89, 012145 (2014). [15] V. R. V. Assis and M. Copelli, Phys. Rev. E 80, 061105 (2009). [16] J. M. Beggs and D. Plenz, J. Neurosci. 23, 11167 (2003). [17] T. Mora and W. Bialek, J. Stat. Phys. 144, 268 (2011). [18] D. R. Chialvo, Nat. Phys. 6, 744 (2010). [19] M. A. Muñoz, Rev. Mod. Phys. 90, 031001 (2018). [20] D. Plenz, T. L. Ribeiro, S. R. Miller, P. A. Kells, A. Vakili, and E. L. Capek, Front. Phys. 9, 639389 (2021). [21] D. R. Chialvo, Acta Phys. Pol. B 49, 1955 (2018). [22] J. O’Byrne et al.,Trends Neurosci. 45, 820 (2022). [23] A. J. Fontenele, N. A. P. de Vasconcelos, T. Feliciano, L. A. A. Aguiar, C. Soares-Cunha, B. Coimbra, L. Dalla Porta, S. Ribeiro, A. J. Rodrigues, N. Sousa et al.,Phys.Rev.Lett.122, 208101 (2019). [24] A. Ponce-Alvarez, A. Jouary, M. Privat, G. Deco, and G. Sumbre, Neuron 100, 1446 (2018). [25] D. Millman, S. Mihalas, A. Kirkwood, and E. Niebur, Nat. Phys. 6, 801 (2010). [26] M. Martinello, J. Hidalgo, A. Maritan, S. di Santo, D. Plenz, and M. A. Muñoz, Phys. Rev. X 7, 041071 (2017). [27] J. M. Cortes, M. Desroches, S. Rodrigues, R. Veltz, M. A. Muñoz, and T. J. Sejnowski, Proc. Natl. Acad. Sci. USA 110, 16610 (2013). [28] M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, PLoS Comput. Biol. 6, e1000846 (2010). [29] O. Kinouchi and M. Copelli, Nat. Phys. 2, 348 (2006). [30] M. Girardi-Schappo, E. F. Galera, T. T. Carvalho, L. Brochini, N. L. Kamiji, A. C. Roque, and O. Kinouchi, J. Phys. Complex. 2, 045001 (2021). [31] M. Girardi-Schappo, L. Brochini, A. A. Costa, T. T. A. Carvalho, and O. Kinouchi, Phys.Rev.Res.2, 012042(R) (2020). [32] T. T. A. Carvalho, A. J. Fontenele, M. Girardi-Schappo, T. Feliciano, L. A. A. Aguiar, T. P. L. Silva, N. A. P. de Vasconcelos, P. V. Carelli, and M. Copelli, Front. Neural Circuits 14, 576727(2021). [33] A. De Candia, A. Sarracino, I. Apicella, and L. de Arcangelis, PLoS Comput. Biol. 17, e1008884 (2021). [34] M. K. Nandi, A. Sarracino, H. J. Herrmann, and L. de Arcangelis, Phys.Rev.E106, 024304 (2022). [35] I. Apicella, S. Scarpetta, L. de Arcangelis, A. Sarracino, and A. de Candia, Sci. Rep. 12, 21870 (2022). [36] R. Corral López, V. Buendía, and M. A. Muñoz, Phys. Rev. Res. 4, L042027 (2022). [37] H. R. Wilson and J. D. Cowan, Biophys. J. 12, 1 (1972). [38] E. Wallace, M. Benayoun, W. van Drongelen, and J. D. Cowan, PLoS ONE 6, 1 (2011). [39] Y. Maruyama, Y. Kakimoto, and O. Araki, Biol. Cybern. 108, 355 (2014). [40] J. D. Cowan, J. Neuman, and W. van Drongelen, J. Math. Neurosci. 6, 1 (2016). [41] F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks, Applied Mathematical Sciences, Vol. 126 (Springer, New York, 1997). [42] R. M. Borisyuk and A. B. Kirillov, Biol. Cybern. 66, 319 (1992). [43] Let us remark that a very similar model has been recently proposed to generalize the standard contact process to include inhibitory units and is described by similar equations [36]. [44] N. Brunel, J. Comput. Neurosci. 8, 183 (2000). [45] N. V. Kampen, Stochastic Processes in Physics and Chemistry (North Holland, 2007). [46] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 3rd ed., Springer Series in Synergetics, Vol. 13 (Springer-Verlag, Berlin, 2004). [47] T. Ohira and J. D. Cowan, in Mathematics of Neural Networks, edited by S. W. Ellacott, J. C. Mason, and I. J. Anderson (Kluwer Academic Publishers, Norwell, 1997), pp. 290–294. [48] P. Glendinning and M. R. Jeffrey, An Introduction to Piecewise Smooth Dynamics (Birkhäuser, Cham, 2019). [49] J. Harris and B. Ermentrout, SIAM J. Appl. Dyn. Syst. 14,43 (2015). [50] M. Kunze, Non-Smooth Dynamical Systems, Lecture Notes in Mathematics, Vol. 1744 (Springer Science & Business Media, Berlin, 2000). [51] E. M. Izhikevich, Dynamical Systems in Neuroscience (MIT Press, Cambridge, 2007). [52] S. H. Strogatz, Nonlinear Dynamics and Chaos with Student Solutions Manual: With Applications to Physics, Biology, Chemistry, and Engineering (CRC Press, Boca Raton, 2018). [53] I. I. Lima Dias Pinto and M. Copelli, Phys. Rev. E 100, 062416 (2019). [54] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz, J. Stat. Mech. (2018) 073402. [55] M. A. Muñoz, R. Dickman, A. Vespignani, and S. Zapperi, Phys. Rev. E 59, 6175 (1999). [56] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models, Collection Alea-Saclay: Monographs and Texts in Statistical Physics (Cambridge University Press, Cambridge, 1999). [57] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature (London) 410, 242 (2001). [58] N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E. Lee DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler, Phys. Rev. Lett. 108, 208102 (2012). [59] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz, Phys. Rev. E 95, 032115 (2017). [60] N. Marshall, N. M. Timme, N. Bennett, M. Ripp, E. Lautzenhiser, and J. M. Beggs, Front. Physiol. 7, (2016). [61] M. A. Muñoz, G. Grinstein, and Y. Tu, Phys. Rev. E 56, 5101 (1997). [62] L. van Veen and M. Hoti, Int. J. Bifurcation Chaos 29, 1950104 (2019). [63] L. Lai, Z. Zhu, and F. Chen, Mathematics 8, 1280 (2020). [64] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature (London) 592, 363 (2021). [65] J. D. Noh and H. Park, Phys.Rev.Lett.94, 145702 (2005). [66] K. Christensen, H. C. Fogedby, and H. Jeldtoft Jensen, J. Stat. Phys. 63, 653 (1991). [67] A. Chessa, H. E. Stanley, A. Vespignani, and S. Zapperi, Phys. Rev. E 59, R12 (1999). [68] R. Dickman and J. M. M. Campelo, Phys. Rev. E 67, 066111 (2003). [69] A. Deluca and Á. Corral, Acta Geophysica 61, 1351 (2013). 034110-17 PIUVEZAM, MARIN, COPELLI, AND MUÑOZ PHYSICAL REVIEW E 108, 034110 (2023) [70] S. Papanikolaou, F. Bohn, R. L. Sommer, G. Durin, S. Zapperi, andJ.P.Sethna,Nat. Phys. 7, 316 (2011). [71] S. R. Miller, S. Yu, and D. Plenz, Sci. Rep. 9, 1 (2019). [72] L. Laurson, X. Illa, S. Santucci, K. Tore Tallakstad, K. J. Måløy, and M. J. Alava, Nat. Commun. 4, 2927 (2013). [73] E. Gudowska-Nowak, M. A. Nowak, D. R. Chialvo, J. K. Ochab, and W. Tarnowski, Neural Comput. 32, 395 (2020). [74] J. Hidalgo, L. Seoane, J. Cortés, and M. Muñoz, PLoS ONE 7, e40710 (2012). [75] J. Almeira, T. S. Grigera, D. R. Chialvo, and S. A. Cannas, Phys. Rev. E 106, 054140 (2022). [76] V. Buendía, P. Villegas, S. di Santo, A. Vezzani, R. Burioni, and M. A. Muñoz, Sci. Rep. 9, 15183 (2019). [77] M. A. Muñoz, R. Juhász, C. Castellano, and G. Ódor, Phys. Rev. Lett. 105, 128701 (2010). [78] P. Moretti and M. A. Muñoz, Nat. Commun. 4, 2521 (2013). [79] G. Odor, Phys.Rev.E94, 062411 (2016). [80] D. T. Gillespie, J. Comput. Phys. 22, 403 (1976). [81] D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977). 034110-18