key: cord-0010890-1u3hvjrt
authors: Li, Shanshan
title: Griffiths phase on hierarchical modular networks with small-world edges
date: 2017-03-06
journal: Phys Rev E
DOI: 10.1103/physreve.95.032306
sha: 76f2ed0058a80e202a3a8031141bbdb2581fcfec
doc_id: 10890
cord_uid: 1u3hvjrt

The Griffiths phase has been proposed to induce a stretched critical regime that facilitates self-organizing of brain networks for optimal function. This phase stems from the intrinsic structural heterogeneity of brain networks, i.e., the hierarchical modular structure. In this work, the Griffiths phase is studied in modified hierarchical networks with small-world connections based on the 3-regular Hanoi network. Through extensive simulations, the hierarchical level-dependent inter-module wiring probabilities are identified to determine the emergence of the Griffiths phase. Numerical results and the complementary spectral analysis of the relevant networks can be helpful for a deeper understanding of the essential structural characteristics of finite-dimensional networks to support the Griffiths phase.

The well-known criticality hypothesis suggests that biological systems operate at the borderline between the sustained active and the sustained inactive state. This has been observed in various processes such as gene expression [1] , cell growth [2] , and neuronal avalanches [3] . In different contexts, the critical point enables optimal transmission and storage of information [4, 5] , maximal sensitivity to stimuli [6] , and optimal computational capabilities [7] . Empirical studies on brain networks [8] [9] [10] , however, exhibit a stretched critical region. The Griffiths phase, characterized by generic power laws over a broad region in the parameter space, provides an alternative mechanism for critical behavior in brain networks without fine tuning [11, 12] . It is confirmed numerically and analytically that the structural heterogeneity induces the Griffiths phase, which eventually enhances the self-organization mechanism of brain networks.

Brain networks have been found to be organized into modules across hierarchies [13] [14] [15] . Modules in each hierarchy are grouped into larger modules, forming a fractallike structure. Previous work models brain networks with finite-dimensional hierarchical modular networks (HMNs) [11, 12] and successfully confirms the existence of the Griffiths phase using dynamical models, such as the susceptibleinfected-susceptible (SIS) model and the contact process. The essential characteristic of previous network models is the hierarchical level-dependent intermodule wiring probabilities. It is conjectured that plain modular networks are not able to support the Griffiths phase, and disorder in different scales significantly influences the properties of critical behaviors [11] . In this work, the idea of a Griffiths phase is extended to other hierarchical structures encountered in previous studies of dynamical processes on complex networks.

Certain hierarchical networks, with a self-similar structure and small-world connections, have been shown to exhibit novel dynamics [16] [17] [18] [19] [20] [21] . In this work, the hierarchical models are designed based on one such example with a finite topological dimension, the 3-regular Hanoi network [16, 18, 22, 23] . To tune the modular feature that is present in brain networks, * shanshan.li@emory.edu a single node of the original network is modified into a fully connected clique with a varying size. By introducing different classes of intermodule connections, the essential heterogeneous connectivity pattern is explored to induce the Griffiths phase on finite-dimensional networks. It turns out that the hierarchical level-dependent intermodule wiring probabilities play an important role, affecting the properties of the phase transition at criticality.

As a complement to the computational approach, the spectral analysis of the adjacency matrix and the Laplacian matrix of networks is conducted. A localized principal eigenvector of the network adjacency matrix indicates the network heterogeneity, which has been used to quantify the localization of activity on networks above the critical propagation rate in the dynamical model [24] . This concept has been applied to analytically explain the emergence of rare regions and the Griffiths phase [11, 12, 25] . However, the observation that a localized principal eigenvector is not necessarily the fingerprint of the Griffiths phase has been reported in highly connected networks with intrinsic weight disorder or finite-size random networks with power-law degree distributions [12, 26] . As an extension to finite-dimensional models, a class of networks is found here where the Griffiths phase is absent although their principal eigenvectors are localized. As the second approach of the spectral analysis of the network connectivity matrix, the Lifschitz tails in the spectral density of the Laplacian matrix is proposed to predict the Griffiths phase analytically. Lifshitz tails have been related to the Griffiths singularity in equilibrium systems [27] . For synchronization and spreading dynamics on networks, simulation and quenched mean-field (QMF) approximation indicate a connection between the Lifshitz tail and the slow dynamics [11, 28, 29] . In this study, the tail distribution of the Laplacian eigenvalues is presented to test how well it predicts the Griffiths phase in the SIS model. This paper is organized as follows: Section II describes the structural properties of hierarchical modular networks, in which the SIS model and its critical behavior are studied; Sec. III reviews the SIS model and the spectral analysis of the adjacency matrix and the Laplacian matrix, and the analytical tools are applied to all the networks considered in this work; Sec. IV presents the numerical results for the SIS model evolving in the networks; Sec. V concludes by highlighting The network features a regular geometric structure, in the form of a onedimensional backbone, and a distinct set of recursive small-world links. The node degree is uniformly 3. the significance of hierarchical level-dependent intermodule wiring probabilities for the emergence of the Griffiths phase and discussing the essential structural characteristics of the Griffiths phase.

