Topological dimension tunes activity patterns in hierarchical modular network models

Connectivity patterns of relevance in neuroscience and systems biology can be encoded in hierarchical modular networks (HMNs). Moreover, recent studies highlight the role of hierarchical modular organization in shaping brain activity patterns, providing an excellent substrate to promote both the segregation and integration of neural information. Here we propose an extensive numerical analysis of the critical spreading rate (or"epidemic"threshold) --separating a phase with endemic persistent activity from one in which activity ceases-- on diverse HMNs. By employing analytical and computational techniques we determine the nature of such a threshold and scrutinize how it depends on general structural features of the underlying HMN. We critically discuss the extent to which current graph-spectral methods can be applied to predict the onset of spreading in HMNs, and we propose the network topological dimension as a relevant and unifying structural parameter, controlling the epidemic threshold.


INTRODUCTION
A surprising variety of biosystems display hierarchical modular architectures, an incomplete list of which includes the human brain structural network, or connectome [1,2], metabolic and regulatory networks [3,4] and fiber networks in connective tissue [5]. All such examples point to a common conjecture: hierarchical modular organization is capable of steering function and activity into localized patterns that provide an excellent tradeoff between segregation of functions (located at diverse moduli) and their integration across very diverse scales, constituting an elegant solution of the so-called segregation-integration dilemma [6][7][8][9][10][11][12]. Further examples of such remarkable dynamic signatures in fields beyond the biological realm have been recently pointed out for parallel processing [13], information retrieval [14], as well as random and quantum walks [15].
From the perspective of modeling techniques, hierarchical modular organization is often encoded in simple mathematical models of synthetic hierarchical modular networks (HMNs) [16]. Song et al. first noted that in HMNs the mutual distance between highly connected hubs is higher than in scale-free (SF) networks [17]. As a consequence, hubs and their neighborhoods are not easily exposed to functional overloads; even if load increases in some location, it mostly remains localized in a subnetwork and does not greatly affect activity at other hubs. Song et al. [17] proposed that the hierarchical (or fractal) organization of functional modules (e.g. metabolic networks) is an evolutionary constraint imposed by the need for network robustness. The concept of increased distances and path lengths in HMNs -as opposed to standard scale-free and small-world models-was later formalized through the observation that HMNs are endowed with a finite topological dimension, D, by Gallos et al. [18]. Let us recall that D quantifies how the number N r of nodes in the local neighborhood of an arbitrary node grows with the distance, r from it, N r ∼ r D ): lower values of D amount to higher distances between network hotspots, while, formally, networks with the small-world property have D → ∞.
The use of simple dynamical models have proven effective as a probing tool to understand the propagation of information or "activity" in biological networks. For instance, paramount features of brain activity have been first highlighted borrowing models from quantitative epidemiology [16,19]; in the context of activity spreading in cortical networks or connectomes -which are well represented by HMN architectures-it has been recently shown that simple "epidemic" processes of activity propagation such as susceptible-infected-susceptible (SIS) model are very convenient [19,20]. In SIS dynamics, nodes can be either active (infected I) or inactive (susceptible S); in terms of neural dynamics an active node corresponds to an active region in the brain, which can activate an inactive neighboring region with a given probability per unit time λ, and which are deactivated at rate µ (which we set to 1 without loss of generality) due to exhaustion of synaptic resources. In the standard critical point scenario, the steady-state average fraction of active nodes ρ ∞ (or prevalence) is zero for low λ and non-zero for high λ, these two regimes being separated by a critical value λ = λ c , at which scale invariant dynamic patterns are recorded. In HMNs instead, there is a whole range of λ values where scale-invariance is observed, i.e. a Griffiths phase [19,21]. The origin of such generic (i.e. occurring in an extended region in parameter space) scale-invariant behavior is believed to be mostly structural: hierarchical modular organization promotes the emergence of rare regions where activity tends to remain localized for large times even if finally it becomes extinct owing to fluctuations. Diagram of the construction method adopted to generate synthetic hierarchical modular networks (HMNs). Densely connected modules of size M0 represent the seed of the network structure (fully connected here). At each level i, modules of size Mi = 2 i M0 is formed as follows: pairs of modules of size Mi−1 are linked, by establishing random connections between their constituent nodes, with probability πi = α/4 i , where the connectivity strength, α, is a key parameter controlling the resulting architecture. This is just an example of a general mechanism by which structural disorder can induce Griffiths phases, with generic critical-like scale-invariant features in complex networks [22,23]. The structural origin of critical-like behavior is in agreement with renormalization arguments, which show how in hierarchical networks even standard percolation produces generic power-law distributions of connected component sizes, a feature otherwise normally ascribed to the critical point (or percolation threshold) [24,25].
Here, we propose the numerical study of the onset of spreading -the epidemic threshold-for the SIS in synthetic HMNs, and its relationship -employing spectral graph analyses-with structural parameters controlling the HMN architecture. We will highlight the topological dimension, D, as the relevant feature that can tune and control the value of the epidemic threshold λ c .

