Griffiths phases and localization in hierarchical modular networks

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 mean-field theory (QMF), focusing on structures with a connection probability that decays exponentially with the level index. Such networks can be embedded in two-dimensional 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 small-world networks originally proposed to model brain connectivity, the topological heterogeneities are not strong enough to induce deviations from mean-field 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, power-law 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 small-world networks with high connectivity in case of link disorder.


I. 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 self-tuning is well known in statistical physics; the paradigm of self-organized 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 absorbing-state 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 rare-region (RR) effects [12] that smear the phase transitions.These effects can make a discontinuous transition continuous, or generate so-called Griffiths phases (GP) [13], in which critical-like power-law dynamics appears over an extended region around the critical point.
In these regions, moreover, non-Markovian, bursty behavior can emerge as a consequence of a diverging correlation time.The inter-communication times of the nodes (which possess no internal memory) follow a fat-tailed distribution [14].Thus, heterogeneities widen the critical region and can contribute to the occurrence of power laws.This provides an alternative mechanism for critical-like 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 finite-dimensional systems [15,16] and can be a reason for the working memory in the brain [17].Although many networks exhibit the small-world 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 [18].
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 [19].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 large-scale units such as visual or sensory-motor cortex.At each level, connections within modules are denser than between different modules [20][21][22][23].Although empirical data confirm this modular organization on some scales [24], the detailed organization of brain networks is not yet experimentally accessible.
A particular kind of hierarchical modular network (HMN) model was proposed and investigated numerically and analytically in [25].On a large-world HMN, which implies a finite topological dimension D, models of the spread of activity exhibit power-law dynamics and rare region effects.However, these power laws are system-size 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 [20] and the synchronization networks of cortical neurons [27] suggest small-world topology [26].The brain network modules of [25] are weakly coupled in such a way that these HMNs are near the percolation threshold, as in the case of the models introduced in [16].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 D 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 tree-like 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 [25].In this case even weaker rare-region 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 sub-modules 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 where n i denotes the number of direct links interconnecting the k i nearest neighbors of node i, exhibits different scaling behavior in modular networks than in unstructured networks.
While in SF networks it decays as C(N) ∼ N −3/4 and in random networks as C(N) ∼ N −1 , in hierarchical modular networks (HMNs) it is constant.Furthermore, while in SF networks the clustering is degree-number independent, in HMNs it decays as C(k) ∼ k −β , 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 4-6 orders of magnitude; most (including the majority of long-range connections) are actually quite weak [34].
The relation of RR effects to localization in the steady state has been studied recently [18,35,36].Here we discuss localization results obtained via a quenched mean-field (QMF) approximation.
In addition to serving as models of brain connectivity, hierarchical modular networks are abundant in nature.They occur in social relations [37], in the WWW [38], metabolic [39] and information systems [40] for example.Thus the results reported here should find application in these and related areas.

II. HIERARCHICAL MODULAR NETWORKS
Recent studies indicate that activity in brain networks persists between the extremes of rapid extinction and global excitation.This so-called limited sustained activity (LSA) does not occur in spreading models defined on structureless, small-world 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 [41].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 long-distance inhibitory connections within the cerebral cortex [42].Kaiser and Hilgetag found that LSA can be found in a larger parameter range in HMNs, as compared with random and non-hierarchical small-world networks.The optimal range of LSA was found in networks in which the edge density increased from the top level (l = l max ) to the bottom (l = l 1 ).Such topologies foster activity localization and rare-region 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 power-law dynamics, we need networks of finite topological dimension.
To generate a hierarchical modular network, we define l max levels on the same set of N = 4 lmax nodes; on the l th level we define 4 l 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, two-dimensional (2d) lattice (see Fig. 1).The probability p l that two nodes are connected by an edge follows p l ≃ b2 −sl as in [41], where l is the greatest integer such that the two nodes are in the same module on level l.After selecting the number of levels and nodes we fix the average node degree, b = k /2, and generate the adjacency matrix A 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 A. 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 lowest-level modules with an edge, chosen randomly from nodes of l 1 .This provides a short-linked 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 axon-length distributions in several real neural networks show a large 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 long-range connections [46].
Typical estimates indicate that ≃ 90% of all cortical connections are formed at the local level (i.e., within a radius of 0.5mm); 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 lowest-level modules are fully connected, and there is a single link among the nodes of modules on level l 2 .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 Appendix.

