Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic M/EEG results

We revisit the CROS (“CRitical OScillations”) model which was recently proposed as an attempt to reproduce both scale-invariant neuronal avalanches and long-range temporal correlations. With excitatory and inhibitory stochastic neurons locally connected in a two-dimensional disordered network, the model exhibits a transition where alpha-band oscillations emerge. Precisely at the transition, the fluctuations of the network activity have nontrivial detrended fluctuation analysis (DFA) exponents, and avalanches (defined as supra-threshold activity) have power law distributions of size and duration. We show that, differently from previous results, the exponents governing the distributions of avalanche size and duration are not necessarily those of the mean-field directed percolation universality class (3/2 and 2, respectively). Instead, in a narrow region of parameter space, avalanche exponents obtained via a maximum-likelihood estimator vary continuously and follow a linear relation, in good agreement with results obtained from M/EEG data. In that region, moreover, the values of avalanche and DFA exponents display a spread with positive correlations, reproducing human MEG results.


Introduction
The critical brain hypothesis [1] has emerged in the last decades as a potential framework for theoretically addressing many intriguing questions that have challenged neuroscientists. In light of the nonlinear nature of individual neuronal dynamics and our rudimentary understanding of the collective phenomena that emerge when they interact, these questions include, for instance, the issue of segregation versus integration, optimization of dynamical repertoire, and response to external stimuli, among others (see e.g. Refs. [1][2][3] for reviews). But if the brain is indeed critical, what is the phase transition? In the following, we briefly review our assessment of the current status of the field, highlight the limitations of a class of models that has dominated the literature in the last decade, and offer an alternative analysis of a more recent model that copes better with experimental results.
A seminal work by Beggs and Plenz in 2003 reported neuronal avalanches experimentally recorded in vitro [4], lending support to the criticality hypothesis and effectively laying the groundwork for a vast field of research that has involved neuroscientists and physicists alike [5]. In their original setup, spontaneous local field potentials (LFPs) were recorded from cultured slices of the rat brain [4]. The term "neuronal avalanche" seemed like a natural choice, since the observed bursts of suprathreshold activity were interspersed by silence and showed a clear separation of time scales (their duration being much shorter than the interavalanche intervals). Moreover, their size s, defined as either the number of electrodes with suprathreshold activity or the sum of the potentials, were shown to be statistically distributed according to a power law, P(s) � s −1.5 [4]. Such scale-invariant statistics are one of the hallmarks of a critical system [6,7].
Other signatures of criticality in the brain have been proposed, such as long-range temporal correlations, which have been observed both at macroscopic and microscopic levels. For instance, detrended fluctuation analysis (DFA) [8] performed in electroencephalographic and magnetoencephalographic signals have revealed temporal correlations and power-law scaling behaviour in spontaneous oscillations of the normal human brain during large time scales [9,10]. At a much smaller scale, a similar technique showed that spike avalanches recorded in different regions of freely-behaving rats also have long-range temporal correlations [11].
A natural next step from the modeling perspective would be to conjugate both ideas. The search for a model that can produce both power-law distributed avalanches and long-range temporal correlations, however, is likely to face some theoretical challenges. Let us briefly discuss why.
It is important to remember that one of the appeals of the 3/2 exponent experimentally observed by Beggs and Plenz is that it coincides with the critical exponent for a branching process [12], which has therefore become a theoretical workhorse in the field. Or, if one extends the idea to rather general networks [13], 3/2 coincides with the avalanche size critical exponent of any model belonging to the universality class of directed percolation (DP) in dimension larger than or equal to its upper critical dimension (d c = 4). In other words, τ = 3/2 is a mean- field exponent for avalanche-size distributions of a class of models in which a continuous phase transition occurs between an absorbing and an active phase [7,14]. A minimum model of this universality class consists of a large (ideally infinite) number of units ("neurons" or "regions of interest") which can be either "on" (active) or "off" (inactive) [7]. If a unit is "off", it may switch to an "on" state with a probability that increases with the number of "on" neighbors and some coupling, say, λ. If a unit is "on", it spontaneously switches to "off" at some constant rate (which accounts for the typical duration of a spike, for instance; variants may include additional states to model the refractory period [15]). The absorbing phase bears its name because, for sufficiently small λ, eventually all units go to the "off" state and the dynamics of the whole network stops. For sufficiently large λ, on the other hand, propagation of the "on" state among units occurs at such a rate that the system stays in an active phase, characterized by stable self-sustained activity, i.e. nonzero time-and ensemble-average density of active sites hρi (the typical order parameter for these models). The boundary in parameter space between those two qualitatively different regimes is the critical point λ = λ c , above which hρi departs continuously from zero as hρi � (λ − λ c ) β , where β is a critical exponent [7]. Precisely at λ = λ c , the system is on the verge of displaying self-sustained activity, so perturbations to the absorbing state (which is still stable) do not have a characteristic time to die out or a characteristic size. In fact, their size s (defined as the number of active sites along the excursion) and duration T (defined as the time between the first and last active site in between epochs of complete quiescence) are power-law distributed at the critical point: P(s) � s −τ and PðTÞ � T À t t . These critical exponents depend on network dimensionality D [14], and for D � 4 the mean-field exponents are precisely τ = 3/2 and τ t = 2, as observed by Beggs and Plenz [4]. The importance of this type of phase transition cannot be overstated. A conjecture put forward by Janssen [16] and Grassberger [17] states that a model with a continuous transition to a single absorbing state and no further symmetries (parity conservation etc) should fall into the DP universality class. In other words, "continuous transitions to an absorbing state belong generically to the directed percolation universality class" [7].
In the model, these perturbations are called "avalanches", among other reasons, because their dynamics are subject to an infinite separation of time scales by construction: once an avalanche is over (all sites "off"), the next one will not start unless the system is arbitrarily perturbed again. This unassuming detail, however, poses a challenge if one wants to explore this class of models to reproduce the above mentioned long-range temporal correlations which were observed experimentally [9,11]. Since consecutive avalanches are, by definition, separated by returns to the absorbing state, it is not apparent how inter-avalanche correlations could emerge (self-organizing mechanisms are a potential candidate, yet to be tested [18]).
Given this state of affairs, an interesting perspective was put forward by the CROS ("CRitical OScillations") model proposed by Linkenkaer-Hansen et al [19,20]. The model has both excitatory and inhibitory stochastic spiking neurons arranged in an L × L square lattice, each neuron interacting with a controlled random fraction of its neighbors within a small ℓ × ℓ region (Fig 1A, see Methods for details). All excitatory neurons receive an independent constant Poisson input which ensures the model formally does not have an absorbing state. Changing the model parameters controlling excitation and inhibition (r E and r I ), Linkenkaer-Hansen et al reported a transition where alpha-band oscillations appear, as marked by the emergence of a peak in the Fourier spectrum of the global network activity. At the transition line in the (r E , r I ) plane, long-range temporal correlations were observed with DFA exponents α ' 1, in contrast with uncorrelated fluctuations far from the transition (α ' 0.5). Remarkably, at the same transition line Linkenkaer-Hansen et al also managed to obtain avalanches with power law distributions of size [19] and duration [20], in either case with the respective meanfield exponent τ = 3/2 and τ t = 2. Importantly, in the absence of an absorbing state, their definition of an avalanche required the imposition of an ad hoc threshold θ, defined as a fraction Γ of the median network activitym. In the CROS model, therefore, it is the crossing of an arbitrary threshold upwards and downwards that marks respectively the beginning and end of an avalanche (see Fig 2A and Methods). Differently from a DP-like model with an absorbing state, the system dynamics itself controls when the threshold will be crossed, raising the question whether that is the main ingredient giving rise to long-range temporal correlations.
This thresholding procedure is often deemed artificial in DP-like models with an absorbing state, among other reasons because it introduces an additional parameter on which results depend, namely, the threshold itself. However, it is important to keep in mind that, when dealing with empirical data, the avalanche definition inevitably requires the introduction of a threshold (albeit usually a local one, differently from the CROS model). When avalanches are understood as the spatio-temporal patterns of "on" sites between "silences", as discussed above, a local threshold is required to operationally define what "on" means in the first place. This issue has been present since the first results emerged for LFPs [4], and reappear in M/ EEG [21][22][23] and fMRI studies [24]. Once one has a point process with discrete events, a second parameter is needed to determine what "silence" means: the duration of the bin with which the time series is parsed. The important point here is that exponents governing avalanche distributions in principle depend on these parameters and often show high variability [21][22][23]. Frequently, an additional criterion is used to choose between a range of values. For instance, if one works under the hypothesis that the phase transition governing brain criticality is that of a branching process, then one can require the branching ratio to be near its critical value of one [4,21]. On the one hand, when found together, a critical branching ratio and a critical exponent τ = 3/2 are results that consistently reinforce each other from a theoretical point of view. On the other hand, this consistency is constrained by the initial assumption of the underlying universality class. Relaxing this assumption and exploring the full variability of the exponents observed experimentally provides a challenge to any model.
The results of the model proposed by Linkenkaer-Hansen et al are remarkable in that they seem to reconcile the emergence of long-range temporal correlations with avalanche distributions which in principle appear in models of a very different nature, as discussed above. The reconciliation is even more appealing given that the coexistence of these two phenomena has been reported experimentally in M/EEG recordings [21,22]. However, several questions remain, which we set out to investigate here.
First, we delve deeper into the following question: to which extent is mean-field DP a good model for the transition observed in the CROS model? To do so, we introduce a tentative order parameter to better characterize the putative phase transition, while also simulating larger system sizes (L = 300) than in the original results (L = 50) [19]. We also point out that it is highly counterintuitive to have a two-dimensional model exhibiting mean-field exponents. That could be due to an insufficiently small ratio ℓ/L (which is addressed here, again, by increasing L). We argue that, if the onset of oscillations in a two-dimensional network is to be reconciled with DP critical exponents, then the DP exponents for D = 2 should also be tested.
Second, we explore parameter space in more detail around the transition line, relaxing the assumption of a DP universality class, allowing the exponents to be adjusted by a maximumlikelihood estimator (MLE) and testing whether the model satisfies other scaling relations. In so doing, we subject the putative phase transition of the CROS model to more stringent tests [25,26]. Moreover, by not imposing τ = 3/2 and τ t = 2, the model reveals exponents that vary continuously within a transition region, opening the door for comparison with experimental results.
Finally, we revisit some experimental results from M/EEG recordings where both avalanche and DFA exponents have shown high variability [22,23]. We explore the extent to which the model is able to reproduce the spread of exponents and the correlations among them. The size of an avalanche can be defined as s g , the total number of spikes, or s θ , the number of spikes minus the threshold value. The avalanche duration T is the time that the fluctuations stay above θ. (B) Detail of the power spectrum around the region where a peak appears. The order parameter φ is given by the ratio between the power peak area ϕ p and the total area ϕ p + ϕ u . Red symbols represent the frequencies that bound the φ area: f min (diamond), f peak (circle) and f max (triangle). https://doi.org/10.1371/journal.pcbi.1006924.g002