HIERARCHICAL MODULAR NETWORKS (HMNS)
We recall the general definition of a HMN, as a network in which smaller, more densely connected modules are clustered into larger and less densely connected super-modules, in an iterative fashion that spans several scales, or hierarchical levels. Diverse algorithms to generate synthetic HMNs have been proposed in the literature [16,19,26,27]. Here, we consider the model proposed in [19] motivated by previous works on optimal HMN architectures [28]. Networks of this type will be undirected and unweighted, although generalizations to the weighted and directed cases can be introduced straightforwardly.
A schematic description of this method is shown in Fig. 1. We call α the connectivity strength of a HMN, as it will appear as a chief parameter controlling the emerging topology. The network is organized in densely connected modules of size M 0 , which represent the level 0 of a hierarchy of links. At each hierarchical level i > 0, super-modules of size M i = 2 i M 0 are formed, joining sub-modules of size M i−1 by wiring their respective nodes with probability π i = α/4 i : the average number of links between two modules at level i will thus be n i = π i M 2 i−1 = α(M 0 /2) 2 , i.e. proportional to α, regardless of the value of i.
It was shown that this construction method ensures the scalability of the network structure, so that the average degree and the topological dimension D reach asymptotic values for large N [19]. In the N → ∞ limit, the effect of lowest-level modules becomes negligible (relegated to transient time scales) and the time asymptotics are dominated by the hierarchical organization: in this regime α becomes the only relevant construction parameter. We remark that, being n i the number of links between any two modules of size M i−1 , the maximum value of n i is given by M 2 0 , so that α can take values in the interval 4/M 2 0 ≤ α ≤ 4. Fig. 2 shows numerical measurements of the topological dimension D, which is found to increase linearly with α, in the limit of large system sizes N . While this result is only approximate for smaller values of N , where higher dimensions are complex to measure due to the size constraint, the linear behavior seems to take over for sizes around N ≈ 10 6 and above. This finding helps contextualize previous results, which pointed to a quasi-linear growth of D with α [19]. While that result was initially dismissed as a possible effect of a too-limited α window and a non-exhaustive size-scaling analysis, our current results, extending up to N = 2 24 ≈ 1.7 × 10 7 , seem to robustly confirm the conjecture of a D ∼ α dependence. While we have not been able to prove this result analytically so far, we will show in the rest of the paper its remarkable implications for activity spreading.