Hanoi networks [16, 18, 22, 23] are based on a simple geometric backbone, a one-dimensional line of n = 2 g nodes. Each node is at least connected to its left and right nearest neighbors on the backbone. To construct the hierarchy to the gth generation, consider parameterizing any node x < n (except for 0) uniquely in terms of two integers (i,j ), i 1 and 1 j 2 g−i , via

Here, i denotes the level of hierarchy, whereas j labels consecutive nodes within each hierarchy. This parametrization raises a natural pattern for long-range small-world edges that are formed by the neighbors x = 2 i−1 (4j − 3) and y = 2 i−1 (4j − 1) for 1 j 2 g−i−1 , as shown in Fig. 1 . Eventually, this procedure constructs a finite-dimensional hierarchical network with a uniform finite node degree 3 and a diameter of ∼ √ n, which is denoted HN3 [16, 18, 22] . To construct a hierarchical structure that models the modular property of real-world brain networks, each single node x in HN3 is replaced by a fully connected clique that contains a finite number m of nodes, thus forming a network of size n × m. Maintaining the structural properties of HN3, the self-similar structure, and the small-world connections, I propose two connectivity patterns between modules in the same hierarchy. In the first paradigm, the single edge in the original HN3 is now formed by two randomly chosen nodes in different cliques, which I denote HMN1. The second paradigm is inspired by previous hierarchical modular models [11, 12] . To distinguish it from HMN1, I denote it HMN2. Previous models share common features, hierarchical construction of modules and level-dependent wiring probabilities. Modules are grouped to form larger modules on the next level. They are connected in either a stochastic way with a level-dependent probability p i or a deterministic way with a level-dependent number of edges.

Since an infinite-dimensional network is predicted not to support the Griffiths phase [11] , to maintain a finite fractal dimension, the number of intermodule connections is stable across hierarchical levels. In this work, I use the stochastic scheme to construct HMN2. As the size of modules in the hierarchical level i increases exponentially with i, the intermodule wiring probability decreases exponentially. In HMN2, for the first hierarchy, modules on this level are cliques themselves. Starting from the second hierarchy, the clique labeled 2(2j − 1) is grouped with the neighbor clique 2(2j − 1) − 1 and 2(2j − 1) + 1, forming a module. For the third hierarchy, the clique labeled 2 2 (2j − 1) is grouped with three left neighbor cliques up to the clique 2 2 (2j − 1) − 3 and three right neighbor cliques up to the clique 2 2 (2j − 1) + 3. Repeating this procedure g generations, the size of the module at the ith generation is m(2 i − 1). The number of all possible stochastic connections between two modules is m 2 (4 i − 2 i+1 + 1). Thus, to ensure that at least one edge exists between them, the level-dependent probability p i is bounded by 1/(m 2 (4 i − 2 i+1 + 1)).

Certain fundamental dynamical models, the susceptibleinfected-susceptible model and the contact process, have been used to model the activity propagation in brain networks [11, 12] . Previous studies focus on the emergence of the Griffiths phase in general complex networks using these simplified models. Quenched disorder, either intrinsic to nodes or topological, has been shown to smear the phase transition at critical points and generate the Griffiths phase. The essential disorder may stem from a node-dependent propagation rate [30, 31] . Recent results also present evidence that the Griffiths phase emerges due to the quenched disorder on the edges, such as a correlated weight pattern in tree networks [32] and an exponentially suppressed weight scheme in random networks [25] .