CROS network model
We employed essentially the same CROS (CRitical OScillations) model as proposed by Poil et al [19]. Excitatory (75%) and inhibitory (25%) neurons were randomly arranged in a bidimensional L × L square lattice, with open boundaries (Fig 1A). Each neuron has a finite local range of connections limited to a square of size ℓ × ℓ centered around it (Fig 1A). We used ℓ = 7 like in the original model [19], which means that a typical neuron (far away from the borders) could connect to a maximum of 48 other neurons in their neighborhood. Excitatory (r E ) and inhibitory (r I ) connectivity are the two free parameters of this model, being set between 2-60% and 30-90%, respectively, of the total number of neurons within the local range. Each neuron sends its outgoing synapses to other neurons at a distance r with probability PðrÞ ¼ Ce À r=a 0 , where a 0 = 1 is the nearest-neighbor distance, and C a normalization constant such that ∑ neighbors P(r) equals r E or r I for presynaptic excitatory and inhibitory neurons, respectively. Notice that the networks have disorder, in the sense that, for a given point (r E , r I ) in parameter space, different synaptic connections are possible (corresponding to how the specific grey neurons in Fig 1A are randomly selected within the ℓ × ℓ interaction range). We will return to this point in the discussion of the results. Five networks for each combination of r E and r I were created and simulated for 2 20 time steps (dt) of 1 ms. Each connection had a weight (W ij ) depending on the nature (excitatory or inhibitory) of the presynaptic and postsynaptic neurons ( Fig 1B).
The neuron dynamics starts with a neuron i integrating all the received input from the connected neighbors and updating its synaptic current I i , which is also subject to an exponential synaptic decay: In Eq (1), S(t) 2 {0, 1} N is a binary vector denoting whether the presynaptic neurons fired in the previous time step, whereas the connection weights are fixed at W EI = 0.011, W EE = 0.02 and W IE = W II = −2 and the time constant is τ I = 9 ms [19]. Given the synaptic current I i we compute the spiking rate variable R Si , which is also subject to an exponential decay: where τ P (excitatory) = 9 ms, τ P (inhibitory) = 12 ms and P 0 is the background spiking probability with P 0 (excitatory) = 10 −6 and P 0 (inhibitory) = 0 [19]. At each time step, each neuron spikes with probability where Θ is the Heaviside function. Then S is updated for the next time step. If a neuron spikes we reset R S to −2 (excitatory neurons) or −20 (inhibitory neurons) [19].

