Griffiths phases and localization in hierarchical modular networks
Abstract
We study variants of hierarchical modular network models suggested by Kaiser and Hilgetag [Front. in Neuroinform., 4 (2010) 8] to model functional brain connectivity, using extensive simulations and quenched meanfield theory (QMF), focusing on structures with a connection probability that decays exponentially with the level index. Such networks can be embedded in twodimensional Euclidean space. We explore the dynamic behavior of the contact process (CP) and threshold models on networks of this kind, including hierarchical trees. While in the smallworld networks originally proposed to model brain connectivity, the topological heterogeneities are not strong enough to induce deviations from meanfield behavior, we show that a Griffiths phase can emerge under reduced connection probabilities, approaching the percolation threshold. In this case the topological dimension of the networks is finite, and extended regions of bursty, powerlaw dynamics are observed. Localization in the steady state is also shown via QMF. We investigate the effects of link asymmetry and coupling disorder, and show that localization can occur even in smallworld networks with high connectivity in case of link disorder.
Introduction
In neuroscience, the criticality hypothesis asserts that the brain is in a critical state, at the boundary between sustained activity and an inactive regime. Theoretical and experimental studies show that critical systems exhibit optimal computational properties, suggesting why criticality may have been selected in the evolution of the nervous system [1]. Although criticality has been observed in cell cultures [2, 3], brain slices and anesthetized animals [4, 5], the debate regarding criticality in alert animals and humans continues [6, 7, 8]. Thus the criticality hypothesis remains controversial in brain science; for a review see [9].
Normally, for a system to be at criticality, certain control parameters need to be tuned precisely, raising the perennial question of how such tuning is achieved in the absence of outside intervention. The possibility of selftuning is well known in statistical physics; the paradigm of selforganized criticality (SOC) has been studied since the pioneering work of [10]. Simple homogeneous models such as the stochastic sandpile exhibit criticality with power laws both in statics and dynamics. This has been understood as the result of a control mechanism that forces the system to an absorbingstate phase transition [11]
Real nervous systems, however, are highly inhomogeneous, so that one must take into account the effects of heterogeneities. It is well known in statistical physics that disorder can cause rareregion (RR) effects [12] that smear the phase transitions. These effects can make a discontinuous transition continuous, or generate socalled Griffiths phases (GP) [13], in which criticallike powerlaw dynamics appears over an extended region around the critical point. In these regions, moreover, nonMarkovian, bursty behavior can emerge as a consequence of a diverging correlation time. The intercommunication times of the nodes (which possess no internal memory) follow a fattailed distribution [14]. Thus, heterogeneities widen the critical region and can contribute to the occurrence of power laws. This provides an alternative mechanism for criticallike behavior without fine tuning, although attaining the GP does require some rough tuning, of the sort that is not difficult to find in biological systems.
It was shown recently that the topological heterogeneity of the underlying networks can result in GPs in finitedimensional systems [15] and can be a reason for the working memory in the brain [16]. Although many networks exhibit the smallworld property and so have an infinite topological dimension, naturally occurring networks are always finite, exhibit cutoffs, and therefore GPs can be expected as a consequence of inhomogeneous topology [17].
Many real networks can be partitioned seamlessly into a collection of modules. Each module is expected to perform an identifiable task, separate from the function of others [18]. It is believed that the brain is organized hierarchically into modular networks across many scales, starting from cellular circuits and cortical columns via nuclei or cortical areas to largescale units such as visual or sensorymotor cortex. At each level, connections within modules are denser than between different modules [19, 20, 21, 22]. Although empirical data confirm this modular organization on some scales [23], the detailed organization of brain networks is not yet experimentally accessible.
Two particular kinds of hierarchical modular networks (HMN1,HMN2) model were proposed and investigated numerically and analytically in [24]. On largeworld HMNs, which imply a finite topological dimension , models of the spread of activity exhibit powerlaw dynamics and rare region effects. However, these power laws are systemsize dependent, so that true GP behavior has not yet been proven. These authors also simulated spreading on empirical brain networks, such as the human connectome and the nervous system of C. elegans. Available empirical networks are much smaller than the synthetic ones and deviations from power laws are clearly visible. Both anatomical connections [19] and the synchronization networks of cortical neurons [27] suggest smallworld topology [26]. The brain network modules of [24] are weakly coupled in such a way that these HMNs are near the percolation threshold, as in the case of the models introduced in [15]. Note, however that requiring the network to be near percolation again raises tuning and robustness issues. Having weaker ties would lead to fragmented networks, while stronger ties result in infinite and the absence of a GP. In the present work we do not assume such fine tuning: we maintain a high density of short edges, rendering the networks well connected. Another way of preserving the integrity of HMN networks with finite dimension involves the random treelike structures studied here.
To study synchronization [28], commonly expected in brain dynamics, the Kuramoto model [29] has been implemented [30] on the same networks as studied in [24]. In this case even weaker rareregion effects were found, resulting in stretched exponential dynamics. One of the main purposes of our study is to delineate conditions a GP. While the networks studied are of finite size, repeating the process on many network realizations and averaging over them, we clearly see convergence towards GP dynamics. We assume that multiple random network realizations may occur in the brain over time, due to reconfigurations of the synapses or as a consequence of weakly coupled submodules and changes in overall intensity of connections between different brain regions due to neuromodulation.
In addition to the dimension, other topological features, which have yet to be studied systematically, may also influence RR effects. The clustering coefficient, defined as
(1) 
where denotes the number of direct links interconnecting the nearest neighbors of node , exhibits different scaling behavior in modular networks than in unstructured networks. While in SF networks it decays as and in random networks as , in HMNs it is constant. Furthermore, while in SF networks the clustering is degreenumber independent, in HMNs it decays as , with in the range 0.75  1 [31]. This means that the higher the node degree, the smaller is its clustering coefficient, possibly reducing its infectiousness in an epidemic process. This suggests an enhancement of localization as in SF models with disassortative weighting schemes [32].
In this work we study spreading models defined on HMNs embedded in 2D space, and investigate how various inhomogeneities (topological or intrinsic) contribute to the occurrence of localization and GPs. We begin by considering purely topological heterogeneity; later we study cases with disorder in the connections as well. It is well known that intrinsic link disorder can appear as the consequence of asymmetry, connecting nodes in the cortex (see [33]). It has also been demonstrated that weights can vary over 46 orders of magnitude; most (including the majority of longrange connections) are actually quite weak [34].
The relation of RR effects to localization in the steady state has been studied recently [17]. Here we discuss localization results obtained via a quenched meanfield (QMF) approximation.
In addition to serving as models of brain connectivity, hierarchical modular networks are abundant in nature. They occur in social relations [35], in the WWW [36], metabolic [37] and information systems [38] for example. Thus the results reported here should find application in these and related areas.
Hierarchical modular networks
Recent studies indicate that activity in brain networks persists between the extremes of rapid extinction and global excitation. This socalled limited sustained activity (LSA) does not occur in spreading models defined on structureless, smallworld networks. On the other hand, brain networks are known to exhibit hierarchical modular structure. This motivated Kaiser and Hilgetag (KH) to perform numerical studies on such networks, to investigate topological effects on LSA [39]. Their hierarchical model reflects general features of brain connectivity on large and mesoscopic scales. Nodes in the model were intended to represent cortical columns rather than individual neurons. The connections between them were taken as excitatory, since there appear to be no longdistance inhibitory connections within the cerebral cortex [40].
Kaiser and Hilgetag generated networks beginning at the highest level, and adding modules to the next lower level, with random connectivity within modules. They explored hierarchical networks with different numbers of levels, and numbers of submodules at each level. The average degree (over all nodes in the network) was set to , motivated by experimental studies. They investigated different topologies by varying the edge density across the levels. All the networks studied by KH are of smallworld type, i.e., they have an infinite topological dimension.
The spreading model investigated by KH is a twostate threshold model, in which, at each time step, inactive nodes are activated with probability if at least of their neighbors are currently active, and active nodes deactivate spontaneously with probability . This model is very similar to reactiondiffusion models known in statistical physics [41, 42, 43], with a synchronous cellular automaton (SCA) updates. Starting from an initial state in which a localized set of nodes is active, these authors followed the density of active sites and their localization up to time steps on networks of sizes .
Kaiser and Hilgetag found that LSA can be found in a larger parameter range in HMNs, as compared with random and nonhierarchical smallworld networks. The optimal range of LSA was found in networks in which the edge density increased from the top level () to the bottom (). Such topologies foster activity localization and rareregion effects.
In this paper we investigate HMNs which possess increasing edge density from top to bottom levels, as did KH, but with finite topological dimension. We will show that although localization is seen in small world networks, to observe GPs, with powerlaw dynamics, we need networks of finite topological dimension.
To generate a hierarchical modular network, we define levels on the same set of nodes; on the level we define modules. We achieve this by splitting each module into four equal sized modules on the next level, as if they were embedded in a regular, twodimensional (2d) lattice (see Figure 1). The probability that two nodes are connected by an edge follows as in [39], where is the greatest integer such that the two nodes are in the same module on level . After selecting the number of levels and nodes we fix the average node degree, , and generate the adjacency matrix by proceeding from the highest to the lowest level. We fill the submatrices with zeros and ones and copy them to appropriate diagonal locations of . We allow unidirectional connections, which is more realistic. In fact, disorder associated with the link orientations turns out to be necessary to observe GPs. Finally, we connect the lowestlevel modules with an edge, chosen randomly from nodes of . This provides a shortlinked base lattice that guarantees connectivity.
In a preliminary study, the extra edges, corresponding to the base lattice, were not added. The resulting networks typically consist of a large number of isolated connected components; the GP effects observed in these structures are a consequence of fragmented domains of limited size. Measurements of axonlength distributions in several real neural networks show a large peak at short distances followed by a long flat tail. Thus there is a dense set of local edges in addition to a sparse network of longrange connections, best fit by an exponential function [44]. Typical estimates indicate that of all cortical connections are formed at the local level (i.e., within a radius of mm); only the remainder leave the local gray matter and project to remote targets. Therefore, we also investigate a variant of HMN2d networks in which the lowestlevel modules are fully connected, and there is a single link among the nodes of modules on level . To broaden further the range of structures, we study hierarchical modular trees (HMTs), which possess the minimum number of edges required for connectivity, and so have no loops. Construction of HMTs is described in the Supplementary material.
Relation to BenjaminiBerger networks
Due to the embedding, there is a correspondence with BenjaminiBerger (BB) networks [45]. BB networks are generated on Euclidean lattices, with short links connecting nearest neighbors. In addition, the network contains long links, whose probability decays algebraically with Euclidean distance :
(2) 
Here we consider modified BB networks, in which the long links are added level by level, from top to bottom, as in [39]. The levels : are numbered from the bottom to top. The size of domains, i.e., the number of nodes in a level, grows as in the 4module construction, related to a tiling of the two dimensional base lattice. Due to the approximate distancelevel relation, , the longlink connection probability on level is:
(3) 
Here is related to the average degree of a node .
It is known that in a onedimensional base lattice the BB construction results in a marginal case, , in which the topological dimension is finite and changes continuously with . For we have small world networks, while for the the topological dimension is the same as the base lattice () in BB networks. The HMN1 networks studied by Moretti and Muñoz [24] are similar to the BB model, without the 1d base lattice, but with the inclusion of a HMN topology. These authors simulated spreading models on HMN1 with finite topological dimension at and found GPs. Given the above mapping, this result is not very surprising, because this HMN corresponds to the case, staying close to the percolation threshold of the network. To ensure connectivity, Moretti and Muñoz added extra links to the large connected components [24]; they assert that the topological dimension remains finite, despite the additional links. In networks with finite sizes this is certainly true, however in the infinite size limit such a system is also infinite dimensional, so that we don’t expect true GPs. For the HMN2, the number of connections between blocks at every level is a priori set to a constant value. In this case numerics of [24, 25] suggest GP effects.
Here we embed the HMNs in 2d base lattices, which is closer to the real topology of cortical networks. In this case the threshold for marginality (i.e., continuously changing dimension and exponents) is expected to be at . We confirm this by explicit measurements of the topological dimension. Furthermore, we study twodimensional HMNs with different values and find GP effects for finite topological dimensions, both at the percolation threshold (for ), and for , where the tail distribution of the link lengths decays very rapidly, in agreement with empirical results for axon lengths.
In addition to the CP (), we also consider threshold models (), which are expected to be more realistic for neuronal systems.
Dimension measurements
To measure the dimension of the network we first compute the level structure by running a breadthfirst search from several, randomly selected seeds. We count the number of nodes with chemical distance or less from the seeds and calculate averages over the trials. We determined the dimension of the network from the scaling relation , by fitting a powerlaw to the data for .
At we observe a percolation threshold near , for which the topological dimension is finite: (see Figure 2). Note that the curves with large exhibit saturation, corresponding to a finitesize effect. The case appears to correspond to a marginal dimension, for which, above the percolation threshold , one observes a continuously varying topological dimension (see Figure 3). A more detailed, local slopes analysis via the effective topological dimension
(4) 
is shown in the inset of Figure 3, where and are neighboring measurement points. Saturation to dependent dimensions can be read off the plot for .
Random hierarchical trees also exhibit a finite topological dimension. Figurere 4 shows for a set of ten independent random trees. The average of over the set of ten trees follows a power law to good approximation, with . (Note that in this case is an average over all nodes.) For regular trees, by contrast, grows exponentially with , as shown in Figure 5. Thus the regular tree construction results in a structure with infinite topological dimension.
Dynamic simulations
Contact process
The contact process (CP) [46, 47], a Markov process defined on a network, can be interpreted as a simple model of information spreading in social [48], epidemic spreading [49, 50, 51], or in brain networks [24]. In the CP sites can be active or inactive. Active sites propagate the activity to their neighbors at rate , or spontaneously deactivate with rate (for simplicity).
We perform extensive activitydecay simulations for HMN2d models with and maximum levels: . In these networks we use directed links between nodes, similar to real nervous systems. We follow in runs for each independent random network sample, starting from fully active initial states, and averaging over thousands of independent random network samples for each . In the marginal longlink decay case at , a clear GP behavior is found (see Figure 6).
We show that a GP can occur for more general parameters than those studied in [24], by following the density decay in networks with and (i.e., at the percolation threshold). Sizeindependent, nonuniversal power laws can be seen in Figure 7.
When we increase the average degree , the GP shrinks to a smaller range of values, as in [15], tending toward a simple critical phase transition. However, it is hard to determine at precisely which value of this happens. We note that the connectivity of the network is not assured, without the addition of extra links. Strictly speaking, however, with this extension networks become infinitedimensional in the limit.
More importantly, GPs are also found in networks connected on the base level via short edges. We study the activity decay for , which corresponds to fast decaying tail distribution for the long links, preserving finite topological dimension . In this construction the average degree at the bottom level is , decreasing systematically with , so that for . Note that that the ratio of short to long links is , in agreement with results for real neural networks [44]. Simulations again yield sizeindependent powerlaw decay curves, confirming GP behavior, as shown in Figure 8.
We also determined the effective decay exponent in the usual manner (see [43]), via
(5) 
where and are neighboring data points. The critical point can be located at , showing meanfield scaling: . Above this threshold powerlaws can still be seen for smaller sizes (), up to , but corrections to scaling become stronger; the curves for exhibit saturation. This is surprising, because in a disordered system with short ranged interactions one would expect a critical phase transition with an ultraslow, logarithmic scaling. However, recent studies of the CP in higher dimensions find meanfield criticality and GP [52]. Our result suggests that the topological heterogeneity generated by the long edges is not strong enough to induce an infiniterandomness fixed point. Otherwise we must assume very strong corrections to scaling in this case.
Disorder due to the randomly chosen orientations of the links turns out to be relevant. For symmetric links our simulations yield nonuniversal power laws, which appear to be sensitive to the system size, suggesting the lack of a true GP in the infinitesize limit.
Burstyness in the CP
We study the distribution of interevent times in the CP on asymmetric, HMN2d networks with , in a manner similar to that described in [14]. Starting from fully active states, on networks of size , we measured between successive activation attempts. We followed the evolution up to Monte Carlo steps (MCs), averaging over independent random networks with runs for each. Througout this paper time is measured by MCs. As Figure 9 shows, fattailed distributions, , emerge in the GP for , while decays exponentially outside this region. The slopes of the decay curves are determined via leastsquares fits in the window . For we find , while for the exponent is slightly larger: . These values are close to the autocorrelation exponent of the critical dimensional CP as in [14], but exhibit deviations due to the heterogeneities. This is understandable, since the effective dimension of this HMN2d is close to one. The scaling variable exhibits logperiodic oscillations. This is the consequence of the modular structure of the network. Furthermore, as in other GP models, logarithmic corrections to scaling are expected.
We expect similar nonuniversal, control parameter dependent tails in in the GPs exhibited by the other HMN2d networks. Furthermore, as shown in [14], powerlaw distributions should arise for other initial conditions, such as localized activity. This suggests that bursty intercommunication events in brain dynamics arise spontaneously near the critical point, in the GP.
CP on random hierarchical trees
We simulate the CP with symmetric links on random hierarchical trees (RHTs) of 262144 nodes. We first perform quasistationary (QS) simulations [53], of the CP on a single RHT. For one structure this yields ; for another, independently generated structure of the same kind, we find .
Our principal interest is in the decay of activity starting from all nodes active. The decay of activity on a single RHT appears to follow the scenario familiar from the CP on regular lattices: powerlaw decay is observed at but not at nearby values, as illustrated in Figure 10. The randomness associated with a single RHT appears to be insufficient to generate a Griffiths phase.
By contrast, evidence of a GP is found if we average over many RHTs. The activity averaged over a large set (  ) of independent realizations, each on a different RHT, shown in Figure 11, decays asymptotically as a powerlaw over a fairly broad range of subcritical values. The decay exponent extrapolates to zero at , as shown in the inset. We note that the average is over all realizations, including those that become inactive before the maximum time of MCs. If we instead restrict the average of to trials that survive to time (or greater), the result is a stretched exponential, , where is a constant and the exponent varies with . For , for example, .
Threshold model simulations
Threshold models represent an attempt to capture the activation threshold of real neurons, by requiring at least active neighbors associated with incoming links to activate a node with probability . In case of active nodes spontanous deactivation occurs with probability written in the reactiondiffusion model notation as:
 (a)

