Topological dimension tunes activity patterns in hierarchical modular networks

A surprising variety of biosystems display hierarchical modular architectures. An incomplete list of these 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 overarching principle: 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–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 nodes is higher than in scale-free (SF) networks [17]. As a consequence, highly connected nodes 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. in 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 SF 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}\sim {r}^{D}$): lower values of D amount to higher distances between network hotspots, while, formally, networks with the small-world property have $D\to \infty$.

The use of simple dynamical models has 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 structural cortical networks (connectomes)—which are well represented by HMN architectures—it has been recently shown that simple 'epidemic' processes of activity propagation such as the susceptible-infected-susceptible (SIS) model can be 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 rate λ, and which can be deactivated at rate μ (which we set to 1 without loss of generality) due to exhaustion of synaptic resources. In the standard scenario, the steady-state average fraction of active nodes ${\rho }^{\infty }$ (or prevalence) is zero for low λ and non-zero for high λ; at the edge of these two regimes there is a critical value $\lambda ={\lambda }_{{\rm{c}}}$, at which scale invariant dynamic patterns are recorded. Instead, in HMNs there is a whole range of λ values where scale invariance is observed, i.e. what is called 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 mostly localized for large times, even if it finally becomes extinct owing to fluctuations.

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–24]. In the case of HMNs, structural disorder arises from the stochastic nature of the processes by which modules at diverse scales are connected. 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) [25, 26].

Here, we propose the numerical study of the onset of spreading—the epidemic threshold—for SIS dynamics in synthetic HMNs, and—employing spectral graph analyses—we investigate its relationship with structural parameters controlling the HMN architecture. We will highlight the network topological dimension, D, as the relevant feature that can tune and control the value of the epidemic threshold ${\lambda }_{{\rm{c}}}$.

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, 27, 28]. Here, we consider the model proposed in [19] motivated by previous works on optimal HMN architectures [29]. 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 figure 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₀, which represent the level 0 of a hierarchy of links. At each hierarchical level $i\gt 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 ${\pi }_{i}=\alpha /{4}^{i}$: the average number of links between two modules at level i will thus be ${n}_{i}={\pi }_{i}{M}_{i-1}^{2}=\alpha {({M}_{0}/2)}^{2}$, i.e. proportional to α, regardless of the value of i.

Zoom In Zoom Out Reset image size

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]. For a generic dynamic process running on a HMN in the $N\to \infty$ 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}_{0}^{2}$, so that α can take values in the interval $4/{M}_{0}^{2}\leqslant \alpha \leqslant 4$. Figure 2 shows computational 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⁶ 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}\approx 1.7\times {10}^{7}$, seem to robustly confirm the conjecture of a $D\propto \alpha$ 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.

Zoom In Zoom Out Reset image size