Neuronal avalanches
Given that some balance of r E :r I produce networks that are continuously active, a threshold is necessary to define the start and end of an avalanche. Poil et. al. [19] proposed a threshold defined as 50% of the median spike activity of the network. Here we defined the threshold θ as: where Γ > 0 is a number andm the median activity (so the original CROS model used Γ = 0.5). A neuronal avalanche starts and ends when the fluctuation of the summed activity of the network A(t) � ∑ i S i (t) crosses θ. If an avalanche starts at time t i and ends at time t f , its duration is T = t f − t i (Fig 2A). We studied two definitions for the size of an avalanche: i) s g ¼ P t f t i AðtÞ, i.e. the total area under the A(t) curve, as originally proposed by Poil et al. [19]; and ii) s y ¼ P t f t i ðAðtÞ À yÞ, i.e. the area above the threshold, as recently proposed by Del Papa et al. [27].
Importantly, we note that the avalanche definition of the CROS model differs from the one normally employed in the analysis of experimental data. In LFP [4], MEG [21][22][23] and fMRI studies [24], usually the activity at each site is thresholded to define a local "binary" event. The whole time series of events is then parsed with a time bin and an avalanche is defined as the spatio-temporal sequence of events between silent bins. When we compare model results with those obtained experimentally, this difference in the avalanche definitions should be taken into account. κ index. We used the κ index proposed by Shew et al [28] to obtain a first approximate assessment of how close an observed (simulated) distribution P obs is from a known theoretical distribution P th . We focused on power-law distributions of avalanche size and duration (both referred to simply as x in what follows). If we restrict the analysis of a power law distribution P th (x) � Cx −μ to an interval [x min , x max ], in the continuum limit one can calculate both the nor- R h x min P obs ðxÞ dx be the observed cumulative function. The κ index, which quantifies the difference between an observed and a theoretical distribution, is defined as where h j are b = 10 observables logarithmically spaced between x min and x max . With this definition, one obtains κ ' 1 for an observed distribution which is close to the theoretical power law, κ < 1 for a subcritical (e.g. exponential-like) distribution and κ > 1 for a supercritical distribution (say, with a characteristic bump at larger values of x). In our case μ stands for τ and τ t , critical exponents from the Direct Percolation (DP) universality class for size and duration of an avalanche, respectively. These exponents depend on the dimension of the underlying model. As mentioned previously, τ = 3/2 and τ t = 2 stand for DP models at or above the upper critical dimension (d � d c = 4), i.e. are mean-field exponents. Since the CROS model studied here is in principle two-dimensional, we also employ τ = 1.268 and τ t = 1.450, which are the critical exponents for two-dimensional DP [14]. We denote the indices to quantify the scale invariance of avalanche size s g or s θ and avalanche duration T as κ g , κ θ and κ T , respectively. Each of these three indices come in two variants, one employing mean-field (MF) exponents for the theoretical distributions and another employing their twodimensional (2D) counterparts. We present heat maps in parameter space of the different κ indices, averaged over 5 realizations of the disorder.
Maximum likelihood estimator. Estimating the quality of a power-law fit with the κ index has the drawback that the exponent μ must be known a priori. In order to loosen this constraint (while complying with more firmly grounded criteria for fitting power-law distributions [29]), we employ the maximum likelihood estimator (MLE) as made available in the powerlaw python package [30]. The exponents for size and duration of avalanches were estimated from a power law with an exponential cutoff, P(x) � x −μ exp(−x/x 0 ), with μ and x 0 as free parameters and fixing minimum values s min = 10 and T min = 4, respectively.
The package supports different probability distributions, which allowed us to use the loglikelihood ratio test [31] to compare the goodness of the fit of the exponentially truncated power-law with that of other possible distributions. We made comparisons with pure powerlaw, lognormal and exponential distributions. For our data, the exponentially truncated power-law was always a better fit than those three other distributions, for size as well as for duration of neuronal avalanches.
Scaling relation between avalanche size and duration. At criticality, one expects a power-law relation between the average size hsi of an avalanche of a given duration T: where σ, ν and z are standard critical exponents respectively governing the cutoff in the size distribution, the typical length of the largest avalanche and the duration of an avalanche as a function of its spatial extent [32]. Furthermore, at criticality one expects avalanche shape collapse, i.e. that the mean temporal profiles of avalanches can be described by a single function by means of scaling [25]. Let s(t, T) be the number of firings at time t of an avalanche with duration T. Then at criticality one should observe where F is a scaling function and 1/(σνz) is the same combination of critical exponents as in Eq 6. The MATLAB package NCC [33] was used to check whether shape collapse occurs in the model, at the same time yielding an independent estimate of 1/(σνz). For the aforementioned analysis only durations longer than 10 ms with at least 10 samples were used. Finally, at criticality the exponent in Eq 6 is related to the exponents ruling the distributions of avalanche size and duration [25,32,34]: These scaling relations were observed in in vitro and ex vivo experiments [25,35]. We probe the model by independently testing whether Eqs 6, 7 and 8 hold near its transition line.

