Anomalous Lifshitz dimensions in hierarchical networks of brain connectivity

Network models of neural connectivity and function often invoke the ability of the brain to localize activity in distinct modules simultaneously. The propensity of a network to do the opposite instead, that is to transmit and diffuse information homogeneously, is measured by its spectral dimension, a quantity that is easily accessible through analyses of random walks, or equivalently diffusion processes. Here we show that diffusive dynamics in hierarchical modular network models, representing brain connectivity patterns, exhibit a strongly anomalous features, pointing to a global asymptotic slowdown at large times and to the emergence of localization phenomena. Using theoretical modeling and very-large-scale computer simulations, we demonstrate that the spectral dimension is not defined in such systems and that the observed anomalous dynamical features stem from the existence of Lifshitz tails in the lower spectral edge of the Laplacian matrix. We are able to derive the correct scaling laws relating the spectral density of states and anomalous dynamics, emphasizing the fundamental role played by the Lifshitz dimension. Our work contributes to establishing a theoretical framework for anomalous dynamical features, such as activity localization and frustrated synchronization in hierarchical and hierarchical-modular networks and helps contextualize previous observations of sub-diffusive behavior and rare-region effects in brain networks. More in general, our results, help shedding light on the relation between structure and function in biological information-processing complex networks.

Understanding the interplay between dynamical processes and the architecture of the networks embedding them is a fundamental problem in diverse fields including material science, genetic regulation and neuroscience.Dynamical features and patterns of activity are often affected or controlled by key structural features of the underlying network, such as the degree distribution, degree correlations, modular organization, k-core structure, etc. [1][2][3][4][5].However, given that such features are usually not independent, a more systematic way to tackle the problem of the interplay between structure and dynamics relies on the use of spectral-graph characterizations of the network architecture [6,7] and, importantly, the network dimension.Statistical mechanics teaches us that dynamical processes such as diffusion, vibrational excitations, and critical properties near second order phase transitions exhibit universal behavior, which depends crucially on the lattice (Euclidean) dimension [1-3, 8, 9].The case of heterogeneous networks is more complex, since multiple and diverse generalizations of the concept of dimension have been proposed [10][11][12].Nevertheless, compelling pieces of evidence show that dimensionality measures are effective determinants of dynamics and activity in networked complex systems.The simplest example is provided by networks with the small-world property [13][14][15], which exhibit diameters that grow only logarithmically with the network size N and, consequently, with diverging Hausdorff dimension.We recall that the Hausdorff dimension d H , also called "topological dimension" in the literature [16,17], can be computed easily starting from the number u i (r) of nodes within distance r from node i: if u i (r) ∼ r d (where .stands for the average over all nodes in the network), implying d H = d Diverging d H in these systems implies enhanced transmission, signal propagation and synchronizability.A somewhat more complex example of how dimensionality controls activity patterns in networks is provided by hierarchical-modular networks, as models e.g. for brain connectivity [18,19].It was first pointed out that the hierarchical-modular organization of brain regions results in network models of finite Hausdorff dimension d H and, at odds with small-world graph topologies, intrinsically large diameters [20].The large-world property resulting from finite d H in hierarchical-modular networks, a purely structural feature of the network, has been linked to signatures of anomalous activity patterns in brain network models, including among others: sustained activity [21], sub-diffusive dynamics [22], localization phenomena and stretched criticality in the form of Griffiths phases [17,23,24], broad avalanche distributions [17,25], states of localized and "frustrated" synchronization [26][27][28][29][30], rounding of first-order phase transitions [31], and ergodicity breakdown [32,33].Importantly, some of these anomalous dynamical traits are, in fact, considered essential to the ability of brain networks and of brain-inspired hierarchical architectures to achieve an optimal balance between segregation and integration [18], allowing them to conduct multiple tasks simultaneously [34] and entailing optimal computational capabilities.Let us also note that a significant part of the above-mentioned phenomenology uses concepts, such as Griffiths phases, first introduced to study localization phenomena in quantum systems described, e.g., by a random tight binding Hamiltonian [35], and later extended, for instance, to the Laplacian matrix of a graph [36,37].
In this paper, we aim at providing a theoretical foundation for the phenomenological observations of anomalous behavior and localization effects obtained so far, establishing for the first time a clear link between the emergence of localized dynamical patterns of activity and synchrony in hierarchical-modular networks and spectral properties of their underlying graphs.The fundamental concept, allowing us to develop our approach is the spectral dimension d s of a graph -as well as an important extension of it-which can be defined and measured by simple random walk (RW) analyses [38][39][40].Given the probability P ij (t) that a random walker starting at time t 0 = 0 from node i arrives at node j after t steps, one can compute the average return probability as , d s is defined as the (average) spectral dimension of the network [40].While for infinitely large networks further complications arise due to the possibility of transient random walks [39], here we focus on networks of finite, albeit very large size N , so that the definition above is intended to hold asymptotically, in the limit of large N and large t.The spectral dimension, not unlike the Hausdorff dimension introduced above, is a generalization of the concept of dimension in a Euclidean continuum.Actually, in discrete lattices d s = d H , so that both generalizations agree with the Euclidean dimension in a continuum system.This equivalence does not hold in general in heterogeneous networks, nor does it in deterministic fractals [8].We notice in particular that while d H is a purely structural measure, d s is an observable of a diffusion process operating on the network, and as such it provides us with a probing tool for dynamical signatures of localization and slowing down, and a first approximation in cases, like that of brain activity, in which more complex dynamics are at play.It was initially conjectured that Griffiths phases in heterogeneous networks can only occur in the finite d H case [16,17].This view was recently challenged by Millán et al., who found that even networks with infinite d H can exhibit similar dynamical regimes, provided that d s is finite instead [41].The relationship between different definitions of a network dimension and their use to predict the emergence of anomalous dynamical patterns thus remains an open question to be fully clarified.
With these considerations in mind, we analyse the spectral dimension in hierarchical-modular network models of brain connectivity [17], with the objective of quantifying how the basic traits of brain activity localization, which are captured by these simple network models, are reflected by d s and stochastic diffusion null models.To this end, we conducted very-large-scale RW simulations in hierarchical-modular network models of tunable Hausdorff dimension d H and computed the average return probabilities R(t).As for the synthetic networks considered for computational analyses, we chose to work with the model proposed in [17,42], which comes with a single effective parameter α (the connectivity strength) and an additional parameter s fixing the number of hierarchical levels.In this model the average network degree k and (asymptotically) the Hausdorff dimension d H are both proportional to α so that sparser networks have smaller Hausdorff dimension.While other models proposed in the literature might differ in the choice of parameters [21,25,43,44], we believe that the conclusions of the current work remain unchanged.In order to capture the large network size limit of the system, we performed random-walk computer simulations on networks of sizes up to 2 25 ≈ 3 × 10 7 , and for time windows large enough as to ensure that all walkers return to the starting node (which, we recall, is possible because N is finite in our case).Fig. 1 shows the return probabilities for a typical choice of parameters (N = 2 25 , s = 24), and for increasing values of α, and thus increasing Hausdorff dimension d H .One can immediately see the anomaly in the asymptotic behavior of R(t): while for intermediate nt the curves develop a heavy tail, resembling a power law, the slope of such tails apparently decreases in absolute value upon increasing aα, and later develops a non trivial large-t bump, significantly before the finite size cutoff appears.In cases in which the spectral dimension is defined, we would expect the slope to increase in absolute value with α, implying that d s increases with d H and one would expect that behavior to be the asymptotic (t → ∞) one.In the present study, instead, one is forced to conclude that dynamical slowing down is so radical that the asymptotics are given by the excess returns at very large t (the bump in the curves above) and the average spectral dimension is, as a consequence, undefined.The absence of a well-defined spectral dimension is an already interesting result within the study localization and dynamical slowing down in models of brain connectivity.Suppression of diffusion and free flow are considered signatures of anomalously slow dynamical regimes, which ensure the balance between global integration and functional modularity of brain activity [20,22,45].To proceed, let us first elucidate the nature of the asymptotic behavior of the return probabilities.Fig. 2 reveals that the large t dependence of R(t) is dominated by a stretched exponential behavior R(t) ∼ e −t β , governed by a non trivial positive β < 1 (for ease of notation, we measure t in dimensionless units).We call β the anomalous exponent, as its value quantifies the dynamical slowing down with respect to the standard scenario where a spectral dimension is defined.To confirm this view, Fig. 2 also shows the dependence of β on α and thus on the Hausdorff dimension of the network.Sparser networks exhibit a more strongly anomalous behavior, with the anomaly exponent loosely proportional to the Hausdorff dimension.
In order to make further progress, we exploit methods of spectral graph theory [6,7].Let us write the exact master equation, describing a random walk as a timecontinuous Markov process on a generic undirected and unweighted network, encoded in an adjacency matrix A as follows: where q(t) is the column vector, whose generic element q i (t) represents the probability of the random walker to reach node i at time t, and L RW is the random walk Laplacian matrix with elements L RW ij = δ ij − A ij /k j with i, j ∈ 1, 2, ...N .Here A ij is the generic element of A, equal to 1 if nodes i and j are linked and 0 otherwise, and k j = N i=1 A ij is the degree of node j.In order to compute the solution to Eq.1, we introduce the normalized Laplacian L, defined by the similarity transformation , where D is the (diagonal) degree matrix of generic element D ij = δ ij k j .L is symmetric and diagonalizable, and by virtue of their similarity, L and L RW have the same spectrum of eigenvalues, although different eigenvectors.We choose to express the solutions of Eq.1 through the eigen-decomposition of L, q j (t 0 ), in which we introduced the heat kernel K(t) of generic element [46]: where λ m is the m-th eigenvalue of L, V im is the i-th component of the eigenvector of L associated with λ m , and t 0 = 0 without loss of generality.This well-known identity allows us to connect the spectral perspective with random-walk simulation results: one can easily see that the average return probability R(t) is related to the trace of K(t) (or heat trace) through the simple relationship The eigenvalues λ m are all real and non-negative and, under the assumptions that the network is undirected and connected, the 0 eigenvalue is unique and one can always choose the labeling 0 = λ 1 < λ 2 ≤ λ 3 ≤ • • • [6].As a consequence, a random walk always reaches a steady state above a time scale given by the interplay of the smallest nonzero eigenvalues [6].By making a continuum spectrum approximation for λ ≥ λ 2 , one can introduce the density of states (eigenvalue density distribution) ρ(λ), so that Eq.3 can be approximated by its continuum limit where the integral is dominated by the contribution of the lower spectral edge.Under these assumptions, the density of states ρ(λ) and the return probability are related through a simple Laplace transform operation.
In lattices, which are endowed with an integer spectral dimension, this result is well known and leads to ρ(λ) ∼ λ (ds/2)−1 , which shows the relationship between the Laplacian spectra and the spectral dimension [40].This result is also well known in terms of vibrational frequencies ω ∝ √ λ in lattices and deterministic fractals, leading to a density of states ρ(ω) ∼ ω ds−1 [8].
While the approximation in Eq.4 holds in the general case, we expect the conclusions regarding d s to be radically different in case of hierarchical-modular networks, for which our numerical results reveal anomalous stretched-exponential tails of the return-probability function R(t).This anomaly is indeed reflected in the lower spectral edges and in particular in ρ(λ) as we show below.It was hypothesized in the past [17] that low eigenvalues of L in such networks form a continuous spectral tail, a Lifshitz tail, in analogy with the random Hamiltonian operators in a tight-binding Schrödinger equation [37].In particular, it was noted that Lifshitz tails may be relevant to assess the subcritical dynamics in models of epidemic spreading, in the framework of a linearized quenched mean-field approximation [23].The randomwalk problem that we study here, instead, can always be mapped exactly to the quantum problem, with the energy eigenvalues E m of the quantum problem being replaced by the Laplacian eigenvalues λ m [37].Under the hypothesis of Lifshitz tails, the integrated density of states of L (the cumulative eigenvalue density distribution) is expected to exhibit a tail of the general form [35] where λ 0 is the lower bound of the continuum spectrum, which we can set equal to 0 in the case of hierarchicalmodular networks, as they possess vanishing spectral gaps (0 < λ 2 1) [17], and c 1 and c 2 are constants.The real number d L coincides with the space Euclidean dimension in the original Lifshitz argument for continuum quantum problems.In the present discrete classical case, since that we cannot yet provide a Lifshitz-like argument, we simply name d L the Lifshitz dimension of the problem.In the light of the above considerations, the density of states ρ(λ) = dN /dλ is dominated by the following low-λ tail Observe that Eq.6 differs significantly from the lattice prediction above ρ(λ) ∼ λ (ds/2)−1 .Is such a difference a spectral signature of the anomalous behaviour (i.e. the stretched-exponential tail of the return probabilities R(t) and the lack of a well-defined spectral dimension d s ) encountered in RW simulations?As the two representations connect through Eq.4, one can compute R(t) as in Eq.4, using the hypothesis of a Lifshitz tail from Eq.6.
While the integral involved in the calculation, of the form w(t) = e g(λ,t) dλ with g(λ, t) = c 2 λ −dL/2 − λt, is highly non-trivial in the case of real d L , here we are only interested in its asymptotic t → ∞ behavior, which is captured by the values of λ for which g(λ, t) is maximum.This is readily obtained through the saddle-point approximation w(t) ≈ e g(λ * ,t) , with λ * the location of such maximum, leading to the final result Eq.7 remarkably recovers a stretched exponential tail behavior, as reported above for computer simulations, confirming for the first time that the dynamical slowing and Using the rescaling from Eq.9, Lifshitz tails appear as straight lines.By choosing values of dL obtained from the RW simulation results through Eq.8, Lifshitz tails collapse in a single curve, confirming that the anomalous dynamical behavior has its origin in the spectral properties of L. lack of spectral dimension can be attributed to the existence of Lifshitz tails and providing us with an interpretation of the anomalous exponent in terms of the Lifshitz dimension: We can conclude that in the present hierarchical-modular network model, not only Lifshitz tails explain the anomalous dynamics and the lack of a well-defined spectral dimension, but also d L provides us with a meaningful dimensionality measure, generalizing the behavior of quantum systems in the continuum, where the Lifshitz dimension identifies the spatial dimension.So far, we have only hypothesized that the lower spectral edge of L exhibits a Lifshitz tail of the form given by Eqs.5 and 6.Now we corroborate such a hypothesis by verifying, not only that the integrated density of states N (λ) obeys the tail behavior in Eq.5, but also that the exponent d L /2 governing it generates the anomalous exponent β of the dynamical simulation, as predicted by the result in Eq.8.To this end, we notice that according to the prediction above, one expects for large 1/λ.We obtain d L = 2β/(1 − β) from Eq.8, and using the values of β obtained from the initial random walk simulations one can easily verify Eq.9 by computing the lower spectral edges of hierarchical-modular networks of the same type.Computational results, shown in Fig. 3 clearly confirm the linear dependence predicted by Eq.9 for small values of lambda.In other words, the prediction based on the Lifshitz tails assumption is correct: the spectra of hierarchical-modular networks exhibit Lifshitz tails, with an associated Lifshitz dimension d L , and their concomitant anomalous dynamical behavior is controlled by d L .This behavior can be called localization, in analogy with the case of quantum systems, where localization stands for absence of diffusion.This behavior is often quantified by looking at the corresponding eigenvectors, whose components almost vanish everywhere except in specific locations in the system.This property, known as eigenvector localization [35] has been long observed in models of epidemic spreading on network classes of techno-social relevance [23,47,48], as well in problems inspired by brain connectivity [17,44] and biological materials [49].
It is noteworthy that the emergence of classical Lifshitz tails in network spectra has been rigorously proved in Erdős Rényi graphs below the percolation threshold [37], i.e. for networks that have not yet developed a giant connected component.Our results suggest a remarkable property of hierarchical-modular networks, which exhibit Lifshitz tails while being connected, i.e. while possessing a single connected component.We believe that Lifshitz tails, and the resulting anomalous dynamical behavior, may be observed in general in network models exhibiting similar localization properties, such as hierarchical trees displaying patchy percolation [50], and dense hierarchical Dyson networks [34], where ergodicity is known to break down in the thermodynamic limit [32,33].In fact, we propose the emergence of Lifshitz tail and their associated Lifshitz dimensions as a criterion for the existence of strong localization in networks.
Our focus on the spectral dimension and its undefined nature in hierarchical networks allows us to establish a connection between network structural properties and anomalous dynamics in systems such as brain networks, where the ongoing structure vs. function debate has long dealt with the issues of relating activity patters to specific anatomical arrangements, or alternatively presenting them as emergent, or self-organized.Our results clearly show a connection between dynamical slowingdown, localization properties, and Laplacian spectra; let us emphasize that such a correspondence is clean-cut because of the simplicity of the random walk model, and possibly because of the simplifying assumptions in the choice of our network model.While diffusion -lacking any form of non-linearity-is arguably a very crude simplification of neural dynamics on the structurally complex human connectome [19], it has been found that the eigenvectors of a diffusion problem on the connectome are relevant in predicting functional patterns of neural activity [51,52].More in general, the Laplacian matrix provides the linearization of oscillator models near a synchronization transition [53] and the relevance of its spectrum in the problem of brain synchronization has been discussed in the literature and related to the observation of frustrated synchronization [26][27][28]41].While some of those results were based on the hypothesis of the existence of Lifshitz tails, here we are able to prove that hypothesis and rationalize those results within a proper theoretical framework, where the slowing down of synchronization processes is governed by the Lifshitz dimension d L , which effectively tunes the dynamical anomalies.Let us also note that, while the approach here resorts to unweighted networks, the Laplacian formalism lends itself to the introduction of weights [32,33], a quantity that in neuroimaging encodes the number of connections between pairs of brain regions.Beyond the case of pairwise interactions, recent advances in integrating the concepts of diffusion and spectral dimension within the broader field of algebraic topology and simplicial complexes [54][55][56] provide a promising avenue to strengthen the theoretical framework for the localization phenomena that we discuss here to describe, for instance, systems with higher-order interactions between their components/nodes.
In conclusion, we have established a theoretical framework for the prediction of anomalous dynamics in hierarchical network models of interest in brain modelling.To our knowledge, hierarchical-modular networks constitute the first heterogeneous network model displaying Lifshitz tails above the percolation threshold, and the first not to exhibit a power-law behavior for the average return probability and a well-defined spectral dimension.Being able to connect these two singular features allows us to rationalize previous experimental observations of activity localization in the brain and their numerical models, where spectral anomalies and Lifshitz tails where only hypothesized.We believe that these results will stimulate interest and further work, in e.g.computational neuroscience, as a way to advance the knowledge on how the brain achieves an optimal balance between segregation and integration.At the same time, we are confident that the present framework will provide us with more powerful tools for the tunability and controllability of network models exhibiting strong localization, relevant in the design of synthetic networks for brain-inspired neuromorphic computing.

FIG. 1 .
FIG.1.Return probabilities for hierarchical-modular networks with N = 225 , s = 23 and increasing α.Even at large sizes, no clear power law decay is visible and the asymptotic behavior is dominated by a stretched exponential tail, before the finite size cut-off takes over.Similar results have been found for smaller sizes and different choices of s.

24 FIG. 2 .
FIG.2.Quantitative analysis of the stretched exponential behavior found in Fig.1.Left: semi-logarithmic-scale plot of the return probability.For each curve (each value of α) the value of β is chosen, which makes the stretched exponential tail appear as a straight line.Color scheme is as in Fig.1.Right: the values of β from the left panel (N = 225 , s = 23, blue line) and the same values for a smaller system (N = 224 , s = 22, gray line).

FIG. 3 .
FIG.3.Lower spectral edge of L. Using the rescaling from Eq.9, Lifshitz tails appear as straight lines.By choosing values of dL obtained from the RW simulation results through Eq.8, Lifshitz tails collapse in a single curve, confirming that the anomalous dynamical behavior has its origin in the spectral properties of L.