A. Relation to Benjamini-Berger networks
Due to the embedding, there is a correspondence with Benjamini-Berger (BB) networks [47].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 R: Here we consider modified BB networks, in which the long links are added level by level, from top to bottom, as in [41].The levels : l = 1, ..., l max are numbered from the bottom to top.The size of domains, i.e., the number of nodes in a level, grows as N l = 4 l in the 4-module construction, related to a tiling of the two dimensional base lattice.Due to the approximate distance-level relation, R ≃ 2 l , the long-link connection probability on level l is: Here b is related to the average degree of a node k .
It is known that in a one-dimensional base lattice the BB construction results in a marginal case, s = 2, in which the topological dimension is finite and changes continuously with b.For s < 2 we have small world networks, while for s > 2 the the topological dimension is the same as the base lattice (d = 1) in BB networks.The networks studied by Moretti and Muñoz [25] 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 HMNs with finite topological dimension at p l ∼ ( 1 4 ) l and found GPs.Given the above mapping, this result is not very surprising, because this HMN corresponds to the s = 2 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 [25]; 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.
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 s = 4.We confirm this by explicit measurements of the topological dimension.Furthermore, we study two-dimensional HMNs with different s values and find GP effects for finite topological dimensions, both at the percolation threshold (for s = 3), and for s = 6, where the tail distribution of the link lengths is exponential.In addition to the CP (m = 1), we also consider threshold models (m ≥ 2), which are expected to be more realistic for neuronal systems.

B. Dimension measurements
To measure the dimension of the network we first compute the level structure by running a breadth-first search from several, randomly selected seeds.We count the number of nodes N r with chemical distance r or less from the seeds and calculate averages over the trials.
We determined the dimension of the network from the scaling relation N ∼ r D , by fitting a power-law to the data for N(r).At s = 3 we observe a percolation threshold near k = 4, for which the topological dimension is finite: D ≃ 1.4 (see Fig. 2).Note that the curves with large k exhibit saturation, corresponding to a finite-size effect.The case s = 4 appears to correspond to a marginal dimension, for which, above the percolation threshold k > 6, one observes a continuously varying topological dimension (see Fig. 3).A more detailed, local slopes analysis via the effective topological dimension is shown in the inset of Fig. 3, where r and r ′ are neighboring measurement points.Saturation to s-dependent dimensions can be read off the plot for r > 100.Random hierarchical trees also exhibit a finite topological dimension.Figure 4 shows N r for a set of ten independent random trees.The average of N r over the set of ten trees follows a power law to good approximation, with D ≃ 0.72.(Note that in this case N r is an average over all nodes.)For regular trees, by contrast, N r grows exponentially with r, as shown in Fig. 5. Thus the regular tree construction results in a structure with infinite topological dimension.

III. DYNAMIC SIMULATIONS A. Contact process
The contact process (CP) [48,49], a Markov process defined on a network, can be interpreted as a simple model of information spreading in social [50] or in brain networks [25].In the CP sites can be active or inactive.Active sites propagate the activity to their neighbors at rate λ/k, or spontaneously deactivate with rate ν = 1 (for simplicity).
We perform extensive activity-decay simulations for HMN2d models with 2b = k = 4 and maximum levels: l max = 8, 9, 10.In these networks we use directed links between nodes, similar to real nervous systems.We follow ρ(t) in 10−100 runs for each independent random network sample, starting from fully active initial states, and averaging ρ(t) over thousands of independent random network samples for each λ.In the marginal long-link decay case at s = 4, a clear GP behavior is found (see Fig. 6).System sizes l max = 8, 9, 10 (thin, medium, and thick lines, respectively).Size-independent power laws are evidence of a Griffiths phase.
We show that a GP can occur for more general parameters than those studied in [25], by following the density decay in networks with s = 3 and k = 4 (i.e., at the percolation threshold).Size-independent, nonuniversal power laws can be seen in Fig. 7.
When we increase the average degree k , the GP shrinks to a smaller range of λ values, as in [16], tending toward a simple critical phase transition.However, it is hard to determine at precisely which value of k 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 s ≤ 4 networks become infinite-dimensional in the N → ∞ limit.
More importantly, GPs are also found in networks connected on the base level via short edges.We study the activity decay for s = 6, which corresponds to an exponential tail distribution for the long links, preserving finite topological dimension D. In this construction the average degree at the bottom level is k ≃ 6, decreasing systematically with l, so that k ≃ 1 for l max .Note that that the ratio of short to long links is ∼ 0.11, in agreement with results for real neural networks [46].Simulations again yield size-independent power-law decay curves, confirming GP behavior, as shown in Fig. 8. FIG.7: CP on asymmetric HMN2d networks as in Fig. 6, but for s = 3.
We also determined the effective decay exponent in the usual manner (see [45]), via where ρ(t) and ρ(t ′ ) are neighboring data points.The critical point can be located at λ c = 2.53 (1), showing mean-field scaling: ρ(t) ∼ t −1 .Above this threshold power-laws can still be seen for smaller sizes (l max = 8, 9), up to λ ≃ 2.55, but corrections to scaling become stronger; the curves for l max = 10 exhibit saturation.This is surprising, because in a disordered system with short ranged interactions one would expect a critical phase transition with an ultra-slow, logarithmic scaling.However, recent studies of the CP in higher dimensions find mean-field criticality and GP [51].Our result suggests that the topological heterogeneity generated by the long edges is not strong enough to induce an infinite-randomness 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 infinite-size limit.