Order parameter, power spectrum and DFA analysis
The frequency spectrum was obtained via the Fast-Fourier Transform (FFT), normalized by the network size (L × L) and the total time of simulation. The spectra were subsequently smoothened by averaging over nonoverlapping windows of 500 frequency steps (each step with width of 9 × 10 −4 Hz, resulting in one smoothened point every �0.45 Hz). The order parameter φ was defined as the ratio between the area of the power peak ϕ p divided by the total area ϕ p + ϕ u (Fig 2B): where ϕ u is the area under the power peak. φ is only different from zero if we can detect a peak in power-spectrum. Operationally, once we have the peak frequency f p , we detected the local minimum f − below it (between 4-18 Hz), which yields Δf � f p − f − . Using a symmetrical interval around f p , we calculated an upper frequency f + � f p + Δf. Given these three points, we fitted a line between f − and f + : the area above the line is ϕ p and the area under it is ϕ u . DFA analysis was performed in two different ways: i) on the raw time series A(t) and ii) the band-filtered time series A(t). For the former we followed the standard procedures described in Refs. [8,36,37]. For the latter, we followed the procedures described in Refs. [10,20]. Briefly, the occurrence of a self-similarity exponent α in the time series is a reliable indicator of long-range time correlation if 0.5 < α < 1, whereas α � 0.5 results from uncorrelated fluctuations. If the power spectrum decays as S(f) � 1/f v , then v = 2α − 1 [38]. So, when α = 1 we have the presence of 1/f noise. We estimated α from a fixed time range between 4 s and 14 s. The DFA exponents corresponding to the raw and band-filtered data are denoted by α raw and α, respectively.
Implementation of band-filtering relied on a Butterworth band-pass filter of 5th order over the frequencies 8-16 Hz. In order to extract the amplitude envelope we took the absolute value of the Hilbert transform. The preprocessing of the data was made in python through the SciPy software using the functions scipy.signal.butter and scipy.signal.hilbert for the Butterworth filter and the Hilbert transform, respectively.