with probability ,
 (b)

with probability .
We use SCA updating, as in [39]. Since there is no spatial anisotropy (which might generate activity currents), we can assume that the SCA follow the same asymptotic dynamics as the corresponding model with random sequential updating. Thanks to the synchronous updates, we can speed up the simulations by times by distributing the nodes among multiprocessor cores; by swapping the random number generation on parallel running graphics cards we obtain a further reduction of decrease in the simulation times. The spatiotemporal evolution of the HMN2d network with , in a single network realization, starting from a fully active state, is shown in Figure 12. After a sharp initial drop in activity, due to spontaneous deactivation of nodes with few neighbors, one observes domains (modules) on which activity survives for a long time, suggesting rare region effects.
Results for threshold models
We begin by discussing the case, for which extensive simulations are performed. We generate networks with average degree and . Note that in this case values of higher than the percolation threshold are needed to avoid modules having separated activities. As before, we averaged over hundreds of independent random networks and thousands of independent realizations. We followed the density up to MCs, starting from a fully active initial condition.
In this case the mean activity density decays more rapidly than in meanfield theory, and is sizeindependent (see Figure 13). For large branching probabilities, seems to take a constant value, suggesting an active steady state, but at late times some deviation is observed, possibly due to finitesize effects.
Homogeneous triplet creation models (i.e., ) are expected to exhibit a meanfieldlike discontinuous phase transition in two or more dimensions(see for example [54]). Disorder induces rounding effects, producing continuous phase transitions or GPs [55, 25]. We study a threshold model with and , using . Our results are similar to those for : nonuniversal powerlaws, again suggesting a GP.
Quenched Meanfield approximation for SIS
Heterogeneous meanfield theory provides a good description of network models when fluctuations are irrelevant. This approximation is attractive because is can be solved analytically in many cases; the results agree with simulation [56, 57]. This analysis treats nodes with different degrees as distinct, but finally averages over all degree values, providing a homogeneous solution for the order parameter. To describe the effects of quasistatic heterogeneities in a more precise way, the socalled Quenched MeanField (QMF) approximation was introduced [58, 59]. For SIS models this leads to a spectral analysis of the adjacency matrix of the network. The susceptibleinfectedsusceptible (SIS) [60] model is similar to the CP, except that branching rates are not normalized by , leading to symmetric governing equations. A rate equation for the SIS model can be set up for the vector of activity probabilities of node at time :
(6) 
Here the are weights attributed to the edges. For large times the SIS model evolves to a steady state, with an order parameter . Since this equation is symmetric under the exchange , a spectral decomposition can be performed on a basis of orthonormal eigenvectors .
(7) 
Nonnegativity of the matrix assures a unique, real, nonnegative largest eigenvalue .
In the longtime limit the system evolves into a steady state and we can express the solution via as
(8) 
Stability analysis shows that above a threshold , with an order parameter . In the eigenvector basis equation (8) can be expanded by the coefficients as
(9) 
and near the threshold we can express with the principal eigenvector. In the QMF approximation and the order parameter can be approximated by the eigenvectors of the largest eigenvalues
(10) 
where with the coefficients
(11) 
As proposed in [59], and tested on several SIS network models [17], localization in the active steady state can be quantified by calculating the inverse participation ratio (IPR) of the principal eigenvector, related to the eigenvector of the largest eigenvalue of the adjacency matrix as
(12) 
This quantity vanishes in the case of homogeneous eigenvector components, but remains finite as if activity is concentrated on a finite number of nodes. Numerical evidence has also been provided for the relation of localization to RR effects, slowing down the dynamics to follow powerlaws [17].
We study localization of the SIS in the HMN2d models introduced in the previous section*s. We analyze the eigenvectors of for (), varying , using system sizes ranging from to . In particular, we performed a FSS scaling study of the IPR in these systems. First we determined the behavior on the small world network of [39] corresponding to . As one can see in Figure 14, the localization at small sizes disappears as . By contrast, for and , the graphs are finitedimensional, and the IPR increases with , tending to a finite limiting value, suggesting localization.
One might question the relevance of network models with small connectivity to mammalian brains, in which is on the order of . To answer this we study models with higher connectivity, i.e., . As Figure 15 shows, the sign of localization, which is weak but nonzero for , now disappears. Next, we add random weights , distributed uniformly over the interval , to the links of the networks. A consequence of this explicit disorder is localization even in highly connected networks, as shown in Figure 15. This result is in agreement with the expectation of limited sustained activity in brain networks [39], meaning that the link disorder prevents overexcitation of a network of high connectivity.
Localization suggests rareregion effects, thus a dynamic GP. Nevertheless, simulations of the CP on such networks do not show extended regions of powerlaws for , but rather a nontrivial (nonmeanfield) continuous phase transition (Figure 16). Decay simulations for yield at , albeit with a cutoff due to the finite system size. This agrees with our result for 1D BB networks with [15]. Spreading simulations starting from single, randomly placed seed result in the survival probability scaling at this critical point: (that is, the symmetry holds to within uncertainty). Here the density increases initially as . One can understand these results, noting that the localization values are rather small, , so that RR effects are too weak to generate observable GP effects in the dynamics.
Naturally, for networks, where already the topological disorder is strong enough to create a GP, inhomogeneities in the interaction strengths amplify the RR effects. This induces GPs even for larger values of , as well as in the threshold models. Further studies of disorder effects are under way.
Conclusions
We investigate the dynamical behavior of spreading models on hierarchical modular networks embedded in twodimensional space, with long links whose probability decays as a powerlaw with distance. This corresponds to an exponentially decreasing connection probability, as a function of the level in the hierarchical construction. The aim of this study is to understand the effects of intrinsic and topological disorder, in particular, extended critical regions in the parameter space, without any (self) tuning mechanism.
If we eliminate the underlying lattice structure we observe powerlaw dynamics for networks near the percolation threshold, when the effective dimension is finite. However, sizeindependent powerlaws, corresponding to GPs, are only seen if we have directed links. Since connectivity of these structures is not guaranteed, we also study random hierarchical trees, with full connectivity. GPs are observed upon averaging over many independent network realizations. The relation to brain networks can be envisaged in the largesize limit, if we regard independent realizations as (almost decoupled) submodules of the entire brain.
When we ensure connectivity via short edges on the lowest level, we find GPs for rapidly decaying longlink probabilities. Both of these network assumptions are in accordance with empirical brain networks. Above the GP, at the critical point, we find meanfieldlike decay of activity, in agreement with results on the CP with quenched disorder on higherdimensional, regular lattices. We have also shown that bursty behavior arises naturally in the GPs, due to autocorrelations which decay via a powerlaw.
We perform a quenched meanfield study of the SIS model on these networks; in the SIS, nodes are connected symmetrically. Finite size scaling of the inverse participation ratio suggests localization on largeworld networks and delocalization on small world structures. However, when we add intrinsic weight disorder, localization can be seen even on smallworld networks. Weight disorder is again to be expected in real brain networks, since the strength of couplings varies over many orders of magnitude. Despite this, we saw no GP effects in the dynamics of weighted CPs with . Instead, we find a nontrivial critical scaling, as has already been observed in other networks. The possibility of a narrow GP in this case is an open issue.
In conclusion, we believe our synthetic HMN2ds are closer to experimental brain networks than previously proposed models, and find numerical evidence for GPs in extended phases in simple models with spreading dynamics. Although we eliminate any selftuning mechanisms, we still find nontrivial slow dynamics as well as localization of activity, which is crucial for understanding real brain network data.
Acknowledgments
We thank R. Juhász and C. C. Hilgetag for useful discussions. Support from the Hungarian research fund OTKA (Grant No. K109577) and the European Social Fund through project FuturICT.hu (grant no.: TAMOP4.2.2.C11/1/KONV20120013) is acknowledged. R.D. is supported in part by CNPq, Brazil. The SCA simulations are accelerated by grapics cards donated by NVIDIA within the professor partnership programme. Géza Ódor thanks the Physics Department at UFMG, where part of this work was done, for its hospitality.
Author contributions statement
Géza Ó. wrote, ran and analysed the CP and SIS simulations on HMN2ds and performed the QMF analysis. R. D. wrote, ran and analyzed the HMT simulations. Gergely Ó. wrote and ran the HMN2 threshold model simulations and dimension measurements on them. Géza Ó. and R. D. wrote the main manuscript text. Géza Ó. prepared figures 13, 69, 1316, R. D. prepared figures 4,5,10,11,1720. Gergely Ó. prepared figure 12. All authors reviewed the manuscript.
Additional information
The authors declare no competing financial interests.
References
 1. Legenstein R. and Maass W., New Directions in Statistical Signal Processing: From Systems to Brain (eds S. Haykin, J. C. Principe, T. Sejnowski, J. McWhirter) 127 V154 (MIT Press, 2008).
 2. Beggs J. and Plenz D., Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 1116711177 (2003).
 3. Tetzlaff C, Okujeni S, Egert U, Wörgötter F, Butz M, SelfOrganized Criticality in Developing Neuronal Networks. PLoS Comput. Biol. 6, e1001013 (2010).
 4. Hahn G et al., Neuronal avalanches in spontaneous activity in vivo. J. Neurophysiol. 104, 33123322 (2010).
 5. Ribieiro T. L. et al., Spike Avalanches Exhibit Universal Dynamics across the SleepWake Cycle. PLoS ONE 5, e14129 (2010).
 6. Bédard C, Kröger H., and Destexhe A, Does the 1/f Frequency Scaling of Brain Signals Reflect SelfOrganized Critical States? Phys. Rev. Lett. 97, 118102 (2006).
 7. Dehghani N et al., Avalanche Analysis from Multielectrode Ensemble Recordings in Cat, Monkey, and Human Cerebral Cortex during Wakefulness and Sleep. Front. Physiol. 3, 302 (2012).
 8. Priesemann V, Valderrama M, Wibral M and Le Van Quyen M., Neuronal Avalanches Differ from Wakefulness to Deep Sleep VEvidence from Intracranial Depth Recordings in Humans. PLoS Comput. Biol. 9, e100298 (2014).
 9. Beggs J. M. and Timme N., Being Critical of Criticality in the Brain. Front. Physiol. 3, 163 (2012).
 10. Bak P., Tang C. and Wiesenfeld K., Selforganized criticality. Phys. Rev. A 38, 364 V374 (1988).
 11. Pruessner G., Self Organized Criticality, (Cambridge University Press, 2012).
 12. Vojta T., Rare region effects at classical, quantum and nonequilibrium phase transitions J. Physics A: Math. and Gen. 39, R143R205 (2006).
 13. Griffiths R. B., Nonanalytic Behavior Above the Critical Point in a Random Ising Ferromagnet. Phys. Rev. Lett. 23, 1719 (1969).
 14. Ódor, G. Slow, bursty dynamics as a consequence of quenched network topologies. Phys. Rev. E 89, 042102 (2014).
 15. Muñoz M. A., Juhász R., Castellano C., and Ódor G., Griffiths Phases on Complex Networks. Phys. Rev. Lett. 105, 128701 (2010).
 16. Johnson S., Torres J. J., and Marro J., Robust ShortTerm Memory without Synaptic Learning. PLoS ONE 8, e50276 (2013).
 17. Ódor G., Spectral analysis and slow spreading dynamics on complex networks. Phys. Rev. E 88, 032109 (2013).
 18. Newman M. E. J., Networks: An Introduction, (Oxford Univ. Press, 2010).
 19. Sporns O, Chialvo D. R., Kaiser M. and Hilgetag C. C, Organization, development and function of complex brain networks. Trends Cogn. Sci. 8, 418 V425 (2004).
 20. Sporns O., Networks Of The Brain, (MIT Press, 2010)
 21. Kaiser M, A tutorial in connectome analysis: topological and spatial features of brain networks. NeuroImage 57, 892 V907 (2011).
 22. Meunier D., Lambiotte R. and Bullmore E., Modular and hierarchically modular organization of brain networks. Front. Neurosci. 4, 200 (2010).
 23. Hilgetag C.C., Burns G.A., O’Neill M.A., Scannell J.W., Young M.P. Anatomical connectivity defines the organization of clusters of cortical areas in the macaque monkey and the cat. Phil. Trans. R. Soc. Lond. B 355, 91100 (2000).
 24. Moretti P., Muñoz M. A., Griffiths phases and the stretching of criticality in brain networks Nat. Commun. 4, 2521 (2013).
 25. Villa Martin P. M., Moretti P. and and Muñoz M. A., Rounding of abrupt phase transitions in brain networks. J. Stat. Mech. P01003 (2015).
 26. Humpries M. D., Gurney K. and Prescott T.J., The brainstem reticular formation is a smallworld, not scalefree, network. Proc. Biol. Sci. 273, 503 V511 (2006).
 27. Yu S., Huang D., Singer W., Nikolic D., A Small World of Neuronal Synchrony. Cereb. Cortex 18, 2891 V2901 (2008).
 28. Rosenblum M. G., Pikovsky A. and Kurths J., Synchronization V A Universal Concept In Nonlinear Sciences (Cambridge University Press, 2001).
 29. Kuramoto Y., Selfentrainment of a population of coupled nonlinear oscillators. Lect. Notes Phys. 39, 420422 (1975).
 30. Villegas P., Moretti P., Muñoz M. A., Frustrated hierarchical synchronization and emergent complexity in the human connectome network. Sci. Rep. 4, 5990 (2014).
 31. Ravasz E. and Barabási A., Hierarchical organization in complex networks., Phys. Rev. E 67, 026112 (2003).
 32. Ódor G. and PastorSatorras R., Slow dynamics and rareregion effects in the contact process on weighted tree networks. Phys. Rev. E 86, 026117 (2012).
 33. Felleman D. J. and Van Essen D. C., Distributed hierarchical processing in the primate cerebral cortex. Cereb. Cortex 1, 147 (1991).
 34. Markov N. T. et al., Weight consistency specifies regularities of macaque cortical networks. Cereb. Cortex 21, 125472 (2010).
 35. Granovetter M. S., The strength of weak ties. Am. J. Sociol. 78, 13601380 (1973).
 36. Flake G. W., Lawrence S., and Giles C.L., Efficient identification of web commnunities Proceedings of the Sixth International Conference on Knowledge Discovery and Data Mining [150160] (ACM, Boston, MA 2000).
 37. Ravasz E., Somera A. L., Mongru D. A., Oltvai Z. N. and Barabási A.L., Hierarchical organization of modularity in metabolic networks. Science 297, 15511555 (2002).
 38. Vazquez A., PastorSatorras R. and Vespignani A., Largescale topological and dynamical properties of the Internet. Phys. Rev. E 65, 066130 (2002).
 39. Kaiser M. and Hilgetag C. C., Optimal hierarchical modular topologies for producing limited sustained activation of neural networks. Front. in Neuroinform. 4, 8 (2010).
 40. Latham P. E. and Nirenberg S., Synergy, Redundancy, and Independence in Population Codes, Revisited. Neural Comp. 16,, 13851412 (2004).
 41. Marro J. and Dickman R., Nonequilibrium Phase Transitions in Lattice Models, (Cambridge University Press, Cambridge, 1999).
 42. Ódor, G., Universality classes in nonequilibrium lattice systems. Rev. Mod. Phys. 76, 663724 (2004).
 43. Ódor G., Universality in Nonequilibrium Lattice Systems, (World Scientific, Singapore, 2008).
 44. Kaiser M., Hilgetag C. C. and van Ooyen A., A simple rule for axon outgrowth and synaptic competition generates realistic connection lengths and filling fractions. Cereb. Cortex 19, 30013010 (2009).
 45. Benjamini I. and Berger N.,The diameter of longrange percolation clusters on finite cycles. Rand. Struct. Alg. 19, 102111 (2001).
 46. Harris T. E., Contact Interactions on a Lattice, Ann. Prob. 2, 969988 (1974).
 47. Liggett T. M., Interacting Particle Systems, (SpringerVerlag, 1985, Berlin).
 48. PastorSatorras R. and Vespignani A., Evolution and Structure of the Internet: A Statistical Physics Approach (Cambridge University Press, Cambridge, 2004).
 49. van Ballegooijen W. M. and Boerlijst M. C., Emergent tradeoffs and selection for outbreak frequency in spatial epidemics. PNAS. USA 101, 1824618250 (2004).
 50. Sun G.Q. et al., Influence of infection rate and migration on extinction of disease in spatial epidemics. J. of Theor. Bio. 264, 95103 (2010).
 51. Sun G.Q. et al., Phase transition in spatial epidemics using cellular automata with noise. Ecol. Res. 26, 333340 (2011).
 52. Vojta T., Igo J., Hoyos J. A., Rare regions and Griffiths singularities at a clean critical point: The fivedimensional disordered contact process. Phys. Rev. E 90, 012139 (2014).
 53. Oliveira M. M. and Dickman R., How to simulate the quasistationary state Phys. Rev. E 71, 016129 (2005).
 54. Ódor, G. Phase transition classes in triplet and quadruplet reactiondiffusion models. Phys. Rev. E 67, 056114 (2003).
 55. V. Martin P. M., Bonachela J. A. and Muñoz M.A., Quenched disorder forbids discontinuous transitions in nonequilibrium lowdimensional systems. Phys. Rev. E 89, 012145 (2014).
 56. Boguñá M., Castellano C., and PastorSatorras R., Langevin approach for the dynamics of the contact process on annealed scalefree networks. Phys. Rev. E 79, 036110 (2009).
 57. Mata A. S., Ferreira R. S., Ferreira S. C., Heterogeneous pairapproximation for the contact process on complex networks. New J. Phys. 16, 053006 (2014).
 58. Chakrabarti D., Wang Y., Wang C., Leskovec J., and Faloutsos C., Epidemic thresholds in real networks. ACM Trans. Inf. Syst. Secur. 10, 126 (2008).
 59. Goltsev A. V., Dorogovtsev S. N., Oliveira J. G., and Mendes J. F. F., Localization and Spreading of Diseases in Complex Networks. Phys. Rev. Lett. 109, 128702 (2012).
 60. Bailey N. T. J., The Mathematical Theory of Infectious Diseases 2nd ed. (Griffin, London, 1975).
