Hallmarks of criticality, such as power-laws and scale invariance, have been empiricaly found in cortical networks and it has been claimed that operating at criticality entails functional advantages, such as optimal computational capabilities, memory, and large dynamical ranges. An apparent paradox is that criticality requires a high degree of fine tuning to emerge and hence, self-organizing mechanisms need to be invoked. Here we show that, taking explicitly into account the complex hierarchical-modular architecture of cortical networks, the standard critical point picture is replaced by an extended critical-like region which corresponds -in the jargon of statistical mechanics -to a Griffiths phase. Using computational and analytical approaches, we find Griffiths phases in synthetic hierarchical networks and also in empirical neural networks such as those of C. elegans and the human connectome. Stretched critical regions, stemming from structural disorder, yield optimal functionality in a generic way, facilitating the task of self-organizing, adaptive and evolutionary mechanisms selecting for criticality. Empirical evidence that living systems can operate near critical points is flowering in contexts ranging from gene expression patterns [1] , to optimal cell growth [2] , bacterial clustering [3] , or flocks of birds [4] . In the context of neuroscience, synchronization patterns have been shown to exhibit broadband criticality [5] and critical avalanches of spontaneous neural activity have been consistently found both in vitro [6] [7] [8] and in vivo [9] suggesting that criticality -with its concomitant power laws and scale invariance -might play a relevant role in intact-brain dynamics [8, 10] . At variance with inanimate matter -for which the emergence of generic or self-organized criticality in sandpile models, type-II superconductors, or solar flares is relatively well understood, [11] [12] [13] -criticality in living systems can be conjectured to be the result of evolutionary or adaptive processes, which for reasons to be understood select for it. The criticality hypothesis [8, 10, 14] states that biological systems can perform the complex computations that they require to survive only by operating at the "edge of chaos" (or critical) i.e. at the borderline between the chaotic phase in which noise propagates unboundedly-thereby corrupting all information processing or storage-and an ordered phase in which perturbations readily fade away, hindering the ability to react and adapt [15, 16]. Critical dynamics provides a delicate trade-off between these two impractical tendencies, and has been argued to imply optimal transmission and storage of information [7, 8, 17], optimal computational capabilities [18], large network stability [16], maximal variety of memory repertoires [19] , and maximal sensitivity to stimuli [20] . Such a "delicate" balance occurs just at a singular or critical point, requiring a precise fine tuning. A very recent functional magnetic resonance imaging (fMRI) analysis of the brain at its resting state, reveals that the brain spends most of the time wandering around a broad region near a critical point, rather than just sitting at it [21] . This suggests that the region where cortical networks operate is not just a critical point, but a whole extended region around it. No satisfactory explanation of this finding has been proposed yet. Here, inspired by this empirical observation as well as by some recent findings in network theory and neuroscience [22] [23] [24] we scrutinize the dynamics of simple models of neural activity propagation when the structural architecture of brain networks is explicitly taken into account. Using a combination of analytical and computational tools, we show that the intrinsically disordered (hierarchical and modular) organization of such networks dramatically alters the dynamics by inducing the emergence -in the jargon of Statistical Mechanics -of a Griffiths phase (GP) [22, 25] . This phase, which stems from the presence of disorder (structural heterogeneity here) is characterized 2 by generic power-laws extending over broad regions. Furthermore, functional advantages usually ascribed to criticality, such as a huge sensitivity to stimuli, are reported to emerge generically all along the GP. Remarkably, not only do we find GPs in stylized models of brain architecture, but also in real neural networks such as those of the C. elegans and the human connectome. Our conclusion is that, as a consequence of the intrinsically-disordered architecture of brain networks, critical-like regions are extended from a singular point (a sharp "edge-of-chaos") to a broad or stretched region (the "antechamber of chaos"), much as evidenced in recent fMRI experiments. The existence of Griffiths phases severely facilitates the task of self-organizing, adaptive, or evolutive mechanisms seeking for critical-like attributes with all its alleged functional advantages. We claim that the unavoidable structural heterogeneity of cortical networks calls for a change of paradigm from the critical/edge-of-chaos to a new one, relying on extended critical-like Griffiths regions. Our work also raises a series of questions worthy of future pursuit. For instance, is there any connection between GPs and the empirically reported generic power-law decay of short-time memories? Do our results extend to other hierarchical architectures such as those encountered in metabolic or technological networks? I. RESULTS The cortex network has been the focus of attention in neuroanatomy for a long time, but only recently, the development of high-throughput methods has allowed the unveiling of its intricate architecture or connectome [26, 27] . Brain networks have been found to be structured in a hierarchical fashion across many scales [28] [29] [30] . Moduli exist at each hierarchical level: "cortical columns" arise at the lowest level, "cortical areas" at intermediate ones, and "brain regions" emerge at the systems level [28] [29] [30] [31] . Each modulus is characterized by having a much denser connectivity within it than with elements in other modules at the same layer. To be able to perform systematic analyses, we design synthetic hierarchical and modular networks (HMNs) with s hierarchical levels, N nodes/neurons and L links/synapses, whose structure can be tuned to mimic real networks. Different approaches can be employed to do so. We employ two different models based on a "bottom-top"approach; in the first, local fullyconnected moduli are constructed and then recursively grouped by establishing new inter-moduli links between them in either a stochastic way with a level-dependent probability p as sketched in Figure 1 (HMN-1) or in a deterministic way with a level-dependent number of connections (HMN-2). Similarly, "top-down" models can also be designed [32] . A way to encode key network structural information is the topological dimension, D, which measures how the number of neighbors of any given node grows when moving 1, 2, 3, ..., r steps  Disorder is well-known to radically affect the behavior of phase transitions (see [25] and refs. therein). In disordered systems, there exist regions characterized by parameter values which differ significantly from their corresponding system averages. Such rare regions can, for instance, induce the system to be locally ordered, even if globally it is in the disordered phase. In this way, in models of activity propagation activity can linger for long times in rare-regions even if the system is in the quiescent phase. In the particular case in which broadly different rare-regions exist -with broadly distinct sizes and time-scales -the overall system behavior, determined by the convolution of those different contributions, becomes anomalously slow (see MM). In contrast with standard critical points, systems with rare-region effects have an intermediate broad phase, separating order from disorder: a Griffiths phase (GP) with generic power-law behavior and other anomalous properties (e.g. [25] and MM). Remarkably, it has been very recently shown that structural heterogeneity can play in networked systems a role analogous to that of standard quenched disorder [22] . In particular, simple dynamical models exhibit GPs when running upon networks with a finite D; instead, in small-worlds, D = ∞, local neighborhoods are too large -quickly covering the whole network -as to be compatible with the very concept of rare (isolated) regions [22] . Therefore, the existence of finite a topological dimension D is an excellent indicator of eventual rare-region effects and Griffiths phases. To model the propagation of neuronal activity, we consider highly-simplified dynamical models running upon HMNs. More realistic models of neural dynamics with additional relevant layers of information could be considered; but we do not expect them to significantly affect our conclusions here. Our approach here consist in modeling activity propagation in a minimal way. It is also noteworthy that in some of the cases we study the network nodes are not neurons but coarse neural regions, for which effective models of activity propagation are expected to work well. 1 Our dynamics is as follows: every node (or "neuron") is endowed with a binary state variable, σ, representing either activity (σ = 1) or quiescence (σ = 0). Each active neuron is spontaneously deactivated at some rate µ (we fix µ = 1), while it propagates its activity to other directly connected 1 The usefulness of these minimalistic models in Neuroscience has been ascertained before e.g. [34] . 6 neurons at rate λ. We have considered two different dynamics: in the first one, (Model A) a synapsis between an active an a quiescent node is randomly selected at each time and proved for activation, while in the other variant (Model B) a neuron is selected and all its neighbors are proved for activation. Details of the computational implementation can be found in the SI. In general, depending on parameters these models can be either in an active phase -for which the density ρ of active nodes reaches a steady-state value ρ s > 0 in the large system-size and large-time limit -or in the inactive phase in which ρ falls ineluctably in the quiescent configuration (ρ s = 0). Separating these two regimes, at some value, λ c , there is a standard critical point where the system exhibits power-law behavior for quantities of interest, such as the time decay of a homogeneous initial activity density, ρ(t) ∼ t −θ , or distribution of sizes of avalanche triggered by an initially localized perturbation P (S) ∼ S −τ . θ and τ are critical exponents. This standard critical-point scenario (see Figure 2 , left panel) holds for regular lattices, Erdős-Renyi networks, and many other types of networks (without disorder or with disorder but with an infinite topological dimension). On the other hand, computer simulations of the different dynamical models running upon our complex HMN topologies with finite D reveal a radically different behavior (see Fig. 2 and MM). The power-law decay of the average density ρ(t) -specific to the critical point in pure systemsextends to a broad range of λ values. The existence of a broad interval with power-law decaying activity is supported by finite size scaling analyses reported in the SI. Likewise, as shown in Fig. 2 , avalanches of activity generated from a localized seed have powerlaw distributed sizes, with continuously varying exponents, in the same broad region. These features are fingerprints of a GP and have been confirmed to be robust against increasing system size (up to N = 2 20 ), using different types of HMN (HMN-1 with different values of α and p, HMN-2 models, all with finite D) and dynamical models, and performing measurements of different dynamic quantities (see SI). One of the main alleged advantages of operating at the edge-of-chaos is the strong enhancement of the system's ability to distinctly react to highly diverse stimuli. In the statistical mechanics jargon this stems from the divergence of the susceptibility at criticality [35]. How do systems with broad GPs respond to stimuli? To answer this question we measure the following two different quantities: The dynamic susceptibility, gauges the overall response to a continuous localized stimulus. It is defined as Σ(λ) = N (ρ f (λ) − ρ s (λ)), where ρ s (λ) is the stationary density in the absence of stimuli and ρ f (λ) is the steady-state density reached when one single node is constrained to remain active. As shown in Fig. 3 Σ becomes extremely large in the GP and, more importantly, it grows as a power-law of system-size, Σ(λ ≈ λ c , N ) ∼ N η , implying that there is an extended divergence (with λ-dependent continuously varying exponents) of the system's response across the GP. (ii) The dynamic range, ∆, introduced in this context in [20] (see SI), measures the range of perturbation intensities/frequencies for which the system reacts in distinct ways, being thus able to discriminate among them. We have computed ∆(λ) in the HMN-2 model (see Fig. 3 ) which clearly illustrates the presence of a broad region with huge dynamic ranges (rather than the standard situation for which large responses are sharply localized around criticality). Therefore, if critical-like dynamics is important to access to a broad dynamic range and to enhance system's sensitivity, then it becomes much more convenient to operate with hierarchicalmodular systems, where the edge of chaos -with extremely large responses and huge sensitivitybecomes a region rather than a singular point. To further confirm the existence of Griffiths phases (beyond direct computational simulations), here we present some analytical results. An important tool in the analysis of network dynamics is provided by spectral graph theory in which the network structure is encoded in some matrix and linear algebra techniques are exploited [36]. For instance, the dynamics of simple models as those in previous sections (e.g. Model-B) is governed by the largest (or principal) eigenvalue, Λ max , of the adjacency matrix, A i,j (with 1's as entries for existing connections and 0's elsewhere) which straightforwardly appears in a standard linear stability analysis (see [37] and SI). It is easy to show that (with very mild assumptions) the critical point is such that λ c Λ max = 1, signaling the limit of linear stability of perturbations on the quiescent state. However, it has been recently shown that this result may not hold for certain networks, for which the largest eigenvalue has an associated localized eigenvector (i.e. with only a few non-vanishing components) 2 [37]. In such a case, linear instabilities for λ slightly larger than 1/Λ max lead to activity localized around only few nodes (where the corresponding eigenvector is localized), not pervading the network and not leading to active) the overall response increases (resp. decreases) for larger systems as marked by the red (resp. blue) arrows. In the GP, Σ grows as a λ-dependent power-law of N . Right axes: dynamic range ∆ (green squares) for the same case as above, N = 2 14 (dynamical protocol iv; see MM). Instead of the usual symmetric cusp-singularity at the critical point, there is a whole region of extremely large ∆ with a strong asymmetry around the transition point, as corresponds to the existence of a GP. a true active state, implying that the critical point is shifted to a larger value of λ. Instead these regions of localized activity resemble very much rare-regions in GPs: activity lingers around them until eventually a large fluctuation kills them. Inspired by this novel idea, we performed a spectral analysis of our HMNs (e.g. HMN-2 nets with α = 1), with the result that -for finite D networks -not only the largest eigenvalue Λ max corresponds to a localized eigenvector, but a whole range of eigenvalues below Λ max (even hundreds of them) share this feature, as can be quantitatively confirmed (see SI). In particular, the principal eigenvector is heavily peaked around a cluster of neighboring nodes. We have conjectured and verified numerically that the clusters where the largest eigenvalues are localized correspond to the rare-regions, with above-average connectivity and where localized activity lingers for long time. Also, we numerically found λ c ≈ 0.41 > 1/Λ max ≈ 0.33, confirming the prediction above. The interval between these two values defines the GP. In addition, we also considered large network ensembles and computed the probability distribution of eigenvalues. We found a type of singularity, usually termed Lifshitz tail, typical of disordered systems (see SI). Lifshitz tails are known to be a 9 spectral fingerprint of Griffiths phases [38]. Therefore, the presence of both (i) localized eigenvectors and (ii) Lifshitz tails, confirms the existence of GPs in networks with complex heterogeneous architectures. Analyses of different nature have revealed that organisms from the primitive C. elegans [29, 30, 39, 40] (for which a full detailed map of its about 300 neurons has been constructed) to cats, macaque monkeys, or humans (for which large-scale connectivity maps are known [26] ) have a hierarchically organized neural network. Such structure is also shared by functional brain networks (e.g. from functional magnetic resonance imaging data) [26, 27, [41][42][43]. Do simple dynamical models (as those in previous sections) running upon the C. elegans neuralnetwork exhibit GPs? A word of caution is needed before delving into this question. GPs (as well as standard phase transitions and critical points) do strictly-speaking exist only in the (thermodynamic) N → ∞ limit, while for small sizes they become fuzzy. Given that the C. elegans network comprises only N < 300 nodes, the maximum size of possible rare regions and the associated power-laws are necessarily cut off at small sizes and short times. Nevertheless, as illustrated in (see Fig. 4 , upper inset), simulations of the dynamical models above (Model B) show a significant deviation from the typical standard critical-point scenario, highlighting an interval of values of λ, in which the avalanche size distributions move from subcritical to supercritical behavior in a rather continuous fashion, exhibiting a stretched (power-law like) decay in a whole range of λ values, although for a very limited interval of avalanche sizes. To mitigate finite-size effects, we have also considered the human connectome network, obtained by Sporns and collaborators using diffusion imaging techniques [26, 27] . It is a highly coarse-grained mapping (as opposed to the previous detailed map) of anatomical connections in the human brain, comprising N = 998 brain areas and the fiber tract densities between them, with a hierarchical organization [29, 30] (see SI). In this case, we find that avalanches are clearly distributed as power laws, with smaller finite-size effects, in a broad range of λ-values (see Fig. 4 ). Actually, truncated power-laws of the form P (S) ∼ S −τ e −S/ξ -with λ-dependent values of τ -provide highly reliable fits of the size distributions, P (S), according to the Kolmogorov-Smirnoff criterion (see SI). Here it becomes apparent that the dynamics does not lead to a traditional critical-point scenario but, instead, it supports the picture of a broad critical-like region. This strongly suggests that if it were feasible to run dynamical models upon the actual human brain network (with about 10 More remarkably, we have also provided evidence that GPs appear in actual real networks such as those of the human connectome and C. elegans, even if owing to the limited available networksizes the associated power-laws are necessarily cut-off. Still, the best fits to data are provided by continuously varying power-laws, truncated by finite-size effects. It is noteworthy that disorder needs to be present at very different scales for GPs to emerge (otherwise rare regions cannot be arbitrarily large); this is why a hierarchy of levels is required. Plain modular networks -without the broad distribution of cluster sizes, characteristic of hierarchical structures -are not able to support GPs. We emphasize that GPs in our HMNs are induced merely by the existence of structural disorder: no additional form of neuronal or synaptic heterogeneity has been considered here; adding further heterogeneity would only enhance rare-region effects and hence GPs. For instance, we have also found GPs for the human connectome, when taking into account the relative weight of structural connections. It should be also noticed that all networks in our study are undirected in the sense that links are symmetrical, while synapses in cortical networks are not. Again, this does not jeopardize our conclusions: directness in the connections only strengthens isolation and hence rare-region effects and GPs. The existence of a GP with its concomitant generic power-law scaling -not requiring a delicate parameter fine tuning -provides a more robust and solid basis to justify the ubiquitous presence of scale-free behavior in neural data, from EEG or MEG records to neural avalanches. More in particular, it gives the key to understand why broad regions around criticality are observed in fMRI experiments of the brain resting state [21] . As we have shown, the system's response is extremely large (and diverges upon increasing the 12 network size) in a broad region of parameter space. This strongly facilitates the task of mechanisms selecting for the alleged virtues of scale-invariance and strongly suggests that a new paradigm is needed: a theory of self-organization/evolution/adaptation to the broad region separating order from chaos. In particular, GPs combined with standard mechanisms for self-organization [46][47][48] are expected to account for the empirically found dispersion around criticality [21] . This also invites the question of whether intrinsically disordered HMNs do indeed generically optimize transmission and storage of information, improve computational capabilities [18], or significantly enhance large network stability [16], without the need to invoke precise criticality. Usually, brain networks are claimed to be small worlds; instead we have shown that large-world architectures are essential for GPs to emerge. A solution to this conundrum was provided by Gallos et al. [42], who found that (functional) brain networks consist of a myriad of densely connected local moduli which altogether form a large world structure and, therefore, are far from being smallworld; however, incorporating "weak ties" to the network converts it into a small world preserving an underlying back-bone of well-defined modules (for this, it is essential that networks have a finite Hausdorff or fractal dimension (see [42] and MM), confirming that finite dimensionality is a crucial feature). In this way, cortical networks achieve an optimal balance between local specialized processing and global integration through their hierarchical organization. On the other hand, weak links are not expected to significantly affect the existence of GPs. Therefore, from this perspective, (i) achieving such an optimal balance and (ii) operating with critical-like properties, can be seen as the two sides of the same coin. This stresses the need to develop new models for the co-evolution of structure and function in neural networks. It is noteworthy, that, a mechanism for robust working memory without synaptic plasticity has been put forward very recently [49]. It heavily relies on the existence of heterogeneous local clusters of densely inter-connected neurons where activity (memories) reverberates. Not surprisingly, it leads to power-law distributed fade away times, which has been claimed to be the correlate of "power-law forgetting" [50]. This is a nice illustration of Griffiths phases at work. Given that disorder is an intrinsic and unavoidable feature of neural systems and that neural network architectures are hierarchical, Griffiths phases are expected to play a relevant role in many dynamical aspects and, hence, they should become a relevant concept in Neuroscience as well as in other fields such as systems biology where HMNs play a key role [51]. We hope that our work contributes to this purpose, fostering further research. 13 Topological dimension: Strictly speaking HMN-1 and HMN-2 networks are finite dimensional only for p = 1/4 in which case the number of inter-moduli connection is stable across hierarchical levels. For p > 1/4 (p < 1/4) networks become more a more densely (rarefied) connected as the hierarchy depth (i.e. the network size) is increased. Deviations from p = 1/4 create fractal-like networks up to certain scale, being good approximations for finite-dimensional networks in finite size. In some works (e.g. [42]), the Hausdorff (fractal) dimension D f is computed for complex networks. We have verified numerically that D f ≈ D in all cases. Dynamical protocols: We employ four different dynamical protocols: (i) Decay of a homogeneous state: all nodes are initially active and the system is let to evolve, monitoring the density of active sites ρ(t) as a function of time. (ii) Spreading from a localized initial seed: an individual node is activated in an otherwise quiescent network. It produces an avalanche of activity, lasting until the system eventually falls back to the quiescent state. The avalanche size S is defined as the number of activation events that occur for the duration of the avalanche itself. The process is iterated and the avalanche size distribution P (S) is monitored. (iii) Identical to (ii), except that the seed is kept active throughout the simulation (continuous stimulus). (iv) Identical to (ii), except that the seed node is subsequently reactivated with probability p stimulus = 1 − exp(−r ∆t) for t > 0 (Poissonian stimulus of rate r). Rare regions and Griffiths phases: Quenched disorder strongly influences the physics of phase transitions [25] . For illustration, consider our Model A running upon a regular lattice but with a node-dependent quenched spreading rate, λ(x) (with average valueλ). If we consider a value ofλ in between, λ c (pure)>λ c (disordered) then the disordered system is in its inactive or quiescent phase, but the existence of disorder can (and does!) generate rare active regions characterized by significantly over-average (actually, values of λ above the critical point of the pure model). In these regions, initial activity can linger for very long periods, but still, as such rare regions have a finite size, they ineluctably end up falling into the inactive state. Considering a rare active region of size s, it decays to the quiescent state after a typical time τ (s) which grows exponentially with cluster size, i.e. τ t 0 exp(A(λ)s), where t 0 and A(λ) do not depend on s (Arrhenius law). On the other hand the distribution of s values is also an exponential (very large 14 regions are exponentially rare). Therefore, the overall activity density, ρ(t), decays as the following convolution integral ρ(t) ∼ dsP (s) s exp [−t/(t 0 e A(λ)s )], which evaluated in saddle-point approximation leads to ρ(t) ∼ t −θ with θ varying continuously with the disorder average value,λ). This generic power-laws are characteristic of Griffiths phases. This is just an example of a generic phenomenon, thoroughly studied in classical, quantum, and non-equilibrium disordered systems systems [25] .
Moretti et al. used PEV localization to analyze the brain network dynamics REF .
Griffiths phases and the stretching of criticality in brain networks
Nature Communications 4, 2521 (2013)
Physics Biology Medicine