B. Burstyness in the CP
We study the distribution of inter-event times in the CP on asymmetric, HMN2d networks with s = 6, in a manner similar to that described in [14].Starting from fully active states, on networks of size l max = 10, we measured ∆t between successive activation attempts.We followed the evolution up to t max = 2 18 MCs, averaging over 1000 − 2000 independent random networks with 10 − 100 runs for each.As Fig. 9 shows, fat-tailed distributions, P (∆t) ∼ (∆t) −x , emerge in the GP for 2.46 < λ < 2.6, while P (∆t) decays exponentially outside this region.The slopes of the decay curves are determined via least-squares fits in the window 20 < ∆t < 7000.For λ = 2.5 we find x = 1.753(4), while for λ = 2.52 the exponent is slightly larger: x = 1.81 (1).These values are close to the auto-correlation exponent of the critical 1 + 1 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 P * (∆t) ≡ (∆t) 1.79 P (∆t) exhibits log-periodic 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 P (∆t) in the GPs exhibited by the other HMN2d networks.Furthermore, as shown in [14], power-law distributions should arise for other initial conditions, such as localized activity.This suggests that bursty inter-communication events in brain dynamics arise spontaneously near the critical point, in the GP.

C. CP on random hierarchical trees
We simulate the CP with symmetric links on random hierarchical trees (RHTs) of 262 144 nodes.We first perform quasistationary (QS) simulations [52], of the CP on a single RHT.
For one structure this yields λ c ≃ 2.72; for another, independently generated structure of the same kind, we find λ c ≃ 2.76.
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: power-law decay is observed at λ c but not at nearby values, as illustrated in Fig. 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 (10 3 -10 4 ) of independent realizations, each on a different RHT, shown in Fig. 11, decays asymptotically as a power-law over a fairly broad range of subcritical λ values.The decay exponent α extrapolates to zero at λ = 2.760(2) ≃ λ c , as shown in the inset.We note that the average is over all realizations, including those that become inactive before the maximum time of 2 × 10 6 MCs.If we instead restrict the average of ρ(t) to trials that survive to time t (or greater), the result is a stretched exponential, where C is a constant and the exponent β varies with λ.For λ = 2.70, for example, β ≃ 0.25.

D. Threshold model simulations
Threshold models represent an attempt to capture the activation threshold of real neurons, by requiring at least m active neighbors associated with incoming links to activate a node.We use stochastic cellular automata (SCA) (with synchronous updates), as in [41].
Since there is 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.The spatio-temporal evolution of the HMN2d network with m = 6, in a single network realization, starting from a fully active state, is shown in Fig. 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.In Fig. 12 blue (darker) colors represent higher, while green and yellow