SPECTRAL ANALYSIS OF HMNS
The study of the epidemic threshold often relies on spectral arguments. In the case of SIS dynamics, network spectra refer to the adjacency matrix A, whose generic element A ij is 1 if nodes i and j are linked, and 0 otherwise. In particular, if one considers undirected networks as we do here, A is symmetric and thus diagonalizable. In most cases one is interested in the higher spectral edge of A, including the largest eigenvalues Λ 1 > Λ 2 ≥ Λ 3 ..., where the uniqueness of the largest Λ 1 is ensured by the Perron-Frobenius theorem whenever the network is connected. Provided that the spectral gap Λ 1 − Λ 2 is large, one can prove that λ QMF c = 1/Λ 1 , within the framework of the quenched mean-field (QMF) approximation, as we recall in the following.
The activity state for a SIS dynamic process is given by the column vector ρ, whose generic element ρ i is the probability that node i is in the active state at a generic time t, and its steady state limit ρ ∞ , obtained for t → ∞. Our starting point is the well-known exact result for the time evolution:ρ ≤ −ρ + λAρ, which originates from the full SIS problem by neglecting (negative) quadratic terms as well as dynamic correlation contributions [29]. In the t → ∞ limit (when the time derivative of ρ vanishes), using the eigendecomposition A which expresses the steady state ρ ∞ of a SIS process of spreading rate λ in terms of its projections on the eigenspaces of A, with the scalar product v m · ρ ∞ = v T m ρ ∞ . The prevalence or steady-state density of active nodes is then defined as ρ ∞ = |ρ ∞ |. If Λ 1 Λ 2 , the sum is dominated by the leading term, the scalar product is approximately equal to unity, and a non-trivial steady state may exist only if λ > 1/Λ 1 . In particular, the QMF result λ QMF c = 1/Λ 1 is recovered when the equal sign is taken in Eq. (1).
The generality of this result in SF networks has been the focus of a recent debate [30][31][32]. A crucial aspect here is introduced by the localization property. An eigenvector of A is called localized if it has all vanishing components, except for a small subset of them. A measure of eigenvector localization is provided by the inverse participation ratio [30,33,34]. It was initially proposed [30,31] that in networks displaying a localized v 1 the epidemic threshold should be replaced by an interval of the spreading rate λ, in which active states are short-lived, ideally a Griffiths phase [31], as unique unstable but localized modes do not suffice to create a global a b where inactive and active nodes are colored in green and red respectively. Patterns of activity are obtained through numerical simulations, as described in the main text. The network size is purposely kept small in order to make the figure clearer. Larger sizes are considered instead for the results described in the rest of the paper. (a) Activity for a metastable configuration that will eventually end in the absorbing ρ ∞ = 0 state, obtained for a subcritical value of λ < λc; activity is localized in some regions, but cannot be sustained dynamically for long times. (b) Activity for a configuration that corresponds to a non-zero steady state ρ ∞ , obtained for a value of λ > λc.
Here too activity appears localized and no spanning (i.e. percolating) cluster arises, however the stability of the steady state suggests that local activity patters interact dynamically -giving rise to an effective coalescence-resulting in a sustained global active state.
state of endemic network activity. While this view has been more recently challenged in the case of SF networks [32], we mentioned above that HMNs are now well-known to exhibit such a Griffiths phase, whose upper bound is given by the actual critical point λ c [19]. Such epidemic threshold does not comply with the QMF prediction, leading to λ c > λ QMF c = 1/Λ 1 [19]. It was argued [19] that this result is related to the fact that in HMNs: (i) spectral gaps are small; and (ii) localization extends to all eigenvectors in the higher spectral edge of A, that is, it is not limited to the principal eigenvector v 1 . We show here the extent to which the basic assumptions of QMF are altered in HMNs, thus making the estimate of λ c a formidable task.
The existence of a large spectral gap is equivalent to saying that the steady state is approximately given by the principal eigenvector v 1 : as argued above, the sum in Eq. (1) would be dominated by the m = 1 term and, taking the equal sign (QMF approximation), ρ ∞ would be simply given by v 1 . Such a strong statement works remarkably well in networks with strong small world properties [33] or entangled connectivity patterns [35,36], but fails to describe systems like ours in which spectral gaps are small [19,37]. We propose a way to relax the above assumption, by considering a candidate steady state σ ∞ as a linear combination of the first m * eigenvectors, that is In particular, we can always choose our eigenvector basis in such a way that v m · σ ∞ ≥ 0 for all m ≤ m * . We note that in the worst case scenario m * = N , meaning that all eigenvalues are needed. Under what conditions is the network able to sustain the candidate steady state, so that ρ ∞ = σ ∞ ? By inserting our expansion for the steady state in Eq. (1), and using the linear independence of eigenvectors, we obtain the condition Eq. (2) shows that in a spectral theory of SIS in which the first m * eigenvalues are necessary, also the corresponding eigenvectors play a role. Then, two mechanisms may be responsible for the lack of a non-trivial steady state with ρ ∞ = 0: (a) if the spreading rate is small, i.e. λ < 1/Λ 1 (and thus, λ < 1/Λ m for all m) then the candidate non-trivial steady state σ ∞ is unstable; (b) even if λ ≥ 1/Λ m for some m, the localization of v m induces v m · σ ∞ ≈ 0 (positive but vanishing): the candidate steady state is weakly stable in principle, however dynamic fluctuations are bound to deactivate it over time. In this case, a virtually stable state results in fact metastable and relatively short-lived [39], in agreement with the Griffiths phase picture [19]. Thus, a nonzero steady state is still possible, if one of the following two conditions is met: i) there are some de-localized eigenstates, i.e. a pathological region in the spectrum exists around a given Λ m , characterized by high density of states and/or low localization, such that it can trigger a stable ρ ∞ such that v m · ρ ∞ 0 and λ c ≈ 1/Λ m ; ii) a finite number of localized states can dynamically coalesce sustaining a global active state (see Fig. 3), i.e. no pathological (Λ m , v w ) pairs exist, however if λ is large enough, λΛ m − 1 > 0 for many values of m < m * : a large enough number of potentially short-lived states is generated, large enough to sustain (or coalesce into) a long-lived steady state ρ ∞ > 0, however with no clear criterion to identify what a large enough m * is.
In order to clarify which of the two scenarios holds in the case of HMNs, and how the value of λ c can be related to spectral properties of A, we conducted an extensive numerical study of HMNs of different α, targeting both the SIS dynamics direct simulation and the network spectral properties.
Our numerical results provide a clear indication that the relevant scenario for HMNs is indeed ii) as exemplified graphically in Fig. 3, illustrating the lack of any de-localized state [40]. Similarly, we did not find any correlation between the inverse of any eigenvector and the numerically determined value of λ c (see below), thus indicating that if an epidemic threshold is reached -above which sustained activity exists-it has to emerge owing to a finite number of unstable, though localized, active regions and their dynamical interplay.