Special rare regions emerge in the dynamical process evolving in networks with quenched disorder. Statistically, the active state lingers in these rare regions for a typical time that grows exponentially with their sizes and eventually ends up in the absorbing state [11, 33] . The exponential size distribution of rare regions induces power-law decays with continuously varying exponents, i.e., the Griffiths phase. Not only in the spreading dynamics, rare regions have also been shown to affect the properties in classical phase transitions in quenched disordered systems, such as the randomly diluted Ising model and Ising model with planar defects, and quantum phase transitions in itinerant magnets with Heisenberg spin symmetry, leading to an essential singularity, the Griffiths singularity [34, 35] .

In the absence of quenched disorder, the Griffiths phase can also be a consequence of the structural heterogeneity of finite-dimensional networks, which is expected to play a role similar to that of the quenched disorder [11, 30] . The QMF approximation applies a spectral analysis of the network adjacency matrix that analytically explains emerging rare regions and the Griffiths phase in networks with quenched disorder [24, 25] . This analytical procedure successfully confirms the Griffiths phase in finite-dimensional hierarchical modular networks in previous work [11] . The spectral analysis of the network Laplacian matrix provides another approach to predicting the Griffiths phase, which is confirmed in [11] and [29] . In this section, I focus on the SIS model and apply the spectral analysis of all the relevant finite-dimensional structures.

In the SIS model, each node in networks is described by a binary state, active (σ = 1) or inactive (σ = 0). An active node is deactivated at a unit rate; otherwise, it propagates the active state to its inactive neighbors at a rate λ. The evolution equation for the probability ρ x (t) that node x is active at time t is

in which A is the network adjacency matrix. A xy is 1 if nodes x and y are connected by an edge; otherwise, it is 0. The Laplacian matrix of a graph is defined with the adjacency matrix as

where L xy is equal to −A xy when x = y, and L xx is equal to y =x A xy , i.e., the degree of node x. Denote eigenvalues and eigenvectors of the adjacency matrix A and the Laplacian matrix L, respectively, as A , f A ( A ) and L , f L ( L ). I here briefly introduce the method used to perform the simulation for the SIS model. The large-scale numerical simulation method for the SIS model developed in [36] determines the critical propagation rate λ c efficiently for various networks. This algorithm considers the SIS model in continuous time. At each time step, one randomly chosen active node deactivates with the probability N i /(N i + λN n ), where N i is the number of active nodes at time t and N n is the number of edges emanating from them. With complementary probability λN n /(N i + λN n ), the active state is transmitted to one inactive neighbor of the randomly selected node. Time is incremented by t = 1/(N i + λN n ). This process is iterated after updating the system.

In this subsection, I review the derivation of the spectral analysis of the adjacency matrix and the Laplacian matrix and apply it to all the relevant networks. Considering the adjacency matrix, the criterion for localization of the steady active state in networks is based on the evaluation of the inverse participation ratio (IPR) of the principal eigenvector corresponding to the largest eigenvalue. Following the notations in [24] , eigenvalues of the adjacency matrix are ordered as A

The probabilities ρ x for each node at the steady state can be written as a linear superposition of the N orthogonal eigenvectors [24] ,

If the largest eigenvalue A 1 is significantly larger than all the others, the QMF approximation predicts the critical point λ c as 1/ A 1 , and the steady state probability as

The order parameter ρ(t) is defined as the average 1 N N x=1 ρ x (t) over all the nodes. At the critical λ c , the order parameter ρ at the steady state can be expanded as,

in which = λ 1 − 1 1 with the coefficients

With the dominant largest eigenvalue and the principal eigenvector, the order parameter ρ can be approximated with ρ ∼ a 1 

, and then a 1 ∼ const and ρ ∼ const. The active state extends over a finite fraction of nodes of the network. The localization of eigenvectors is quantified by their inverse participation ratio (IPR) (shown in [24] ),