Results for threshold models
For large m, the probability of finding a Griffiths phase is small; we therefore begin by discussing the m = 2 case, for which extensive simulations are performed.We generate networks with average degree k = 24 and s = 6.Note that in this case values of k 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 t max = 10 5 MCs, starting from a fully active initial condition.
In this case the mean activity density decays more rapidly than in mean-field theory, and is size-independent (see Fig. 13).For large branching probabilities, ρ(t) seems to take a constant value, suggesting an active steady state, but at late times some deviation is observed, possibly due to finite-size effects.
Homogeneous triplet creation models (i.e., m = 3) are expected to exhibit a mean-field- like discontinuous phase transition in two or more dimensions(see for example [53]).Disorder induces rounding effects, producing continuous phase transitions or GPs [54].We study a threshold model with m = 3 and s = 6, using k = 36.Our results are similar to those for m = 2: non-universal power-laws, again suggesting a GP.

IV. QUENCHED MEAN-FIELD APPROXIMATION FOR SIS
Heterogeneous mean-field 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 [55][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 quasi-static heterogeneities in a more precise way, the so-called Quenched Mean-Field (QMF) approximation was introduced [58][59][60][61].For SIS models this leads to a spectral analysis of the adjacency matrix A ij of the network.The susceptible-infected-susceptible (SIS) [62] model is similar to the CP, except that branching rates are not normalized by k, leading to symmetric governing equations.A rate equation for the SIS model can be set up for the vector of activity probabilities ρ i (t) of node i at time t: Here the w ij = w ji are weights attributed to the edges.For large times the SIS model evolves to a steady state, with an order parameter ρ ≡ ρ i .Since this equation is symmetric under the exchange i ↔ j, a spectral decomposition can be performed on a basis of orthonormal eigenvectors e(y).Non-negativity of the matrix B ij ≡ A ij w ij assures a unique, real, nonnegative largest eigenvalue y M .
In the long-time limit the system evolves into a steady state and we can express the solution via B ij as Stability analysis shows that ρ i > 0 above a threshold λ c , with an order parameter ρ ≡ ρ i .In the QMF approximation one can find λ c and ρ(λ) for λ ≃ λ c from the principal eigenvector; in particular, 1/λ c = y M .
As proposed in [61], and tested on several SIS network models [18,35,36], 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 e(y M ) of the adjacency matrix as This quantity vanishes ∼ 1/N in the case of homogeneous eigenvector components, but remains finite as N → ∞ 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 power-laws [18,35].
We study localization of the SIS in the HMN2d models introduced in the previous sections.
We analyze the eigenvectors of B ij for b = 2 ( k ≃ 4), varying s, using system sizes ranging from N = 256 to N = 262144.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 [41] corresponding to s = 3/2.As one can see in Fig. 14, the localization at small sizes disappears as I(N) ∼ 1/N 3 .By contrast, for s = 3 and s = 4, the graphs are finite-dimensional, and the IPR increases with N, tending to a finite limiting value, suggesting localization.
One might question the relevance of network models with small connectivity to mammalian brains, in which k is on the order of 10 3 .To answer this we study s ≤ 4 models with higher connectivity, i.e., k ≃ 50.As Fig. 15 shows, the sign of localization, which is weak but nonzero for k = 4, now disappears.Next, we add random weights w ij , distributed uniformly over the interval (0, 1), to the links of the networks.A consequence of this explicit disorder is localization even in highly connected networks, as shown in Fig. 15.This result is in agreement with the expectation of limited sustained activity in brain networks [41], meaning that the link disorder prevents over-excitation of a network of high connectivity.
Localization suggests rare-region effects, thus a dynamic GP.Nevertheless, simulations of FIG.14: Finite size scaling of the inverse participation ratio in weakly coupled HMN2d models, with maximum levels l max = 4, 5, ..9.The s = 3 (bullets) and s = 4 (boxes) results suggest localization (finite I) in the infinite-size limit.Lines are power-law fits to the data.For s = 1.5, corresponding to the symmetrized, small-world network (model-6 of [41]) no evidence of localization is seen.Lines show power-law fits to the data.In the unweighted case, no localization effect can be seen, and I decays linearly with N .
the CP on such networks do not show extended regions of power-laws for s ≤ 4, but rather a nontrivial (non-mean-field) continuous phase transition (Fig. 16).Decay simulations for l = 9 yield ρ(t) ≃ t −0.50 (1) at λ c ≃ 2.588(1), albeit with a cutoff due to the finite system size.This agrees with our result for 1D BB networks with s = 3 [16].Spreading simulations starting from single, randomly placed seed result in the survival probability scaling at this critical point: P (t) ≃ t −0.50 (1) (that is, the symmetry α = δ holds to within uncertainty).
Here the density increases initially as ρ(t) ≃ t 0.23 (1) .One can understand these results, noting that the localization values are rather small, I ∼ 0.08, so that RR effects are too weak to generate observable GP effects in the dynamics.Naturally, for s > 4 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 k , as well as in the threshold models.Further studies of disorder effects are under way.

V. CONCLUSIONS
We investigate the dynamical behavior of spreading models on hierarchical modular networks embedded in two-dimensional space, with long links whose probability decays as a power-law 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.
When we do not maintain the underlying lattice structure we observe power-law dynamics for networks near the percolation threshold, when the effective dimension is finite.However, size-independent power-laws, 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 large-size limit, if we regard independent realizations as (almost decoupled) sub-modules of the entire brain.
When we ensure connectivity via short edges on the lowest level, we find GPs for rapidly decaying long-link probabilities.Both of these network assumptions are in accordance with empirical brain networks.Above the GP, at the critical point, we find mean-field-like decay of activity, in agreement with results on the CP with quenched disorder on higher-dimensional, regular lattices.We have also shown that bursty behavior arises naturally in the GPs, due to autocorrelations which decay via a power-law.
We perform a quenched mean-field study of the SIS model on these networks; in the SIS, nodes are connected symmetrically.Finite size scaling of the inverse participation ratio suggest localization on large-world networks and de-localization on small world structures.
However, when we add intrinsic weight disorder, localization can be seen even on small-world 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 s=4.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 self-tuning mechanisms, we still find nontrivial slow dynamics as well as localization of activity, which crucial for understanding real brain network data.
FIG. 17: Definition of node labels for a network with L = 2 A hierarchical tree is constructed by first linking the four quadrants via three edges.The link set B i 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 N At level 1 the number of links per node is 3/4; at level j, there are 4 L−j blocks connected via 4 L−j −1 links.The number of links per node at this level is therefore (4 L−j −1)/4 L ≃ 1/4 j .
Let p j denote the probability that a randomly chosen edge linking blocks at level j be present.
At level 1, links connect nodes within four-site modules.Since there are six possible links, and since just three are present within each module, we have p 1 = 1/2.At level 2, a link connects nodes within a 16-node module; at this level a node is linked to one of the 12 nodes outside its basic 4-site module.There are (16 • 12)/2 = 96 possible links, and since only three are present, p 2 = 1/32.At level j, an edge links nodes within a module containing 4 j nodes.The number of possible links at this level is 4 j (3 • 4 j−1 )/2, so that p j = 1/2 4j−3 , and, in general, p j /p j−1 = 1/16.Since a tree represents a minimally connected structure, p j cannot decay faster than this rate, in any connected hierarchical network based on four-node modules.

FIG. 1 :
FIG.1: Two lowest levels of the hierarchical network construction with 4 nodes/module.Dashed lines: l = l 1 , dotted lines: l = l 2 .The solid lines denoted R1 are randomly chosen connections within the bottom level (fully connected) modules, while those denoted R2 provide random connections on the next level.Links can be directed.

4 FIG. 2 :
FIG. 2: Number of nodes within the chemical distance r in the s = 3 HMN2d models with levels l = 9.Different curves correspond to different k /2-s.The dashed line shows a power-law fit for the b = 2, with N (r) ∼ r 1.4 , suggesting a topological dimension of D = 1.4 at the percolation threshold.

FIG. 3 :
FIG. 3: Number of nodes within chemical distance r in HMN2d networks with s = 4 and l = 9 levels.Different curves correspond to different k -s.Inset: local slopes d ef f of the N (r) curves, defined in Eq. 4.

FIG. 4 :FIG. 5 :
FIG.4: Number of nodes within chemical distance r in a set of ten random hierarchical trees with 262144 nodes.The broad yellow curve is an average over the set of ten structures.

FIG. 6 :
FIG.6: CP on asymmetric HMN2d networks with s = 4: decay of activity for λ values as indicated.

23 FIG. 16 :
FIG. 16: Evolution of the activity in the weighted CP on HMN2d networks with s = 4, l = 9, and k = 50 , for λ = 2.585, 2.588, and2.59(bottom to top).The main plot shows the decay in case of a fully active initial state, while the inset shows the growth of activity starting from a single active node for the same values of λ.Dashed lines are power-law fits for λ c ≃ 2.588.