Experimental M/EEG recordings
To compare the model results with those of experimental M/EEG recordings, the values of the exponents from the plots of Refs. [22,23] were extracted with the software WebPlotDigitizer (https://automeris.io/WebPlotDigitizer/).

System-size and threshold dependence
In the CROS model the definition of an avalanche depends on the threshold parameter θ (Fig  2A), whereas other signatures of criticality (DFA exponents, power spectra and order parameter) do not (Fig 2B). We therefore initially investigate these two types of markers of criticality separately, starting by those which do not depend on θ.
We start by addressing the dependence of the results on the system size L. The original results of Poil et al were obtained with L = 50 and ℓ = 7 [19]. As such, it is not clear whether the ratio ℓ/L is sufficiently small for the two-dimensional nature of model to be reflected in the results. We refer to the definition of topological dimension D which, as concisely summarized by Moretti and Muñoz, governs "how the number of neighbours of any given node grows when moving 1, 2, 3, . . ., r steps away from it: N r � r D , for large values of r" [39]. In a previous work on a similar topology, with large two-dimensional structure and disorder at small scale within an interaction range ℓ, results at the scale ℓ were governed by the random-like (disordered) local structure, with mean-field like behavior, whereas exponents compatible with D = 2 were observed at larger scales [40]. Bidimensionality makes sense if distances r are allowed to become much larger than any local interacting range.
To probe the robustness of the results, we therefore kept ℓ = 7 while increasing system size up to L = 300 (S1 Fig). In all cases, we obtained the following scenario, which coincides with the original results of Poil et al [19] and which we summarize for the largest system size in Fig  3. Keeping r I fixed and increasing r E , one observes not only an increase in the average activity hAi, but eventually also the emergence of alpha-band oscillations in A(t) (Fig 3A). These can be characterized by power spectra, in which a peak emerges (Fig 3C) at a transition region (line) in the (r E , r I ) parameter plane (Fig 3B). Alternatively, the order parameter φ can be used to estimate the transition line (Fig 3D), where its fluctuations are also maximal (Fig 3E). For each of these time series, we have assessed LRTCs in two different ways (see Methods). First, by calculating the DFA exponent α raw of the raw time series A(t) (Fig 4A and 4B), similarly to what has been previously done with spiking data [11,41]. Second, by calculating the DFA exponent α of the band-filtered time series (Fig 4C and 4D, see Methods), similarly to the procedure usually applied to M/EEG data [9,10,20,22,23].
Interestingly, the transition region is also characterized by nontrivial (larger than 0.5) DFA exponents (Fig 4E and 4F), as originally reported [19]. However, results for L = 300 suggest an interplay between system size and band-filtering. For L = 50, the effects of band-filtering are mild: both α raw and α reach values close to one at the transition region (Fig 4A, 4C and 4E). For L = 300, on the other hand, the values of band-filtered α are smaller than for L = 50 ( Fig  4C and 4D), while α raw is weakly affected (Fig 4B). Importantly, even for L = 300 and bandfiltered data, α remains above 0.5 in the transition region, thus signaling long-range time correlations.
Therefore, the above threshold-independent markers of criticality (power spectrum, order parameter and DFA exponents) seem to be robust to an increase in system size, suggesting a phase transition at which the system has long-range temporal correlations while oscillations start to emerge. In what follows, all results were obtained for networks with L = 300.
Turning our attention now to avalanches, we investigate the dependence of the results on the threshold y ¼ Gm (see Eq 4). The κ index is a useful tool to summarize how close the size and duration distributions are to power laws in parameter space (as long as one is confident about which exponent to expect, as we will discuss below). We used κ g and κ θ for the original (s g ) and above-threshold-only (s θ ) definitions of avalanche size, respectively (see Methods), with κ T characterizing the distributions of avalanche durations (T). For each of these three, we imposed either the usual mean-field (τ = 3/2, τ t = 2) or the 2D (τ = 1.268, τ t = 1.450) exponents of the directed percolation universality class, in a total of six variants of κ. Exploring their heat map in parameter space for Γ = 0.5, 0.75 and 1 (S2 Fig), we obtain critical regions in parameter space where κ ' 1, pointing to power-law distributions of sizes and durations. Since the results for s θ were more robust than for s g , the former will be used throughout. Except for the largest value of Γ, the phase diagrams also proved to be rather insensitive to this parameter (S2 Fig), so we will use Γ = 0.5 (as originally proposed [19]) unless otherwise stated.

Power-law distributed avalanches: Exponents from 2D or mean-field directed percolation?
With a larger system size (L = 300) and local interactions (ℓ = 7), we turn to the question whether the avalanche distributions produced by the model are compatible with the universality class of two-dimensional directed percolation (2D-DP). Starting with size distributions, if one imposes the 2D-DP exponent τ = 1.268 and calculates the corresponding κ θ , a critical region with κ θ ' 1 is indeed found (Fig 5A). However, if the mean-field directed percolation (MF-DP) exponent τ = 3/2 is used instead, the heat map of the corresponding κ θ index also displays a critical region (Fig 5B) near the 2D-DP one. In both cases, a close inspection of a few representative points in parameter space (Fig 5C and 5D) shows clear subcritical distributions (consistent with κ θ < 1) as well as supercritical distributions with a mild bump (κ θ > 1; moving just a bit further into the supercritical region, one obtains essentially a few giant avalanches-often just one). More importantly, fitting the exponent of a power-law with exponential cutoff (with the MLE method) at the putative critical region, one obtains an exponent reasonably close to the value for 2D-DP at one point in parameter space ( Fig 5C) and another exponent close to that of MF-DP at a different point (Fig 5D). Both power laws hold for at least three decades at the points with κ θ � 1. A very similar scenario is observed when the distributions of avalanche duration P(T) are studied. Heat maps for the κ T indices for 2D-DP and MF-DP exponents (Fig 6A and 6B, respectively) both show κ T ' 1 in putative critical regions, with subcritical and supercritical behaviors off these regions. The MLE-fitted exponents are again close to the theoretical 2D-DP and MF-DP values, perhaps with a better fit in the first case in comparison with the second (approximately three decades in Fig 6C versus two in Fig 6D).
Given the current absence of a proper theoretical connection between a transition to collective oscillations and the DP universality class, the above presented results suggest that the κ index is not the most appropriate tool for clarifying this issue, since by construction it relies on a priori knowledge of the distribution exponent. Henceforth we therefore relax this constraint, taking a more agnostic approach towards the values of the exponents and letting the MLE method determine them.

Maximum-likelihood estimator and continuously varying exponents
If one fixes r I = 0.6 and measures the order parameter φ as the excitatory connectivity r E changes, the characteristic plot of a second-order phase transition emerges (Fig 7A). Within a transition region (shaded area of Fig 7A), φ departs continuously from zero and, consistently, the DFA exponent peaks (Fig 7A corresponds to a cross section of Figs 3D, 4B and 4D). This region is defined by the fact that within it the distributions of avalanche sizes are better fitted by power laws with exponential cutoffs than by exponentials or lognormals, according to the loglikelihood test (see Methods). Similar results are observed for different system sizes (left column of S3 Fig). We also note that the constant Poisson input P 0 driving the excitatory neurons (Eq 2) is small (see Methods), so in the subcritical regime the average firing rate F (proportional to the density hρi of active sites), while strictly nonzero, is also very small. In the transition region, the rise of φ is accompanied by the rise of F, which for all practical and numerical purposes could also be used as an order parameter in this case (S4 Fig). This seems to be a particularity of the CROS model for this set of parameters, other models usually show different critical points for the onsets of self-sustained activity and oscillations [42,43].  Zooming in this transition region, the MLE method yields exponents for the distributions of size (Fig 7B) and duration (Fig 7C) which vary continuously as parameter space is traversed. Note that these are averages over different realizations of the disorder (see Methods). Eventually, both 2D-DP and MF-DP values are crossed for the truncated power laws which best fit the size and duration distributions.
At this point, it is unclear whether either of these two theoretical conjectures should be chosen as a better description of criticality in the model. We therefore applied a simple test, asking whether there is some combination of parameters for which the values of both τ and τ t come close to those of either of the two conjectured universality classes. From the functions hτi(r E ) and hτ t i(r E ) (Fig 7B and 7C), one can parametrically plot hτi vs hτ t i, either with (Fig 7E) or without (Fig 7D) averaging over the disorder. In the plane (τ, τ t ), simulation results eventually come close to the 2D-DP exponents, but not to the MF-DP exponents. This in principle seems consistent with the fact that the model is two-dimensional, especially since care has been taken to increase the system size.
But how robust are these results? We should now revisit their dependence on the threshold parameter Γ, a crucial ingredient of the model that directly impacts the very definition of an avalanche (Methods, Eq 4). While the κ index proved to be fairly robust against variations in Γ (S2 Fig), its usefulness was shown to be rather limited, because its calculation requires a priori knowledge of the power law exponent. So how does the overall picture of the much more informative Fig 7 change when the threshold is varied?
First, we observe that the absolute values of both avalanche exponents increase with increasing Γ (Fig 8A and 8B). This is interesting, because it relates to a similar result originally obtained by Beggs and Plenz [4]: by decreasing the time bin employed to slice their time series, they observed an increase in the absolute value of τ. In other words, by imposing stricter conditions (shorter silences), larger avalanches became less likely, in the sense that the power law describing the size distribution became steeper. Despite the differences in how avalanches are defined, the CROS model reproduces the same trend (also seen in other experimental and modelling setups [23,41,43,45]), to the extent that imposing larger thresholds is tantamount to decreasing the bin width.
Second, as the values of τ and τ t change with Γ, they keep a reasonably linear relationship with each other (Fig 8C and 8D), which is consistent with the right-hand side of Eq 8 being approximately constant. Moreover, and most importantly, this linear relationship depends on Γ. For increasing Γ, the linear spread of exponents in the (τ, τ t ) plane is gradually displaced, eventually departing from the 2D-DP values and approaching the ones for MF-DP. Therefore, while for Γ ' 0.5 it is possible to simultaneously find both 2D-DP exponents (τ ' 1.268 and τ t ' 1.45), for Γ ' 0.85 it is possible to find both MF-DP exponents (τ = 1.5 and τ t = 2) at the transition region of the model. These results remain qualitatively the same for different system sizes (middle column of S3 Fig).

More stringent tests of criticality
From a theoretical point of view, the presence of power laws with exponents that vary continuously (Figs 7 and 8) poses a challenge to the CROS model: where in parameter space is the phase transition (assuming there is one)? It has long been argued that power law distributions, on their own, are not a sufficient signature of a critical point, and indeed several other approaches have been proposed to substantiate claims of criticality [24][25][26]. Here we ask whether the exponents satisfy three different scaling relations [25].
First, at criticality one expects the average avalanche size hsi of a given duration T to scale as hsiðTÞ � T 1 snz (Eq 6). We have confirmed that this holds within the whole transition region of the CROS model for different system sizes and threshold values (right column of Fig 9).
Second, if s(t, T) is the number of firings at time t of an avalanche of duration T, then at criticality one should be able to observe shape collapse, that is, the rescaling sðt; TÞ � T 1 snz À 1 F ðt=TÞ (Eq 7) of the temporal profiles of avalanches of different durations, where F is a scaling function. At the transition region, the avalanches of the CROS model do satisfy shape collapse (left and middle column of Fig 9), with an exponent that is in good agreement with the one independently estimated by the scaling relation of Eq 6 (compare the exponent values in the middle and right columns of Fig 9).
The fact that the two scaling relations above hold within the whole transition region of the CROS model does not allow one to pinpoint where the critical point is and what its critical exponents are. In fact, Eq 6 was reported to hold for experimental data even away from criticality [25]. Besides, Touboul and Destexhe have recently argued that models which are not critical can satisfy both these scaling relations [26]. There is, however, a third scaling relation connecting the exponent in Eq 6 and those governing the distributions of avalanche size and duration: ðtÀ 1Þ (Eq 8). Satisfying this scaling relation is regarded as a much more stringent criterion for criticality, and we are unaware of non-critical models where it holds [26].  Fig 7). The color code in (C) and (D) represents different values of Γ, as described in (C). Orange points are experimental M/EEG results extracted from Palva et al. [22]. In all figures, error bars represent the standard deviation over 5 runs of the CROS model. https://doi.org/10.1371/journal.pcbi.1006924.g008 We proceed to test Eq 8 in the CROS model by noting that its left-hand and right-hand sides can be independently evaluated, and then compared for consistency. While the exponent in the left-hand side was estimated from fitting Eq 6, the right-hand side employed the avalanche exponents obtained from the maximum-likelihood estimator. Results are summarized in Fig 10. For fixed r I , as the parameter r E is varied within the transition region (shaded area of Fig 7A), the left-hand side stays below the right-hand side, failing to satisfy the scaling relation (Fig 10). Note that this failure is independent of a requirement to fit either conjectured universality class (2D-DP or MF-DP): the left-and right-hand sides of Eq 8 do not cross anywhere for the CROS model. Changing the threshold parameter Γ does not fix this issue. In fact, note that the difference between the left-and right-hand sides reaches a minimum for lower thresholds (Fig 10). While lower thresholds yield avalanche exponents closer to those of 2D-DP ( Fig  8C and 8D