The study of the epidemic threshold often relies on spectral arguments. In the case of SIS dynamics, network spectra refer to the adjacency matrix ${\bf{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, ${\bf{A}}$ is symmetric and thus diagonalizable. In most cases one is interested in the higher spectral edge of ${\bf{A}}$, including the largest eigenvalues ${{\rm{\Lambda }}}_{1}\gt {{\rm{\Lambda }}}_{2}\geqslant {{\rm{\Lambda }}}_{3}...$, where the uniqueness of the largest ${{\rm{\Lambda }}}_{1}$ is ensured by the Perron–Frobenius theorem whenever the network is connected. Provided that the spectral gap ${{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2}$ is large, one can prove that ${\lambda }_{{\rm{c}}}^{\mathrm{QMF}}=1/{{\rm{\Lambda }}}_{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 ${\boldsymbol{\rho }}$, whose generic element ${\rho }_{i}$ is the probability that node i is in the active state at a generic time t; we call ${{\boldsymbol{\rho }}}^{\infty }$ its steady state limit, obtained for $t\to \infty$. Our starting point is the well-known exact result for the time evolution: $\dot{{\boldsymbol{\rho }}}\leqslant -{\boldsymbol{\rho }}+\lambda {\bf{A}}{\boldsymbol{\rho }}$, which originates from the full SIS problem by neglecting (negative) quadratic terms as well as dynamic correlation contributions [30]. In the $t\to \infty$ limit (when the time derivative of ${\boldsymbol{\rho }}$ vanishes), using the eigendecomposition ${\bf{A}}={\sum }_{m=1}^{N}{{\rm{\Lambda }}}_{m}{{\bf{v}}}_{m}^{T}{{\bf{v}}}_{m}$, where ${{\bf{v}}}_{m}$ is the mth eigenvector, we obtain

$$\begin{eqnarray}&&{{\boldsymbol{\rho }}}^{\infty }\leqslant \lambda \displaystyle \sum _{m=1}^{N}{{\rm{\Lambda }}}_{m}({{\bf{v}}}_{m}\cdot {{\boldsymbol{\rho }}}^{\infty }){{\bf{v}}}_{m},\end{eqnarray} \tag{ 1 }$$

which expresses the steady state ${{\boldsymbol{\rho }}}^{\infty }$ of a SIS process of spreading rate λ in terms of its projections on the eigenspaces of ${\bf{A}}$, with the scalar product ${{\bf{v}}}_{m}\cdot {{\boldsymbol{\rho }}}^{\infty }={{\bf{v}}}_{m}^{T}{{\boldsymbol{\rho }}}^{\infty }$. The prevalence or steady-state density of active nodes is then defined as ${\rho }^{\infty }=| {{\boldsymbol{\rho }}}^{\infty }|$. If ${{\rm{\Lambda }}}_{1}\gg {{\rm{\Lambda }}}_{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 $\lambda \gt 1/{{\rm{\Lambda }}}_{1}$. In particular, the QMF result ${\lambda }_{{\rm{c}}}^{\mathrm{QMF}}=1/{{\rm{\Lambda }}}_{1}$ is recovered when the equal sign is taken in equation (1).

The generality of this result in SF networks has been the focus of a recent debate [31–33]. A crucial aspect here is introduced by the localization property. An eigenvector of ${\bf{A}}$ is called localized if it has all vanishing components, except for a small subset of them. A quantitative measure of eigenvector localization is provided by the inverse participation ratio (IPR) $I({{\rm{\Lambda }}}_{m})={\sum }_{i=1}^{N}{v}_{{mi}}^{4}$ of eigenvector ${{\bf{v}}}_{m}$, corresponding to eigenvalue ${{\rm{\Lambda }}}_{m}$ [31, 34, 35]. A perfectly delocalized eigenvector, with components of the order of ${ \mathcal O }({N}^{-\tfrac{1}{2}})$, has $I({\rm{\Lambda }})\sim { \mathcal O }({N}^{-1})$, while for a localized eigenvector $I({{\rm{\Lambda }}}_{m})$ remains constant for $N\to \infty$. It was initially proposed [31, 32] that in networks displaying a localized ${{\bf{v}}}_{1}$ the epidemic threshold should be replaced by an interval ${ \mathcal R }$ of spreading rate values λ, in which active states are short-lived, ideally a Griffiths phase [32], as unique unstable but localized modes do not suffice to create a global state of endemic network activity. In such a situation, a true active state would be achieved at a ${\lambda }_{{\rm{c}}}$ equal to the inverse of the an eigenvalue corresponding to a delocalized eigenvector [31]. The view of a delocalized eigenvector ${{\bf{v}}}_{m}$ inducing the epidemic threshold at ${\lambda }_{{\rm{c}}}=1/{{\rm{\Lambda }}}_{m}$ has been recently challenged in the case of SF networks, and in general in networks with a degree distribution decaying slower than exponentially; for such networks it has been shown that active states arise from the dynamic interplay and mutual reactivation of hubs, rather than from the activation of a precise delocalized eigenvector [33, 36].

How do such considerations apply to HMNs? Are active states promoted by a delocalized eigenvector? Or are they a result of a reactivation mechanism? Are reactivation mechanisms possible, provided that HMNs often exhibit exponential degree distributions [19], thus lacking highly connected hubs? We mentioned above that HMNs are now well-known to exhibit a Griffiths phase, whose upper bound is given by the actual critical point ${\lambda }_{{\rm{c}}}$ [19]. Such epidemic threshold does not comply with the QMF prediction, leading to ${\lambda }_{{\rm{c}}}\gt {\lambda }_{{\rm{c}}}^{\mathrm{QMF}}=1/{{\rm{\Lambda }}}_{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 ${\bf{A}}$, that is, it is not limited to the principal eigenvector ${{\bf{v}}}_{1}$. We show here the extent to which the basic assumptions of QMF are altered in HMNs, thus making the estimate of ${\lambda }_{{\rm{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 ${{\bf{v}}}_{1}$: as argued above, the sum in equation (1) would be dominated by the m = 1 term and, taking the equal sign (QMF approximation), ${{\boldsymbol{\rho }}}^{\infty }$ would be simply given by ${{\bf{v}}}_{1}$. Such a strong statement works remarkably well in networks with strong small world properties [34] or entangled connectivity patterns [37, 38], but fails to describe systems like ours in which spectral gaps are small [19, 39, 40]. We propose a way to relax the above assumption, by considering a candidate steady state ${{\boldsymbol{\sigma }}}^{\infty }$ as a linear combination of the first m^* eigenvectors, that is ${{\boldsymbol{\sigma }}}^{\infty }\approx {\sum }_{m=1}^{{m}^{* }}({{\bf{v}}}_{m}\cdot {{\boldsymbol{\sigma }}}^{\infty }){{\bf{v}}}_{m}$. In particular, we can always choose our eigenvector basis in such a way that ${{\bf{v}}}_{m}\cdot {{\boldsymbol{\sigma }}}^{\infty }\geqslant 0$ for all $m\leqslant {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 ${{\boldsymbol{\rho }}}^{\infty }={{\boldsymbol{\sigma }}}^{\infty }$? By inserting our expansion for the steady state in equation (1), and using the linear independence of eigenvectors, we obtain the condition

$$\begin{eqnarray}&&(\lambda {{\rm{\Lambda }}}_{m}-1)({{\bf{v}}}_{m}\cdot {{\boldsymbol{\sigma }}}^{\infty })\geqslant 0\ \ \ \forall m\leqslant {m}^{* }.\end{eqnarray} \tag{ 2 }$$

Equation (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 ${{\boldsymbol{\rho }}}^{\infty }\ne 0$:

• (a)  
  if the spreading rate is small, i.e. $\lambda \lt 1/{{\rm{\Lambda }}}_{1}$ (and thus, $\lambda \lt 1/{{\rm{\Lambda }}}_{m}$ for all m) then the candidate non-trivial steady state ${{\boldsymbol{\sigma }}}^{\infty }$ is unstable;
• (b)  
  even if $\lambda \geqslant 1/{{\rm{\Lambda }}}_{m}$ for some m, the localization of ${{\bf{v}}}_{m}$ induces ${{\bf{v}}}_{m}\cdot {{\boldsymbol{\sigma }}}^{\infty }\approx 0$ (positive but vanishing): the candidate steady state is weakly stable in principle, however dynamic fluctuations are bound to deactivate it over time. In cases like this, a virtually stable state results in fact metastable and relatively short-lived [31], in agreement with the Griffiths phase picture [19, 32].Thus, a nonzero steady state is still possible, if one of the following two conditions is met:

• (i)  
  There are some delocalized eigenstates, i.e. a pathological region in the spectrum exists around a given ${{\rm{\Lambda }}}_{m}$, characterized by high density of states and/or low localization, such that it can trigger a stable ${{\boldsymbol{\rho }}}^{\infty }$ with ${{\bf{v}}}_{m}\cdot {{\boldsymbol{\rho }}}^{\infty }\gg 0$ and ${\lambda }_{{\rm{c}}}\approx 1/{{\rm{\Lambda }}}_{m}$. In such a case, a spanning cluster of active nodes emerges by activation of a delocalized eigenvector. This mechanism is in essence the same proposed for SF networks in [31].
• (ii)  
  A finite number of localized states can dynamically coalesce sustaining a global active state (see figure 3), i.e. no pathological $({{\rm{\Lambda }}}_{m},{{\bf{v}}}_{w})$ pairs exist, however if λ is large enough, $\lambda {{\rm{\Lambda }}}_{m}-1\gt 0$ for many values of $m\lt {m}^{* }$: a large enough number of potentially short-lived states is generated, large enough to sustain (or coalesce into) a long-lived steady state ${{\boldsymbol{\rho }}}^{\infty }\gt 0$, however with no clear criterion to identify what a large enough m^* is. In such a case, active nodes may not form a spanning cluster. This phenomenology of dynamic coalescence is analogous to the reactivation picture proposed for SF networks [33], although it is obtained here using spectral methods, and highlighting the role of the higher spectral edge of ${\bf{A}}$. In SF networks hubs mutually reactivate each other. In HMNs, which may lack proper hubs, we conjecture that densely connected modules can generate dynamic pattern of reactivation.

Zoom In Zoom Out Reset image size

In order to clarify which of the two scenarios holds in the case of HMNs, and how the value of ${\lambda }_{{\rm{c}}}$ can be related to spectral properties of ${\bf{A}}$, we conducted an extensive numerical study of HMNs of different α, targeting both the SIS dynamics simulations and the network spectral properties. Details on how ${\lambda }_{{\rm{c}}}$ is measured are given in the next section.

Our numerical results provide a clear indication that the relevant scenario for HMNs is indeed (ii), as we found that the value of $1/{\lambda }_{{\rm{c}}}$ does not correspond to any special eigenvalue associated with any singular density of states or a delocalized eigenvector. This observation is exemplified in figure 4, where the higher spectral edge of HMNs is shown, including eigenvalues in the range $1\leqslant {\rm{\Lambda }}\leqslant {{\rm{\Lambda }}}_{1}$, and the corresponding $1/{\lambda }_{{\rm{c}}}$ is highlighted. In order to rule out completely the possibility of delocalized eigenvectors playing an important role, we choose to show the somewhat extreme case of HMNs with n_i = 10 and ${M}_{0}=4$, corresponding to a high value of $\alpha =2.5$: such HMNs display comparatively weak localization and even a moderate number of almost-delocalized eigenvectors. Eigenvalue densities and the corresponding IPRs are shown for varying Λ. While both plotted quantities exhibit a feature-rich behavior, with alternating peaks of density of states and localization, the simulated $1/{\lambda }_{{\rm{c}}}$ does not correspond to any of such features. More importantly, $1/{\lambda }_{{\rm{c}}}$ is not locked to any region in the spectrum: for different values of α (or n_i), $1/{\lambda }_{{\rm{c}}}$ is found to shift independently with respect to spectral peaks, in some cases corresponding to less localized states, in other cases to more localized ones (the latter scenario being the one depicted in figure 4). Even in such a case of a weakly localized HMN as the one in figure 4, the active state does not seem to emerge because of a more weakly-localized or almost-delocalized eigenvector. Of course, these conclusions are all the more valid in the case of more strongly localized HMNs (lower α, lower n_i, lower M₀), where the IPR reaches significantly higher values [19].

Zoom In Zoom Out Reset image size

These observations allow us to rule out the scenario of a delocalized eigenvalue triggering a true active state, and to conclude that active states in HMNs are the result of the dynamic coalescence mechanism discussed above. Simulation data confirm this picture, as exemplified graphically in figure 3, illustrating the lack of any spanning cluster of active nodes, even in the active state: 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 dynamic interplay.

Finally, we provide a qualitative explanation for the multi-peaked shape of the spectral tail, as shown in figure 4. In order to understand its origin, let us consider a network consisting of two densely connected modules, a and b, of the same size and comparable features. Let us start from the situation in which the two modules are isolated from one another: in this case each of them will exhibit a largest eigenvalue, ${{\rm{\Lambda }}}_{1}^{a}$ and ${{\rm{\Lambda }}}_{1}^{b}$ (of comparable values, ${{\rm{\Lambda }}}_{1}^{a}\approx {{\rm{\Lambda }}}_{1}^{b}$), and two large spectral gaps, since ${{\rm{\Lambda }}}_{1}^{a}\gg {{\rm{\Lambda }}}_{2}^{a}$ and ${{\rm{\Lambda }}}_{1}^{b}\gg {{\rm{\Lambda }}}_{2}^{b}$. The resulting network, comprising the two modules, will also have ${{\rm{\Lambda }}}_{1}^{a}$ and ${{\rm{\Lambda }}}_{1}^{b}$ as its two largest eigenvalues. By weakly connecting modules a and b, e.g. by drawing a small number of links, the resulting connected network will still display two largest eigenvalues ${\tilde{{\rm{\Lambda }}}}_{1}^{a}$ and ${\tilde{{\rm{\Lambda }}}}_{1}^{b}$, which are only perturbatively different from the original ${{\rm{\Lambda }}}_{1}^{a}$ and ${{\rm{\Lambda }}}_{1}^{b}$: one still has ${\tilde{{\rm{\Lambda }}}}_{1}^{a}\approx {\tilde{{\rm{\Lambda }}}}_{1}^{b}$, although now the Perron–Frobenius theorem ensures that they are distinct, while smaller eigenvalues remain significantly smaller. The resulting network has a small spectral gap, but the higher eigenvalues are well separated from the rest. Repeating this operation for all pairs of modules and all hierarchical levels as required to build a HMN, as in figure 1, it is easy to see how, every time two super modules are connected, pairs of eigenvalues accumulate at fixed locations in the spectrum, generating the peaks in figure 4. The spread in such peaks is given by the stochastic nature of the building process, namely by the wiring probabilities ${\pi }_{i}$. 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]. We also note that the extreme separation between peaks in figure 4 is a result of the choice of a high n_i = 10. Smaller values of n_i lead to closer and closer peaks, until no peak separation is visible any more [19].

Given the complexity of our dynamic scenario, in which a large but undefined number of eigenvalue-eigenvector pairs is expected to play a role, is it still possible to relate ${\lambda }_{{\rm{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 ${\bf{A}}$ increases accordingly. A simple argument to shed light on this observation is as follows: according to the spectral picture described above 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 ${{\bf{A}}}^{(i)}\sim {c}_{i}\left(\begin{array}{cc}0 & \alpha \\ \alpha & 0\end{array}\right),$ which will contribute a larger eigenvalue proportional to ${c}_{i}\alpha$ to the higher spectral edge of ${\bf{A}}$, where to a first approximation c_i is a factor that depends on the hierarchical level i (roughly speaking, the factor determining the location of the peaks in figure 4). We can corroborate our view by looking at figure 5(a), where it is shown that the value of the generic mth eigenvalue is strictly proportional to α, for a wide range of values of m 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 ${{\rm{\Lambda }}}_{1}$ only, effectively changing the size of the spectral gap [34].

Zoom In Zoom Out Reset image size

The exceptional behavior of HMN spectra, for which a single tuning parameter is able to tune a wide portion of the spectral edge, turns crucial for our study of the epidemic threshold: provided that ${\lambda }_{{\rm{c}}}$ scales as the inverse of some relevant eigenvalues in the higher spectral edge, responsible for the dynamic coalescence mechanism described above, the similarity property of eigenvalue spectra ${{\rm{\Lambda }}}_{m}\sim \alpha$ provides us with a simple prediction for the epidemic threshold, which can be expressed as:

$$\begin{eqnarray}&&{\lambda }_{{\rm{c}}}\sim \displaystyle \frac{1}{\alpha }.\end{eqnarray} \tag{ 3 }$$

In order to verify this conjecture, we measured ${\lambda }_{{\rm{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 ${\lambda }_{{\rm{c}}}$ make their estimate a more delicate task [41]). A detailed description of our procedure is given in figure 6. We performed simulations 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⁵–10⁶ nodes) the dynamics display no appreciable size effects. In particular, we find that steady states for $\lambda \gt {\lambda }_{{\rm{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\to \infty$ limit.

Zoom In Zoom Out Reset image size

The validity of equation (3) is then confirmed by our simulation results for ${\lambda }_{{\rm{c}}}$ in figure 5(b) (main figure and inset): a remarkable scaling property

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{m}\sim \displaystyle \frac{1}{{\lambda }_{{\rm{c}}}}\end{eqnarray} \tag{ 4 }$$

holds for all eigenvalues in the spectral edge and, by virtue of the ${{\rm{\Lambda }}}_{m}\sim \alpha$ relationship, the conjecture in equation (3) is finally verified (figure 5(b), inset). In passing, we note that as all eigenvalues in the range scale with $1/{\lambda }_{{\rm{c}}}$, this must happen, in particular, for ${{\rm{\Lambda }}}_{1}$: while the principal eigenvalue alone cannot justify the value of the epidemic threshold, the standard ${\lambda }_{{\rm{c}}}^{\mathrm{QMF}}=1/{{\rm{\Lambda }}}_{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\propto \alpha$, we can rewrite our estimate as

$$\begin{eqnarray}&&{\lambda }_{{\rm{c}}}\sim \displaystyle \frac{1}{D}.\end{eqnarray} \tag{ 5 }$$

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 ${\lambda }_{{\rm{c}}}$. The picture emerging from our results is that of a complex dynamic scenario, as the one described by equation (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\to \infty$ limit, equation (5) predicts a vanishing epidemic threshold, recovering the well-known QMF result for SF networks, which are intrinsically infinite dimensional. The actual study of the existence of a finite ${\lambda }_{{\rm{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 α. Figure 7 shows the steady state density of active nodes, ${\rho }^{\infty }$, as a function of $(\lambda -{\lambda }_{{\rm{c}}})/{\lambda }_{{\rm{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 ${\rho }^{\infty }\sim {(\lambda -{\lambda }_{{\rm{c}}})}^{\beta }$, 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.

Zoom In Zoom Out Reset image size

We discussed the essential reasons why the QMF approximation does not provide a correct prediction for the epidemic threshold in HMNs. 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, corresponding to active modules. However such 'reservoirs' of instability can create a global active state owing to fluctuations that transiently activate other modes, effectively connecting diverse localized clusters. Such physical picture is analogous to the hub-mediated reactivation mechanism, recently proposed for SF networks, where degree distributions are power-laws [33, 36]. In HMNs, where degree distributions are exponential and high-degree hubs are absent, reactivation occurs at the inter-module level and across hierarchical levels.

We propose a spectral framework, which extends the QMF approach by taking into account the contribution of a whole range of eigenvalues in the higher spectral edge of ${\bf{A}}$, and their corresponding eigenvectors. We find that even though for HMNs the epidemic threshold depends on a (large) number of such eigenvalues, it turns out to be inversely proportional to a unique parameter: the network topological dimension, D. While our numerical results corroborate this finding 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 ${\lambda }_{{\rm{c}}}$ on D ensures that, being brain HMNs finite dimensional, the activity propagation 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, the optimal tuning of a single parameter, e.g. during the development or pruning of neuronal networks, allows the system to settle in the correct dynamic regime.

In a broader context, the relationship between spectral properties of HMNs and their topological dimension D may prove of great interest to the compelling problem of random walks and first-passage phenomena in HMNs, whose study can be addressed harnessing the knowledge of graph spectra [15, 42, 43]. In passing, we note that the absence of a spanning cluster of active nodes in HMNs, as found here for SIS dynamics, is reminiscent of the remarkable breakdown of ergodicity of random walks in Dyson hierarchical graphs [43].

An intriguing question is: how broad in general the eigenvalue range ${ \mathcal R }$ characterized by the scaling property in figure (5) is, and more importantly how ${ \mathcal R }$ it changes with system size. Our current data allows us to conclude that in HMNs of size $N={2}^{17}$, at least the 10³ largest eigenvalues exhibit the scaling behaviour, although the breakdown of that behavior for smaller eigenvalues does not appear as a sharp transition, and a clear change in regime is hard to locate. Considering the size scaling properties shown in figure 4, where the number of eigenvalues in each bin in the spectral tail increases linearly with N, we may only conjecture at this point that also the number of eigenvalues in ${ \mathcal R }$ increases linearly with N. A validation of this conjecture is beyond the scope of this manuscript, and is left for future work.

Finally, the possibility that in certain networks the delocalization of individual eigenvectors may control the emergence of active states, as proposed in [31] and initially conjectured for HMNs in [19], remains an appealing scenario, in which structure and function are tightly interconnected. What requirements a network structure should meet in order to unequivocally yield such behavior is still subject of investigation.

In conclusion, we have studied analytically and 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 D of a HMN. We have shown how α and, equivalently—D are responsible for the tuning of the epidemic threshold ${\lambda }_{{\rm{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.