SCALING OF THE EPIDEMIC THRESHOLD IN HMNS
Given the complexity of this scenario, in which a large but undefined number of eigenvalue-eigenvector pairs is expected to play a role; is it still possible to relate λ c to a single scalar property of the spectrum?
In order to answer this question, we notice that the properties of the higher spectral edge of HMNs can be tuned by acting on α: upon increasing α, the number of non-zero off-diagonal elements of A increases accordingly. A simple argument to shed light on this observation is as follows: at every hierarchical level i, the approximate and coarse-grained connectivity pattern between two modules would be very roughly represented by an effective weighted adjacency matrix structured as A (i) ∼ c i 0 α α 0 , which will contribute a larger eigenvalue proportional to c i α to the higher spectral edge of A. We can corroborate our view by looking at Fig. 4a, where it is shown that the value of the generic i-th eigenvalue is strictly proportional to α, for any value of i in the higher spectral edge. This is unlike what happens in Erdös-Rényi (ER) or SF networks, where increasing connectivity -for instance by tuning the average degree in ER networks or the degree distribution in SF networks -primarily affects the scaling of the principal eigenvalue Λ 1 only, effectively changing the size of the spectral gap [33]. We stress that our approach here is kept purposely simple, while the correct way to address the analytical computation of HMN spectra has been recently discussed for the Laplacian spectrum of the Dyson fully connected hierarchical graph [14,15]. The exceptional behavior of HMN spectra, for which a single tuning parameter is able to tune the entire spectral edge, turns crucial for our study of the epidemic threshold: provided that λ c scales as the inverse of some relevant eigenvalues in the higher spectral edge, the similarity properties of eigenvalue spectra Λ m ∼ α provides us with a simple prediction for the epidemic threshold, which can be expressed as: We measured λ c computationally by employing the standard method of simulating SIS dynamics starting from a homogeneous initial condition of all active nodes, and recording steady states values (more refined techniques relying on susceptibility measurements have been proposed for SF networks, in which vanishing values of λ c make their estimate a more delicate task [38]). Simulations we performed for sizes N = 2 17 and N = 2 20 , recording minimal to no deviations between the two cases. From the perspective of size scaling, it appears that in this size range (10 5 -10 6 nodes) the dynamics display no appreciable size effects. In particular, we find that steady states for λ > λ c are stable against fluctuations, and do not decay as a consequence of finite sizes. Smaller sizes, however, may still allow large enough fluctuations to make steady states unstable. This apparent inconsistency (an unclear distinction between supercritical and subcritical dynamics), affects systems of very limited sizes, and is of course understandable considering that phase transitions are correctly defined only in the N → ∞ limit. The resulting λ c is plotted in the inset of Fig. 4b as a function of α. The validity of Eq. (3) is confirmed by simulation results in Fig. 4b, where it is shown that there is a remarkable scaling property for all eigenvalues in the spectral edge. In passing, we note that as all eigenvalues in the range scale with 1/λ c , this must happen, in particular, for Λ 1 : while the principal eigenvalue alone by itself cannot justify the value of the epidemic threshold, the standard λ QMF c = 1/Λ 1 criterion survives in a weaker form as a scaling law (not an equality), and more importantly, it extends to all eigenvalues that participate in the onset of activity.
Finally, recalling our initial results, that for large enough networks one has D ∼ α, we can rewrite our estimate as The conjecture that the network dimension acts as a structural determinant of spreading activity appears corroborated: D is indeed capable of tuning the value of λ c . The picture emerging from our results is that of a complex dynamic scenario, as the one described by Eq. (2), which can be correctly understood by identifying the topological dimension D as the relevant tuning parameter for epidemic spreading in HMNs.
We also remark that in the D → ∞ limit, Eq. (3) predicts a vanishing epidemic threshold, recovering the wellknown quenched mean field result for SF networks, which are intrinsically infinite dimensional. The actual study of the existence of a finite λ c in SF graphs requires in fact more refined theoretical tools, while here we only stress how our generalization of the QMF approach also contains its original predictions.
Finally, we addressed the question of how the supercritical phase diagram is affected by tuning D, or equivalently α. Fig. 5 shows the steady state density of active nodes, ρ ∞ , as a function of (λ − λ c )/λ c , for different values of α (and D). The striking result that we find is that while the emergence of scaling clearly suggests a critical scenario of the familiar form ρ ∞ ∼ (λ − λ c ) β , even for the lower values of α considered here the system is above its critical dimension and the β exponent is insensitive to dimensionality. A detailed study of the phase transition for even lower dimensional HMNs is beyond the scope of this manuscript, and is being considered for future work. Here, we would like to emphasize that the emergence of dynamic scaling is robust against structural variations, in a range of connection densities and dimensions which we consider relevant in fields such as neuroscience.