Comparison with M/EEG experimental results
Despite the difficulties in reconciling our results with a proper phase transition (DP-like or otherwise), the scale-invariant behavior of the CROS model in its transition region is very rich and lends itself to comparison with experimental data. In particular, we wanted to probe whether the spread of exponents yielded by model is able to reproduce features of phenomena observed experimentally. Two main sources impact the variability of experimental results: first, the dependence on fitting parameters. As we discussed previously, the CROS model captures, for instance, the increase in the slope of the power laws as the threshold (implicit in the avalanche definition) is increased (Fig 8A and 8B). But even if fitting parameters are fixed, a second factor comes into play, namely, inter-subject variability. In the following, we revisit human magneto-and electroencephalographic (M/EEG) results that show a spread of exponents across subjects, and check how the CROS model performs in reproducing them. In this comparison, it is important to keep in mind some caveats: for instance, the data analysis and the model employ different avalanche definitions (see Methods); also, the topology of the CROS model is quite simplified, neglecting known features of cortical interactions, such as long-range connections [46].
In M/EEG recordings of humans performing a threshold-stimulus detection task, Palva et al have found that source-reconstructed data exhibit robust power-law long-range temporal correlations (i.e. well characterized DFA exponents) as well as scale-invariant avalanches (obtained by converting the recordings into a binary point process via a threshold at three standard deviations [22]). Results support the hypothesis of the human brain working near a critical point. But the avalanche and DFA exponents varied across subjects. When replotted in the (τ, τ t ) plane, their data spreads linearly (Fig 8C), a fact originally unreported [22] which, again, is consistent with a constant right-hand side of Eq 8. Moreover, the slope of the spread is reproduced by the CROS model, which achieves a reasonable quantitative agreement with the data for larger values of the threshold parameter (Fig 8C).
Palva et al also showed that avalanche and DFA exponents were not independent, but rather negatively correlated, with similar results observed for recordings of subjects at rest, as well as for the dependence of the DFA exponent of behavioral time series (hits and misses) [22]. To appreciate the implications, it is important to emphasize the difference in time scales: while avalanches can last at most a couple of hundred milliseconds, long-range time correlations in the M/EEG data were observed for tens of minutes [22].
This very interesting experimental result, which connects signatures of criticality at very different time scales, is not well reproduced by the CROS model. If one looks at single runs instead of averages (assuming that the result for a given experimental subject would be modeled by a single realization of the model disorder), simulations within the transition region do indeed give rise to raw-data DFA and avalanche exponents negatively correlated (lower plots of S5 Fig). But a fundamental discrepancy between model and experiments remains, because the DFA procedure was applied to the M/EEG data after band-filtering [22], whereas negative correlations in the CROS model are dominant for unfiltered data only (S5 Fig).
In a subsequent MEG study by the same group, however, Zhigalov et al dedicated special attention to the dependence of the exponent τ (governing avalanche size distributions) on two fitting parameters: the threshold T 0 used to define a binary event at each site, and the time bin Δt used to parse the binary events and define an avalanche [23]. For each point in the (T 0 , Δt) plane, a different distribution of avalanche sizes was obtained. In particular, when they analyzed the region in the (T 0 , Δt) plane where truncated power laws (with an exponential cutoff) were the best fit to the size distributions, they obtained a positive correlation (r = 0.80) between the band-filtered DFA exponents and τ (Fig 11). Moreover, they noticed that in this region of parameter space the size exponent was close to its MF-DP value, τ = 3/2 [23].
To compare the CROS model results with those of Zhigalov et al, we chose Γ = 0.85, a threshold value which yields exponents close to the MF-DP universality class (Fig 8C). We explored the values of band-filtered α and τ (also fit by power laws with exponential cutoff, see Methods) within the transition region. The CROS model results for single runs in the (α, τ) plane show not only positive correlations (r = 0.90), but also a good quantitative agreement with the values of the exponents of the MEG results (Fig 11). The positive correlation is robust with respect to a different choice of Γ, as long as the system size is large (S5 Fig).