A finite IPR value of the principal eigenvector corresponds to a localized eigenvector, while an IPR approaching 0 corresponds to a delocalized principal eigenvector. I apply the concept of the IPR to all the relevant networks to examine whether a localized principal eigenvector exists, which may suggest the emergence of rare regions and the Griffiths phase in the QMF approximation [11, 12] . As shown in Fig. 2(a) , the IPR of the principal eigenvector increases with the clique size m towards a finite value for different maximum generation g of HMN1. This suggests that the principal eigenvectors for HMN1 configurations are localized. Additionally, localized eigenvectors corresponding to the largest eigenvalues are also found, shown in Fig. 2(b) . For HMN2, I focus on level-dependent intermodule wiring probabilities, p i = 4 −(i+1) and p i = 4 −i , in this work. The backbone as well as the first-level intermodule wiring probability is fixed at 1/4, where the modules are the basic cliques described in Sec. II. Values of IPR are shown in Fig. 3(a) . The largest value comes from the network configuration where the single clique contains two nodes, and the probability is p i = 4 −(i+1) . In this case, the network is statistically almost fragmented. Numerical results in Sec. IV indeed show the emergence of the Griffiths phase as a trivial consequence of the network disconnectedness. To examine the Griffiths phase in a connected network with a finite fractal dimension, the network configuration m = 3 and p i = 4 −(i+1) is also chosen for the numerical simulation. For the stochastically constructed HMN2, as the clique size m or level-dependent probability increases while the other factor is kept fixed, modules become more and more connected with other modules in the same level, and the value of IPR decreases, shown in Fig. 3(a) . The regime over the parameter m and the level-dependent p i that possibly supports the Griffiths phase is narrow. It is not surprising to see that the localized principal eigenvector exists for HMN2 network configurations with a finite value of IPR. In Figs. 3(b) and 3(c), I illustrate the localized eigenvectors corresponding to large eigenvalues in two network configurations.

As the second approach of the spectral analysis to further confirm the Griffiths phase, I study the spectral density of the network Laplacian at the lower edge, the Lifshitz tail. The Laplacian matrix is positive semidefinite, i.e., L i 0 and L 1 = 0, following the notations in [29] . The smallest nonzero 

A xy ρ y (t), (9) which can be rewritten using the Laplacian matrix as

A linear stability analysis is performed [29] , similar to the synchronization process [37] . The normal modes of the perturbations above the absorbing state can be written

Using the Laplacian spectrum, ρ x (t) can be expressed as the expansion with the Laplacian eigenvalues and eigenvectors,

and the total density is determined by the lowest eigenvalues of the spectrum,

In the continuum limit, in which L c is the experimentally determined end of the tail value. A power-law distribution P ( ) of the lower edge of the Laplacian spectrum suggests Griffiths phase behavior above the absorbing state [29] :

For comparison I calculated the Lifshitz tails for HMN1 and HMN2, shown in Fig. 4 . In the plot, the probability distribution P ( ) is calculated with the bin size δ = 0.0001 over 100 independent graph realizations. The Lifshitz tail for the HMN2 configuration with g = 13, m = 2, p i = 4 −(i+1) and with (17) and

while the Lifshitz tail for HMN1 slightly deviates from a power law, suggesting the lack of the Griffiths phase according to [29] .

In this section, I present results from numerical study of the SIS model in all the network configurations of HMN1 and HMN2 using the simulation method introduced in Sec. III A. The network is initialized as a fully active graph. The system is updated each step until the maximum time t max (10 6 ) is reached or in the case of activity extinction. Simulations for each propagation rate λ are repeated on 1000-5000 independent network realizations. The order parameter ρ(t) for each λ is the average of all the runs. I also derive the effective decay exponent by fitting critical power laws ρ(t) ∼ t −α eff with the efficient exponent defined as [12, 25] Fig. 7 . However, the Griffiths phase is a trivial consequence of the disconnectedness of HMN2 as discussed in Sec. III B. For network configurations of HMN2 that are connected, I choose the case of m = 3, p i = 4 −(i+1) , at which the corresponding value of IPR is sufficiently large. Since intermodule connections are established stochastically, there is a chance that all the possible intermodule edges fail to be connected. To avoid this case, at least one intermodule connection is enforced to exist by repeating the construction of graphs in the simulation. Numerical results for a connected HMN2 are presented in Fig. 8 . We see a nearly sizeindependent power law in a stretched regime of λ. Comparing Fig. 7 with Fig. 8 , as m increases while p i is kept fixed, and vice versa, the regime in the parameter space of λ for the Griffiths phase is expected to become narrow until it disappears when HMN2 becomes highly connected.

