article: 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.
  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–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 (dc = 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 〈ρ〉 (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 〈ρ〉 departs continuously from zero as 〈ρ〉 ∼ (λ − λ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. 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 (rE and rI), 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 (rE, rI) 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 mean-field 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 activity m ˜. 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–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–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 maximum-likelihood 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.
  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 (rE) and inhibitory (rI) 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) = C e - r / a 0, where a0 = 1 is the nearest-neighbor distance, and C a normalization constant such that ∑neighbors P (r) equals rE or rI for presynaptic excitatory and inhibitory neurons, respectively. Notice that the networks have disorder, in the sense that, for a given point (rE, rI) 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 rE and rI were created and simulated for 220 time steps (dt) of 1 ms. Each connection had a weight (Wij) 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 Ii, which is also subject to an exponential synaptic decay: I i (t + d t) = (I i (t) + ∑ j W j i S j (t) ) (1 - d t τ I). (1) In Eq (1), S (t) ∈ {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 WEI = 0. 011, WEE = 0. 02 and WIE = WII = −2 and the time constant is τI = 9 ms [19]. Given the synaptic current Ii we compute the spiking rate variable RSi, which is also subject to an exponential decay: R S i (t + d t) = [ R S i (t) + I i (t) ] (1 - d t τ P) + P 0 d t τ P, (2) where τP (excitatory) = 9 ms, τP (inhibitory) = 12 ms and P0 is the background spiking probability with P0 (excitatory) = 10−6 and P0 (inhibitory) = 0 [19]. At each time step, each neuron spikes with probability P S = { R S i × Θ (R S i), if R S i < 1 1, if R S i ≥ 1 (3) where Θ is the Heaviside function. Then S is updated for the next time step. If a neuron spikes we reset RS to −2 (excitatory neurons) or −20 (inhibitory neurons) [19]. Given that some balance of rE: rI 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: θ = Γ × m ˜, (4) where Γ > 0 is a number and m ˜ 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 Si (t) crosses θ. If an avalanche starts at time ti and ends at time tf, its duration is T = tf − ti (Fig 2A). We studied two definitions for the size of an avalanche: i) s g = ∑ t i t f A (t), i. e. the total area under the A (t) curve, as originally proposed by Poil et al. [19]; and ii) s θ = ∑ t i t f (A (t) - θ), 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–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. 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): φ = ϕ p (ϕ u + ϕ p), (9) 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 fp, we detected the local minimum f− below it (between 4–18 Hz), which yields Δf ≡ fp − f−. Using a symmetrical interval around fp, we calculated an upper frequency f+ ≡ fp + Δ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/fv, 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. 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/).
  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: Nr ∼ rD, 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 rI fixed and increasing rE, one observes not only an increase in the average activity 〈A〉, 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 (rE, rI) 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 band-filtered 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 θ = Γ m ˜ (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 (sg) 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 sg, 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. 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. If one fixes rI = 0. 6 and measures the order parameter φ as the excitatory connectivity rE 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 P0 driving the excitatory neurons (Eq 2) is small (see Methods), so in the subcritical regime the average firing rate F (proportional to the density 〈ρ〉 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 〈τ〉 (rE) and 〈τt〉 (rE) (Fig 7B and 7C), one can parametrically plot 〈τ〉 vs 〈τt〉, 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). 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–26]. Here we ask whether the exponents satisfy three different scaling relations [25]. First, at criticality one expects the average avalanche size 〈s〉 of a given duration T to scale as 〈 s 〉 (T) ∼ T 1 σ ν z (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 σ ν z - 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: 1 σ ν z = (τ t - 1) (τ - 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]. 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 rI, as the parameter rE 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), the “near miss” in Fig 10 occurs for values of 1/ (σνz) which are far from those of 2D-DP. This scenario is robust with respect to changes in system size (right column of S3 Fig). 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′ 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′, Δt) plane, a different distribution of avalanche sizes was obtained. In particular, when they analyzed the region in the (T′, Δ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). 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 mean-field 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.
summary: Since the first experimental observation of scale-invariant neuronal avalanches, the idea that the brain could be operating near a phase transition has received much attention. But if the brain is critical, what is the phase transition? Experimentally, the sizes and durations of local field potential (LFP) avalanches were power-law distributed with exponents 3/2 and 2. Theoretically, these are the exponents of a critical branching process, a coincidence which has inspired many models with a common feature: a phase transition from a silent (absorbing) to an active phase. These models, however, struggle with another experimental result: brain activity can also exhibit long-range temporal correlations (LRTCs, another fingerprint of a critical system). An attempt to model both phenomena was put forward via the CROS (“CRitical OScillations”) model, in which neuronal avalanches and LRTCs were observed at a transition where alpha-band oscillations emerge. We show that the avalanche exponents of the CROS model do not necessarily agree with those of a critical branching process, but can vary near the transition region, just like the exponents governing LRTCs. Moreover, the spread of exponents observed in the model is in good agreement with experimental results from human MEG data.
keywords: action potentials
  medicine and health sciences
  classical mechanics
  fluid mechanics
  neural networks
  membrane potential
  social sciences
  electrophysiology
  neuroscience
  brain mapping
  network analysis
  neuroimaging
  color codes
  research and analysis methods
  computer and information sciences
  imaging techniques
  animal cells
  behavior
  fluid dynamics
  continuum mechanics
  physics
  cellular neuroscience
  psychology
  percolation
  cell biology
  data visualization
  physiology
  neurons
  biology and life sciences
  cellular types
  physical sciences
  magnetoencephalography
  neurophysiology
title: Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic M/EEG results