CONCLUSIONS
We discussed the essential reasons why the QMF approximation does not provide a correct prediction for the epidemic threshold in hierarchical modular networks. The first eigen-mode to become unstable is localized and thus, it cannot represent an endemic state of activity once fluctuations are taken into consideration. Instead, a finite number of unstable eigen-modes may induce an active state. Remarkably such set of eigenvectors does not create a spanning or percolating cluster expanding through the whole network, but just localized clusters of activity. However such "reservoirs" of instability can create a global active state owing to fluctuations that transiently activate other modes, effectively connecting diverse localized clusters.
Consequently, we introduced a novel framework, in which the epidemic threshold depends on a (large) number of eigenvalues but, remarkably, turns out to be inversely proportional to a unique parameter; the network topological dimension, D. While our results corroborate this result for HMNs, the possibility of extending it to more general network classes can only be conjectured at this stage.
The importance of our results in the description of biological systems endowed with hierarchical modular organization can be better understood considering the example of brain networks. A direct dependence of λ c on D ensures that, being HMNs finite dimensional, the epidemic threshold never vanishes. Should that happen, information units would boundlessly propagate through the network, making it impossible to manage information properly, something usually associated with epileptic forms of activity. On the other hand, by simply tuning a single parameter, e.g. during the development or pruning of neuronal networks, the dynamical regime of the whole network could be established.
In conclusion, we have studied numerically some of the most relevant properties of the onset of spreading in HMNs. We have highlighted a crucial construction parameter, the connectivity strength α -that controls the network fractal-like structure and shown its proportionality to the topological dimension of a HMN. We have shown how α and, equivalently-D are responsible for the tuning of the epidemic threshold λ c in HMNs, by acting directly on the spectral properties of the adjacency matrix. Thus, by slightly modifying a unique network-structure controlling parameter, it is possible to regulate the spreading rate that is necessary to generate sustained activity, i.e. the "epidemic threshold" and, in this way, the overall state of activity can be controlled by the network architecture.
We hope that our work can stimulate further studies in this direction, possibly providing a deeper analytical understanding of the relationship and interplay between structural and dynamic patterns of localization.
AS and PM acknowledge financial support from the Deutsche Forschungsgemeinschaft, under grant MO 3049/1-1. MAM is grateful to the Spanish-MINECO for financial support (under grant FIS2013-43201-P; FEDER funds). We thank P. Villegas and S. di Santo for useful discussions and for a critical reading of the manuscript.