As introduced in Sec. I, one significant advantage of biological systems operating at criticality is the diverging reaction to highly diverse stimuli. From the perspective of statistical mechanics, this is caused by the divergence of susceptibility at criticality. To measure the divergence of response in the Griffiths phase, here I use the concept of dynamic susceptibility, which is applied to gauge the overall response to a continuous localized stimulus in [11] . This dynamic susceptibility is defined as

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. 9 , becomes large in the region of the Griffiths phase, and more importantly it grows with the network size N , which implies a divergent response over an extended region.

The Griffiths phase induced purely by structural disorder suggests an alternative self-organizing mechanism for brain function. Brain networks are shown to have strong modularity, and densely connected modules are organized in a hierarchical pattern. Another important feature found empirically is the small-world topology, which ensures the efficient information transfer between modules. To solve the the conundrum between small-world topology and large-world architecture of brain networks, incorporating progressively weaker intermodule connections while maintaining well-defined modules is proposed as a solution [38] . In this work, I construct two classes of synthetic hierarchical modular networks that incorporate weak intermodule connections in distinguished ways. The first model, HMN1, builds the hierarchical smallworld connections into planar modular networks. The second model, HMN2, organizes the modules into hierarchies, and the intermodule wiring probability is level dependent. Both of them possess a self-similar structure and small-world long-range connections, based on the 3-regular Hanoi network [23] .

I study the Griffiths phase by evolving the fundamental SIS model in the HMNs designed. As a further exploration into the Griffiths phase caused by the structural heterogeneity of networks, I present numerical results for two classes of networks. The results suggest the essential role of the leveldependent intermodule wiring probability in the emergence of the Griffiths phase. The first class of hierarchical networks, HMN1, is not able to support the Griffiths phase, although they satisfy the structural criteria, such as the finite fractal dimension, the modular structure, and the hierarchical heterogeneity. The second class of hierarchical networks, HMN2, is constructed to possess a hierarchical pattern in the intermodule wiring probabilities, which therefore require a delicate tuning to maintain a connected, finite-dimensional network. The hierarchical pattern of intermodule connections results in more heterogeneous HMN2 networks, while HMN1 network configurations are more homogeneous. This can be shown by considering the Lifshitz tails for HMN1 and HMN2 in Fig. 4 . The Lifshitz tail of HMN1 deviates from that of HMN2, and compared to HMn2, the spectrum of HMN1 is closer to the original network HN3 [39] . Although HN3 possesses a hierarchical structure, it is not sufficiently disordered to induce the Griffiths phase. The difference in the hierarchical patterns of HMN1 versus HMN2 significantly affects the phase transition and existence of the Griffiths phase.

As a complement to the computational efforts, the spectral analysis proposed in the quenched mean-field approximation suggests that a finite IPR of the principal eigenvector of the adjacency matrix can be considered to indicate the localization of activity that may result in the emergence of rare regions and the Griffiths phase under certain circumstances. Although all the network configurations of HMN1 prove to have a finite IPR and localized eigenvectors corresponding to the largest eigenvalues, only when the structural disorder of intermodule connections is sufficient as in HMN2 does the Griffiths phase appear. As a counter example to previous finite-dimensional models with localized principal eigenvectors that support the Griffiths phase [11, 12] , a class of finite-dimensional networks with a localized principal eigenvector is found where the Griffiths phase is absent. This raises questions about a more generalized theoretical analysis of the network adjacency matrix that applies to all networks considered previously and currently. Besides the spectral analysis of the adjacency matrix, the Lifshitz tail of the Laplacian spectrum presents a power-law probability distribution at the lower edge of the spectrum of HMN2 networks, while the tail distribution in HMN1 deviates from a power law. Numerical results confirm that the property of the phase transition may be related to this difference in the Lifshitz tails of HMN1 versus HMN2.

Proc. Natl. Acad. Sci. USA 105

Proc. Natl. Acad. Sci. USA

I would like to thank Prof. Stefan Boettcher for helpful discussions. This work was supported by the NSF through Grant No. DMR-1207431, which is gratefully acknowledged.