Conclusions
Given that the connection between critical phenomena and scale-invariant statistics has long been established, power-law distributions of avalanches is a reasonable first signature to look for if one believes a neuronal system is critical. It is important to keep in mind, however, that power-law distributions, on their own, are generally not a sufficient signature of criticality.
Moreover, even if the analysis is restricted to power-law distributions, it is important to remind that critical exponents depend on many factors, such as dimensionality, disorder and the nature of the phase transition. In particular, the avalanche size exponent τ ' 1.5 originally observed by Beggs and Plenz [4] coincides with any model in the MF-DP universality class. But other experimental setups have revealed different exponents [25,35,47]. Likewise, models with different topologies or ingredients in their dynamics might give rise to exponents which differ from those of MF-DP.
Models with disorder and without an analytical solution, such as the one we studied here, can be particularly challenging. The CROS model is two-dimensional, with disorder, and the onset of its oscillations in principle seems physically different from the absorbing-active phase transition of the DP universality class. One of the main results of this manuscript is that avalanche and DFA exponents vary continuously within a transition region of parameter space where oscillations emerge (Figs 7 and 8). Therefore, if one insists on the mean-field DP values and looks for the τ ' 1.5 and τ t ' 2 exponents, one may find one or the other, as we did by using the κ indexes as markers of criticality (Figs 5B and 6B). These, however, proved to be insufficient signatures. For instance, fixing the threshold parameter at Γ = 0.5, the MF-DP exponents could not be obtained simultaneously (Fig 7D and 7E). What we did obtain simultaneously for that threshold value were the size and duration exponents of 2D-DP (Fig 7D and  7E), a result that might seem reassuring given that the model is indeed two-dimensional and assuming that a connection between directed percolation and the onset of oscillations can be established at all. But even that shred of consistency disappears if the threshold is increased: for Γ ' 0.85 we did obtain MF-DP exponents simultaneously. In fact, as shown in Fig 8C, it is possible to cover a significant fraction of the (τ, τ t ) plane by varying two parameters, so it would be no surprise if the model could reproduce exponents of other universality classes as well. Recently, di Santo et al have investigated a different model with a similar claim, namely, that power-law distributed avalanches occur at the edge of synchronization [43]. They also simulated a two-dimensional model, yet exhibited MF-DP avalanche exponents. It would be interesting to check what the strategies presented here would reveal if applied to their model.
In any case, a theoretical understanding of whether DP can be an "effective theory" for the phase transition leading to oscillations is still missing. Consider, for instance, a recent study by Coleman et al that has numerically investigated avalanches in a variant of the globally coupled version of the paradigmatic Kuramoto model (MF-K) [44], a minimal model to understand the onset of collective oscillations. Tuning the large-N model to the slightly subcritical regime, the authors define avalanches as excursions of the Kuramoto order parameter between consecutive zero crossings. Through scaling analysis, they arrive at avalanche exponents τ = 1.15 ± 0.1 and τ t = 1.28 ± 0.05 [44]. Note that, on the one hand, these values are far from those of the MF-DP universality class, despite the mean-field nature of the model. Interestingly, on the other hand, their exponents are compatible with those we observed in the transition region of the two-dimensional CROS model for Γ = 0.5 (Fig 7D and 7E). However, the closest approach of the left-and right-hand sides of Eq 8 for the CROS model occurs for a value of (τ t − 1)/(τ − 1) far from that of MF-K (Fig 10). In light of these results, therefore, it seems theoretically unwarranted to expect the "classical" MF-DP exponents (τ = 3/2 and τ t = 2) at the onset of collective oscillations.
Besides, even though the two-dimensional CROS model can exhibit exponents for avalanche size (τ) and duration (τ t ) which are compatible with DP, we did not obtain a fully satisfactory connection between the two scenarios. At a DP-like phase transition, avalanches are expected to satisfy scaling relations involving their size and duration (Eq 6), shape collapse (Eq 7) and consistency between the exponents involved (Eq 8). The first two scaling relations are often satisfied experimentally [25,35] and we found them to be valid within the transition region of the CROS model as well (Fig 9). However, they are also satisfied by models which are not critical [26]. The third scaling relation (Eq 8), on the other hand, is considered more stringent and hence a more decisive criterion for criticality [26]. By plotting is left-and right-hand sides as function of a coupling parameter, we expected the scaling relation to determine which, among the spread of exponents in the (τ, τ p ) plane, are the truly critical ones. However, we did not find a region in parameter space where Eq 8 held (Fig 10).
The fact that our results for the CROS model fail to satisfy Eq 8 should be interpreted with some caution. It may signal that the model does not have a true critical point, but rather undergoes some other transition phenomenon. If so, this would go in the direction of a number of other works which suggest alternatives to strict criticality as theoretical possibilities to explain the available experimental data. Some of these possibilities are, for instance, self-organized quasi-criticality (oscillations around a critical point [48]), pseudo self-organized criticality (a broad region of parameter space "close" to a critical point [49]), Griffiths phases (continuously varying exponents around a true critical point [39,50]), or even the Boltzmann chaos regime [26]. But there is also the possibility that the model would satisfy Eq 8 if one employed the regular definition of avalanche, with local thresholds to define binary events and a time bin to parse the series of events. Though computationally expensive, this would be a natural next step not only to better understand the model, but also to bring its methodology closer to that employed in the analysis of experimental data, thus improving the quality of the comparison. Along the same line, it would be worth investigating how the CROS model behaves if its synaptic connectivity incorporates more realistic ingredients, such as sparse long-range connections. As is well known, the small-world property is often sufficient to throw models into the meanfield version of their universality class [51]. On the one hand, this might help pinpoint a true critical point in the model. On the other hand, it is not apparent a priori whether the continuously varying exponents of the CROS model that we have obtained in its original formulation would remain.
Finally, we have shown that the CROS model is a very promising starting point for modeling some scale-invariant properties of M/EEG data. Despite the differences in avalanche definition between the model and the data analysis, the model quantitatively reproduces the linear relation between the avalanche exponents obtained by Palva et al [22] (Fig 8C), as well as the positive correlation between τ and the DFA exponent α obtained by Zhigalov et al [23] (Fig  11). In light of the results presented here, future lines of investigation naturally emerge on both the modelling and experimental fronts. For instance, since the agreement with the data occurs near a transition region of the model, it would be worth searching for dynamical mechanisms which could lead the system to self-organize around that narrow region of parameter space. On the experimental side, the more stringent tests of criticality that we have applied to the CROS model might be extended to the experimental data as well, raising the question whether it would be able to reveal, among a spread of exponents, the truly critical ones.