Supplementary information
Although the set of nodes is isomorphic to a square lattice of sites, we shall label the nodes not by their Cartesian coordinates but rather using a base4 notation. At level , the full set of nodes is divided into four quadrants labeled 0, 1, 2, and 3, proceeding counterclockwise from the lower left. At level , each of the quadrants is similarly divided into four subquadrants, labeled in the same manner. Division into subquadrants continues all the way down to level 1, in which the elements are individual nodes. A given site can be specified by the labels of the quadrant, subquadrant, etc. to which it belongs. Denoting the labels at levels 1, 2,…, L, by , the address of a site is (see Fig. S17).
Consider a set of four nodes, 0, 1, 2, and 3. There are six possible links between them: (0,1), (0,2), (0,3), (1,2), (1,3), and (2,3). A tree with four nodes has exactly three links. Of the subsets of three links, sixteen yield a connected tree [for example: {(0,1), (1,2), (2,3)}], while four correspond to a cycle of three nodes, leaving the fourth node isolated [example: {(1,2), (1,3), (2,3)}]. To construct a tree of four nodes, we choose one of the sixteen link sets at random.
A hierarchical tree is constructed by first linking the four quadrants via three edges. The link set is chosen at random. Then, for each link, we choose sites at random within each of the two quadrants connected by the link, to serve as the connected nodes. Now we repeat this process within each of the subquadrants, and so on, until we reach the basic modules of four sites. The latter are again connected by sets of three links, chosen independently from the basic collection of sixteen link sets. At the end, each basic module is connected internally, and to a module at the next level, etc., so that we have a connected graph of nodes and edges. Although we have chosen to begin the construction at level , we could equally have begun at level 1, or some intermediate level: the choices of link sets and nodes at the various levels are mutually independent. These steps are illustrated in the Figs. S18, S19, S20.
At level 1 the number of links per node is 3/4; at level , there are blocks connected via links. The number of links per node at this level is therefore . Let denote the probability that a randomly chosen edge linking blocks at level be present. At level 1, links connect nodes within foursite modules. Since there are six possible links, and since just three are present within each module, we have . At level 2, a link connects nodes within a 16node module; at this level a node is linked to one of the 12 nodes outside its basic 4site module. There are possible links, and since only three are present, . At level , an edge links nodes within a module containing nodes. The number of possible links at this level is , so that , and, in general, . Since a tree represents a minimally connected structure, cannot decay faster than this rate, in any connected hierarchical network based on fournode modules.
Although our principal interest is in random trees, as described above, we may also consider regular trees; these are constructed by using the same link set at each level. The resulting structure turns out to have infinite topological dimension, as discussed in Section ”Dimension measurements”.