Persistence of hierarchical network organization and emergent topologies in models of functional connectivity
Ali Safari, Paolo Moretti, Ibai Diez, Jesus M. Cortes, Miguel Ángel Muñoz
PPersistence of hierarchical network organization andemergent topologies in models of functional connectivity
Ali Safari a , Paolo Moretti a, ∗ , Ibai Diez b,c , Jesus M. Cortes d,e,f , Miguel ´AngelMu˜noz g,h a Institute of Materials Simulation, Friedrich-Alexander-Universit¨at Erlangen-N¨urnberg,Dr.-Mack-Str. 77, D-90762 F¨rth, Germany b Functional Neurology Research Group, Department of Neurology, Massachusetts GeneralHospital, Harvard Medical School, Boston, MA 02115, USA c Gordon Center, Department of Nuclear Medicine, Massachusetts General Hospital,Harvard Medical School, Boston, MA 02115, USA d Biocruces-Bizkaia Health Research Institute, Barakaldo, Spain e Departament of Cell Biology and Histology. University of the Basque Country(UPV/EHU), Leioa, Spain f IKERBASQUE: The Basque Foundation for Science, Bilbao, Spain g Departamento de Electromagnetismo y F´ısica de la Materia, Universidad de Granada,Granada E-18071, Spain h Instituto Carlos I de F´ısica Te´orica y Computacional, Universidad de Granada, GranadaE-18071, Spain
Abstract
Functional networks provide a topological description of activity patterns in thebrain, as they stem from the propagation of neural activity on the underlyinganatomical or structural network of synaptic connections. This latter is wellknown to be organized in hierarchical and modular way. While it is assumedthat structural networks shape their functional counterparts, it is also hypothe-sized that alterations of brain dynamics come with transformations of functionalconnectivity. In this computational study, we introduce a novel methodologyto monitor the persistence and breakdown of hierarchical order in functionalnetworks, generated from computational models of activity spreading on bothsynthetic and real structural connectomes. We show that hierarchical connec-tivity appears in functional networks in a persistent way if the dynamics is setto be in the quasi-critical regime associated with optimal processing capabil- ∗ Corresponding author
Email address: [email protected] (Paolo Moretti)
Preprint submitted to Neurocomputing February 11, 2021 a r X i v : . [ c ond - m a t . d i s - nn ] F e b ties and normal brain function, while it breaks down in other (supercritical)dynamical regimes, often associated with pathological conditions. Our resultsoffer important clues for the study of optimal neurocomputing architecturesand processes, which are capable of controlling patterns of activity and infor-mation flow. We conclude that functional connectivity patterns achieve optimalbalance between local specialized processing (i.e. segregation) and global inte-gration by inheriting the hierarchical organization of the underlying structuralarchitecture. Keywords: brain networks, functional connectivity, hierarchical modular networks,segregation integration
1. Introduction
The recent convergence of neuroscience and network science opens up snewopportunities to approach the study of brain function [1, 2, 3]. A fundamentalissue in this context is how structure and function are related [4, 5, 6, 7, 8] . Inthe context of network neuroscience, structure refers to the network mappingsof the brain, also known as connectomes, as derived from the actual anatomicalconnections between brain regions, also called “connectomes” [9, 10]. Structuralconnectivity (SC) is thus encoded in networks where nodes are coarse-grainedrepresentations of specific brain regions, and the links express the presence ofwhite-matter based connections between pairs of nodes, while weights associatedwith links conventionally measure the number of such connections. Structuralnetworks are then represented by weighted, non-negative and sparse adjacencymatrices, effectively containing the anatomical routing information of a brain.The direction of each connection can currently be recorded only for certaintypes of connectomes (e.g. mice or primates); in most cases, including thatof the human connectome, adjacency matrices are symmetric as directions ofconnections are still not detected.Functional connectivity (FC), instead, is often measured from neural activity2orrelations rather than actual anatomical pathways. As activity correlationscan be measured between any pair of nodes, the matrix representations forfunctional networks differ substantially from those for structural networks inthat they lose the property of being sparse. FC matrices are dense and cor-responding sparse-network representations can only be obtained by applyingarbitrary thresholding procedure. As a further source of complication, FC ma-trices lose the non-negativity of SC matrices too, as activity correlations aresigned, revealing the existence of both correlations and anti-correlations [11].This aspect makes the application of thresholds an even more delicate issue,subject to a large deal of arbitrariness, as in choosing for instance to excludeanti-correlations, or to study correlations and anti-correlations separately.We note that in spite of this complexity, FC data have been recorded foryears now and have led to ground-breaking advances in the the understandingof brain function. At the clinical level, FC data have been successfully relatedto the occurrence of brain pathologies [12, 13, 14, 15]. At the graph theoret-ical level, FC networks have facilitated the introduction of statistical physicsconcepts such as scale invariance, avalanches, criticality and localization in thebrain [16, 17, 18, 19, 20, 21]. Advances in acquisition and analysis of FC datahave allowed the introduction of novel and challenging concepts, such as thatof dynamical FC, i.e., the more ambitious study of the time dependence anddynamics of functional networks, obtained through the recording of multiplefunctional networks, each one over a short time window [22, 23, 24]. In all theseapproaches, the question of what relevant topological information is lost whenthresholds are applied remains open.Structural brain networks [25, 26] are well known examples of biologicalsystems, which exhibit a hierarchical modular structure. Hierarchical modu-lar organization is often described by simple mathematical models of synthetichierarchical modular networks (HMNs). Is it possible to say the same aboutfunctional networks? In principle, the answer to this question is affirmative asalso FC networks are known to exhibit modular and hierarchical patterns ofactivity [27, 10, 23]. But, how general can the answer to the previous question3e if functional network topology can vary so wildly because of arbitrary thresh-old choices? How persistent is the hierarchical organization of FC data, uponapplying different threshold and, more importantly, in different states of neuralactivity? Are structural differences between SC and their corresponding FCnetworks representative of pathological states? Answering these questions is anecessary step towards the understanding of the relationship between structuraland functional connectivity in the brain. In particular, the study of this relation-ship often relies on assessing the similarity between structural and functionalconnectivity. Similarity measures were proposed in the past (see for instance[8]), highlighting how “similarity” itself may be non-trivially dependent on thechoice of thresholds used to extract functional networks. In this work, we followa more fundamental approach, in the attempt to formulate a criterion to de-cide if a functional network is hierarchical and to do so irrespective of thresholdchoices.The interest in the hierarchical organization of brain activity patterns isrooted in the observation that hierarchical (and hierarchical-modular) networksexhibit desirable properties of robustness. Remarkably, hierarchical networks donot exhibit a single percolation threshold, where a giant component emerges to-gether with a scale-free distribution of cluster sizes [28]. Instead, in hierarchicalmodular networks, upon varying the control parameter p (i.e. upon removing afraction 1 − p of links or nodes) power-law distributions of connected componentsizes are encountered for a broad interval of values of p [29, 17]. The functionalcounterparts to this simple structural property are striking: • activity in hierarchical models of the brain is sustained even without anyfine-tuning of regulatory mechanisms [30]; • avalanches of neural activity are power-law distributed in size without theneed to fine tune spreading control parameters [17, 18], rationalizing theexperimental observation of scale invariant activity patterns [31]; • simple dynamical models, normally displaying clear-cut tipping points(phase transitions), exhibit instead extended quasi-critical phase (e.g. Grif-4ths phases [32, 33, 34, 18, 35, 19, 20]), where activity propagates by rareregion effects and states of local coherence emerge as chimera-like states[36, 37, 38].These observations corroborate the view that the hierarchical organizationof the structural network of anatomical connections is responsible for the brainability to localize activity, avoiding the opposing tendencies where active statesdie out (as encountered in advance stages of some neurodegenerative diseases)or invade the system (as typically occurs during pathological epileptic seizures).From the perspective of neurocomputing and information processing, HMNsthus provide an optimal architecture, which ensures the balance of activity segre-gation and integration [39], and results in enhanced computational capabilities,large network stability, maximal variety of memory repertoires and maximaldynamic range [18].In order to analyze the topology of the resulting FC networks, we exploittheir associated spectral properties, allowing us to properly identify hierarchicalfunctional networks . In particular, we focus on the well-known property of van-ishing spectral gaps —which, following the ideas in [18]— we use as a measureof the hierarchical organization. Networks characterized by dense connectivityand high synchronizability exhibit a large separation between the two largesteigenvalues of the adjacency matrix [40], or the two lowest distinct eigenvalues ofthe Laplacian matrix, a property that we refer to as large spectral gaps. Hierar-chical networks of interest in brain modeling, instead, display localized patternsof activity and, as a spectral counterpart, vanishing (although non-zero) spec-tral gaps [18, 36]. While the spectral characterization of HMNs includes manycomplex aspects, the single fact that they exhibit vanishing spectral gaps isquite remarkable and has been directly related with some important dynamicalfeatures such as self-sustained activity and local coherence [18, 36]. These as-pects make hierarchical networks unique and, more importantly, constitute thereason why hierarchical organization is essential in brain networks.In what follows, we generate FC data from Monte Carlo simulations of simple5odels of activity spreading in synthetic HMN models for structural networksand on real SC data. We show that, independently of the choice of thresh-olds, functional networks generated in the quasi-critical regime always displayhierarchical organization (as quantified through spectral gaps). This feature islost as soon as FC data is produced in the super-critical regime – i.e. above quasi-critical regime. Our results offer important insights regarding the topo-logical changes that FC data undergoes in pathological states (the super-criticalregime), providing clues for FC analysis as a diagnostic tool. At the same time,our methodology provides a new tool for the analysis of persistent features infunctional data, whose appearance is not an artefact of an arbitrary thresholdchoice.
2. Materials and methods
Structural networks.
We use computer generated models of HMNs, in orderto mimic realistic structural networks. As different types of algorithms areused in the literature to generate HMNs [30, 41, 18, 20, 21], in the followingwe refer the model proposed in [21] without loss of generality. We call α the connectivity strength of a HMN, as it is the main parameter controlling theemerging topology. The network is organized in densely connected modules ofsize M , which represent the level 0 of a hierarchy of links. At each hierarchicallevel i >
0, super-modules of size M i = 2 i M are formed, each one joiningtwo sub-modules of size M i − by wiring their respective nodes with probability π i = α/ i : the average number of links between two modules at level i is n i = π i M i − = α ( M / , i.e. proportional to α , regardless of the value of i .For a generic dynamic process running on a HMN in the N → ∞ limit, the effectof lowest-level modules becomes negligible (relegated to transient time scales)and the time asymptotics are dominated by the hierarchical organization: inthis regime α becomes the only relevant construction parameter. We remarkthat, being n i the number of links between any two modules of size M i − , thehighest possible n i is given by M , so that α can take values in the interval6 /M ≤ α ≤ Graph spectra.
We identify an unweighted graph through its adjacency matrix A ij , whose generic element ij is 1 if nodes i and j are connected, and 0 otherwise.This choice represents a simplification of the more complex case of a weightedgraph, where each link between i and j comes with a weight W ij . The spectrumof such an unweighted graph is given by the eigenvalues of A ij . Since we dealwith graphs that come in a unique connected component, the Perron-Frobeniustheorem ensures that the eigenvalue of largest modulus λ is unique, real andpositive. The spectral gap g is then defined as the difference in modulus between λ and the second eigenvalue λ . Here we deal with undirected HMNs, so that A ij is symmetric, and g = λ − λ . We can also compute the Laplacian spectrumof a graph, which consists of the eigenvalues of the graph Laplacian L ij = δ ij (cid:80) l k j − A ij —where k i = (cid:80) j A ij is the degree of node i — which providesa discretization of the Laplace-Beltrami operator [42]. If A ij is symmetric, L ij is symmetric and positive semi-definite, and the spectral gap g L is defined asthe modulus of its smallest nonzero eigenvalue. An alternative discretization ofthe Laplace-Beltrami operator is provided by the normalized Laplacian matrix L ij = δ ij − A ij / (cid:112) k i k j , for which the spectral gap g L can be computed as above.Finally, for completeness, let us also mention that the results for L also extendto the random walk Laplacian, which enters the master equation of randomwalks, and has the same spectrum of eigenvalues as L [42]. Dynamical model.
We simulate dynamics in HMN models of the brain using anull model of activity propagation in network —whose use has proven effectiveas a probing tool to understand paramount features of brain activity [30, 18];—in particular, we consider the, so-called, susceptible-infected-susceptible (SIS)dynamics. In this well-known model, nodes can be either active (infected I) orinactive (susceptible S); in terms of neural dynamics an active node correspondsto an active region in the brain, which can activate an inactive neighboring re-gion with a given probability κ , and which can be deactivated at rate µ (which7e set to 1 without loss of generality) due to e.g., exhaustion of synaptic re-sources.In particular, within the so-called quenched mean-field approximation, theequations for our model read [43] ddt ρ i ( t ) = − µρ i + κ [1 − ρ i ( t )] N (cid:88) j =1 A ij ρ j ( t ) , (1)where ρ i ( t ) is the probability that node i is in the active state at time t . Wecall ρ ( t ) = (cid:104) ρ i ( t ) (cid:105) the average activity at time t . In the general case, SISdynamics results in a critical value κ = κ c , above which activity invades thesystem indefinitely reaching a non-zero steady state ρ ( t → ∞ ) > κ < κ c (quasi-critical regime) constitute a Griffiths phase and are associated with rare-regioneffects and e.g. slow activity relaxation, pointing to an effective model fornormal brain function. We consider initial states where all nodes are activeand simulate time evolution for different values of κ . Following this protocol inHMNs, simulations are characterized by an initial transient regime, after which,the asymptotic behavior takes over, namely a steady state for κ > κ c and apower-law time-decay within a finite range of κ below κ c (a Griffiths phase; see[33, 34, 18]). Generation of functional data.
As a proxy for empirically measured FC net-works, in our computational study we determine co-activation matrices [44, 45].In particular, we acquired activity data from a time interval I in our simulationtime series ρ i ( t ). In particular, we choose I to occur after the initial transientregime described above (which can be simply determined by visual inspection).If κ > κ c , I is an interval of the steady state. If κ < κ c , I covers a part of thepower-law-decay regime instead. From the data in I , we generated a matrix C ij with generic ij element equal to the probability that nodes i and j are simulta-neously active in a given time bin ∆ t ( ρ i (∆ t ) = ρ j (∆ t ) = 1). Using these data,we averaged C ij across multiple realizations of the dynamics (at least 10 in allcases), and considered a threshold θ , thus generating an adjacency matrix of8he functional network for each θ by imposing A ij ( θ ) = H( θ − C ij ), where H( )denotes the Heaviside step function. Let us remark that, when recording thesesynthetic functional data, we focused on a single realization of the structuralnetwork hosting the dynamics, mimicking a real-life scenario in which diagnos-tics are conducted on a single subject.We us emphasize that the use of co-activation matrices as a proxy for FCis motivated by the fact that in the simplified dynamic model time series aresequences of binary states (1 and 0), which makes it natural to define simulta-neous activity. Experimental time series, instead, are not constituted by binarystates and require functional data to be computed as Pearson correlation coef-ficients. As both methods consist in computing activity correlations, we believethat our forthcoming conclusions for the synthetic case and its co-activation ma-trices carry over to real systems and correlation matrices. Furthermore, sincethe method proposed here requires only a finite time interval I , it finds a nat-ural application to the study of dynamical FC, where functional networks aregenerated from short time intervals. Mapping to percolation problem.
For every functional network characterized byan adjacency matrix A ij ( θ ), we varied θ continuously and we monitored howthe network structure evolves as measured by two indicators: • the size of the largest connected component s ; • the spectral gap g (as well as its Laplacian counterparts for completeness).For θ ≈ s = N and g = N ( g L = N , g L = 1), while for θ ≈ s = 1and g ≈ g L = g L = 0). Thus, for intermediate θ the network necessarilyundergoes a percolation-like phase transition, signalled by a decrease in s . Ifthe functional network is non-hierarchical, at the same transition point there isa change in terms of g , i.e. its becomes very small (see below). In other words,both s and g act as order parameters of the same percolation-like transitionand share the same critical value of θ . However, this scenario changes for hi-9rarchical networks. Since a hierarchical network must possess a small spectralgap without being fragmented, the two transitions must be well separated in-stead. In particular, we expect a threshold interval in which the spectral gap g decreases by orders of magnitude (the order of N ), while s remains closeto its maximum. Let us remark that the idea of analyzing the dependence of s on θ in functional networks is not new, but it was developed by Gallos andcollaborators [39]. However, here the focus is instead on a similar procedure forboth g and s and comparting the corresponding transition points. Experimental data.
Even if the present work relies on the generation of func-tional networks by simulating activity on synthetic structural networks, in or-der to achieve more general conclusions, we also extended the analyses to thecase in which the underlying structural networks are actual connectomes as de-rived from neuroimaging studies. In particular, we use the connectomes of twohealthy subjects, taken from the results of a broad experimental study, whichis described in the Appendix.
3. Results
Our main goal is to find an effective way to assess if a functional network ishierarchical, thus mimicking the topology of the underlying structural network,or if activity correlations result in a different, emergent topology. To this end, wedefine the concept of hierarchical functional network as follows. Since many ofthe desirable dynamic properties of hierarchical modular networks (rare regioneffects, generic criticality, localization [17, 18, 36, 20]) can be traced back to theirsmall spectral gaps g (as well as g L and g L ), we define a hierarchical functionalnetwork as a network that • i) is generated by activity patterns on a structural HMN; • ii) possesses vanishing spectral gaps for a finite range of θ values, whilestill being connected (i.e. not fragmented).10 .01 0.02 0.03 0.0410 o r d e r p a r a m e t e r a d e n s e h i e r a r c h i c a l f r a g m e n t e d s / Ng / N l a r g e s p e c t r a l g a p s b d e n s e f r a g m e n t e d s / Ng / N Figure 1: Persistence and breakdown of hierarchical organization. a) Largest connected com-ponent size s and spectral gap g in the quasi-critical regime ( κ ≤ κ c ), or Griffiths phase, forco-activation matrices A ij ( θ ), generated by activity propagation dynamics on a single struc-tural HMN. A range of threshold values θ exists, where the resulting functional network isconnected ( s = N ) and has vanishing spectral gaps (by up to three orders of magnitude, theupper limit given by the system size). More precisely, we construct the range as the interval[ θ , θ ], such that, given a fixed f < g ( θ ) /N = s ( θ ) /N = f . Here a we choose f = 0 . θ to signal the first significant drop in largest component size: s /N > f indicatesa network which is still nearly intact; s /N < f signals a significant reduction in giant com-ponent size, thus pointing to the onset of fragmentation. Different structural networks mayproduce different s ( θ ), requiring different choices of f . b) Same as in a) but for simulationsrun in the super-critical regime ( κ > κ c ). No significant drop in the spectral gap is recordedwhile the network is intact, pointing to the loss of hierarchical order. For as long as thenetwork is connected ( s = N ), g drops by negligible amounts, and remains of the order of N . Note that the second requirement ensures that the vanishing spectral gapproperty is not a trivial consequence of a network breaking into many connectedcomponents [42].To produce the functional networks we consider structural HMNs of size2 , which already result on dense co-activation matrices of 2 thus makingthe rest of the study computationally intensive. We note that most experi-mentally acquired connectomes belong in the same size range, allowing for adirect comparison as we will show later on. By running simulations on thesesynthetic HMNs and acquiring functional data as described in Section 2, we11bserved (see Figure 1a) that this definition of hierarchical functional networksdescribes exactly what happens in the quasi-critical regime. Indeed, simulationdata for the quasi-critical case show that, upon increasing the threshold θ , thenetwork undergoes a percolation-like transition (red dotted line), as the largestconnected component size s (red circles) diminishes from N (fully-connectednetwork) to 1 (fully-fragmented network). As we mentionaed in Section 2, ourmain focus is on the spectral gap g and on the way its dependence on θ differsfrom that of s . We find that the spectral gap g (blue squares) too exhibits atransition (blue dashed line), from N (when the functional network is dense andfully connected) to vanishing values. The two transitions however do not coin-cide and a range of control parameter values θ emerges, where the network isconnected ( s = N ), but spectral gaps vanish ( g → θ interval as theone that begins when g starts decreasing significantly and ends when s startsdecreasing by a comparable amount. Thus, within this range, the functionalnetwork exhibits a small spectral gap while still preserving its non-fragmentedstate; i.e. the network is hierarchical, accordingly to the criterion we proposed.The strength of this result resides in its dependence on the dynamical regimethat we simulate. We just showed that the functional network is hierarchical inthe quasi-critical regime, more precisely in the Griffiths phase. But, what hap-pens in the super-critical regime, obtained for activation rates κ > κ c , which isnormally associated with abnormal (excessive) neural activity? Figure 1b showsthat in this regime the two phase transitions almost coincide and, more impor-tantly, there is no range of θ values, where the functional network is connectedand the spectral gaps vanish. The functional network is not hierarchical, sincewhenever it is connected it has large spectral gaps (which only decrease slightly,still maintaining huge values, and not dropping by an amount comparable to N ).The picture of hierarchical functional networks that develop an emergent,non-hierarchical topology in the super-critical regime is further corroborated bythe study of the degree distributions of such networks, as shown in Figure 2.12 k P ( k ) k a = 0.028= 0.029= 0.030= 0.032= 0.033= 0.034 k P ( k ) k b = 0.092= 0.093= 0.094= 0.095= 0.096 Figure 2: Degree distributions of functional networks. a) Quasi-critical regime. Red curvesare obtained for values of θ producing networks with maximum separation between s and g .Larger values of θ (black curves) occur when the network is close to failure ( s → − . P ( k ) ∼ k − . isencountered for higher values of θ . .01 0.02 0.03 0.0410 o r d e r p a r a m e t e r d e n s e f r a g m e n t e d s / Ngg L / N Figure 3: Persistence of hierarchical organization, as from the study of Laplacian matrices.Data is obtained as in Figure 1a, employing spectral gaps from the graph Laplacian L andthe normalized Laplacian L . Data for s is the same as in Figure 1a and is reported here forcomparison. We choose f = 0 . In the quasi-critical phase, and for values of θ within the range highlighted inFigure 1a, degree distributions have power-law tails with continuously varyingexponents [46, 18]. In the super-critical regime instead, degree distributionsdisplay heavier power-law tails with an emergent, limiting exponent close to3 / g applied to theadjacency matrix A ij ( θ ). This choice is motivated by the fact that the func-tional network of adjacency matrix A ij ( θ ) is produced by an SIS process, whichis linearized by the adjacency matrix A ij of the underlying structural network.Nevertheless, one may wonder if the result above extends to Laplacian matrices[42], which are of interest in problems of transmission, diffusion and, more im-portantly in the case of brain dynamics, synchronization [36, 47, 48]. Figure 3shows that our main result carries over to the case of Laplacian matrices, bothin the non-normalized variant L (historically, the Kirchhoff Laplacian) and inthe normalized one L (sharing the same spectrum of eigenvalues as the randomwalk Laplacian). 14 igure 4: SC adjacency matrix of the human connectome. SC matrices of size N = 2514are averaged over 12 healthy patients. The hierarchical modular organization is highlightedby choosing a node labeling scheme based on agglomerative clustering methods [8]. Whilein our study we employ individual SC matrices, here we show a network average for ease ofvisualization. o r d e r p a r a m e t e r a d e n s e h i e r a r c h i c a l f r a g m e n t e d s / Ng / N b l a r g e s p e c t r a l g a p s f r a g m e n t e d s / Ng / N d e n s e o r d e r p a r a m e t e r c f r a g m e n t e dd e n s e h i e r a r c h i c a l s / Ng / N d d e n s e f r a g m e n t e d s / Ng / N l a r g e s p e c t r a l g a p s Figure 5: Persistence and breakdown of hierarchical organization of functional networks gen-erated from two real connectomes. (a) and (b): quasi-critical regime and super-critical regimefor the first subject. (c) and (d): quasi-critical regime and super-critical regime for the secondsubject. Results are as for the synthetic HMN case shown in Figure 1. Subject variabilityintroduces changes in the functional form of the curves for s and g . In particular, in orderto highlight the different intervals of θ , we fix f = 0 . f = 0 . s /N for which the networksstart fragmenting. Our main conclusions are the same as in Figure 1. In the quasi-criticalregime a range of θ emerges, where g decreases by over, at least, an order of magnitude whilethe network is still not fragmented. As for our numerical model, this behavior breaks downin the super-critical case, signaling a pathological state. N = 2514). The hierarchical modularstructure of human connectomes is shown in Figure 4, were we show the SCadjacency matrix, averaged over 12 healthy patients. We ran our simulationson SC adjacency matrices of individual subjects and performed our spectralanalysis of the resulting FC data. Our findings, shown in Figure 5 are strikingas they show how the same fingerprints of hierarchical functional organization inthe quasi-critical case emerge by simulating activity propagation simulations onactual connectomes. Data shown here is produced from the structural networkof two healthy subjects (top and bottom row respectively).As real connectomes encode much greater complexity than the HMN modelsused above, we notice significant changes in the shapes of the s ( θ ) functions,where an initial slow decrease in s anticipates the much faster drop, whichwe identify with the onset of the network fragmentation. We also notice signifi-cant quantitative differences introduced by subject variability. Notwithstandingthese differences, the conclusions of our analysis are the same as for the syntheticHMN case: for simulations run in the quasi-critical regime, one can highlight arange of θ in which s is still in its regime of slow decrease, while spectral gapshave dropped by up to two orders of magnitude. A systematic study of a muchlarger number of subjects is beyond the scope of this paper and is left for futurework.
4. Conclusion and discussion
The relationship between structure and function in brain networks remainsto this day a field of enormous interest. The ultimate goal of devising diagnos-tic tools that leverage concepts of network theory applied to SC connectivitydata begs for a deeper understanding of how dynamic patterns originate from17tructural motifs, or to which extent they self-organize irrespective of thosemotifs. What we introduced here is a minimal-model analysis, which howevercaptures some essential aspects of the structure-function relationship. In thequasi-critical regime, activity patterns are strongly localized and result in a FCnetwork which closely mimics the SC hierarchical patterns. Our criterion toestablish the similarity of SC and FC consists in highlighting the existence of anon-trivial threshold interval, in which FC networks exhibit the small-spectralgap property, known to characterize FC networks. Similarity thus consists inthe inheritance of a basic spectral property. This is not to say that activitycorrelations do not play a role in shaping FC patterns: the FC degree distribu-tions do not replicate those of SC, they in fact are power laws with continuouslyvarying exponents. This form of generic scale invariance has the same struc-tural origin as the one observed in avalanche size measurements in this sameregime [18, 49]: the hierarchical organization induces rare region effects, whichare reflected both in avalanches and in the topology of the functional network.In the super-critical regime, which represents a pathological state akin toepilepsy, the spectral fingerprint of hierarchical networks is lost and functionalconnectivity earns a global and much denser network topology. While it wasshown in the past that networks with localized structures host localized activitypatterns also in the super-critical phase [50] —something we were able to con-firm in HMN models of brain activity [21]— the change in spectral propertiesand degree distributions reported here points to a significant reorganization ofdynamics and correlations, in which modules mutually reactivateOur methodology addresses the analysis of activity correlation data withoutfixing an arbitrary threshold. Thresholds are treated as a control parameter in aphase diagram, which is analyzed as a whole, to show how, in the quasi-criticalregime, threshold values which correspond to non-trivial functional networktopologies (fully connected, fragmented) also lead to the non trivial property ofsmall spectral gaps, which allows us to conclude that the functional networkseffectively inherit the hierarchical organization of the structural ones.We believe that our work serves as a proof-of-concept study, for more com-18lex applications in which multiple spectral indicators are considered, and inwhich FC is inferred from real time series, rather than from numerical simu-lations. For instance, by altering the natural dynamics by drugs able to blockexcitation and/or inhibition, one could analyze how and when those alterationsare reflected in the hierarchical structure of FC networks and in its possiblebreakdown.The simple methodology we adopt here, based on co-activations, has clearadvantages while generating functional networks under minimal assumptions,and is motivated by the simplified nature of the dynamic model in use here.At the same time, our study can also be conducted on FC data extracted asPearson correlation, allowing for instance the separation of correlation and anti-correlation effects. Given the fact that both structural and functional data isavailable to us for a large number of subjects, as detailed in Section 2, a large-scale application of our methodology is planned an the next step of our investi-gation. In passing, we note that since the present results can be obtained withdata acquired in a short time window, our approach also lends itself naturallyto the study of experimental data sets obtained in the form of dynamical FC.Our work of course relies on essential simplifications of the original neuro-science problem. Namely our simulated dynamics is a simplification of neuralactivation processes. Nevertheless, we believe that the results reported hereprovide useful insights regarding the extent to which structure and function areconnected in brain networks, and how pathological states can be accompaniedby radical shifts in functional connectivity. We believe that the type of spec-tral analyses presented here to study persistent and emergent topologies in FCdata sheds light on the relationship between optimal architectures and tuneableperformance in the broader context of neurocomputing.
References [1] Sporns O, Chialvo DR, Kaiser M, Hilgetag CC. Organization, develop-ment and function of complex brain networks. Trends in cognitive sciences19004;8:418–25.[2] Bullmore E, Barnes A, Bassett DS, Fornito A, Kitzbichler M, MeunierD, et al. Generic aspects of complexity in brain imaging data and otherbiological systems. Neuroimage 2009;47:1125–34.[3] Sporns O. Networks of the brain. MIT Press 2010;.[4] Zhou C, Zemanov´a L, Zamora G, Hilgetag CC, Kurths J. Hierarchicalorganization unveiled by functional connectivity in complex brain networks.Phys Rev Lett 2006;97:238103. doi: .[5] M¨uller-Linow M, Hilgetag CC, H¨utt MT. Organization of excitable dynam-ics in hierarchical biological networks. PLoS Comput Biol 2008;4:e1000190.[6] Bullmore E, Sporns O. The economy of brain network organization. NatureRev Neurosci 2012;13:336–49.[7] Sporns O. Structure and function of complex brain networks. DialoguesClin Neurosci 2013;15(3):247–62.[8] Diez I, Bonifazi P, Escudero I, Mateos B, Mu˜noz MA, Stramaglia S, et al.A novel brain partition highlights the modular skeleton shared by structureand function. Scientific reports 2015;5:10532.[9] Craddock R, Jbabdi S, Yan C, Vogelstein J, Castellanos F, Martino AD,et al. Imaging human connectomes at the macroscale. Nat Methods2013;10:524–39. doi: .[10] Fornito A, Zalesky A, Bullmore E. Fundamentals of Brain Network Anal-ysis. Academic Press; 2016. ISBN 9780124081185.[11] Zhan L, Jenkins LM, Wolfson OE, GadElkarim JJ, Nocito K, ThompsonPM, et al. The significance of negative correlations in brain connectivity.Journal of Comparative Neurology 2017;525(15):3251–65.2012] Biswal BB, Mennes M, Zuo XN, Gohel S, Kelly C, Smith SM, et al. Towarddiscovery science of human brain function. Proceedings of the NationalAcademy of Sciences 2010;107(10):4734–9. doi: .[13] Fox M, Greicius M. Clinical applications of resting state functional connec-tivity. Frontiers in Systems Neuroscience 2010;4:19. doi: .[14] Essen DV, Ugurbil K, Auerbach E, Barch D, Behrens T, Bucholz R, et al.The human connectome project: A data acquisition perspective. NeuroIm-age 2012;62(4):2222 –31. doi: https://doi.org/10.1016/j.neuroimage.2012.02.018 .[15] Lee M, Smyser C, Shimony J. Resting-state fmri: A review of methodsand clinical applications. AJNR Am J Neuroradiol 2013;34(10):1866–72.doi: .[16] Haimovici A, Tagliazucchi E, Balenzuela P, Chialvo DR. Brain organizationinto resting state networks emerges at criticality on a model of the humanconnectome. Physical Review Letters 2013;110:178101.[17] Friedman EJ, Landsberg AS. Hierarchical networks, power laws, and neu-ronal avalanches. Chaos 2013;23:013135.[18] Moretti P, Mu˜noz MA. Griffiths phases and the stretching of criticality inbrain networks. Nature Commun 2013;4:2521.[19] ´Odor G. Localization transition, lifschitz tails, and rare-region effects innetwork models. Phys Rev E 2014;90:032110.[20] ´Odor G, Dickman R, ´Odor G. Griffiths phases and localization in hierar-chical modular networks. Sci Rep 2015;5:14451.[21] Safari A, Moretti P, Mu˜noz MA. Topological dimension tunes activ-ity patterns in hierarchical modular networks. New Journal of Physics21017;19(11):113011. URL: https://doi.org/10.1088%2F1367-2630%2Faa823e . doi: .[22] Allen EA, Damaraju E, Plis SM, Erhardt EB, Eichele T, Calhoun VD.Tracking whole-brain connectivity dynamics in the resting state. Cerebralcortex 2014;24:663–76.[23] Cabral J, Kringelbach ML, Deco G. Functional connectivity dynami-cally evolves on multiple time-scales over a static structural connectome:Models and mechanisms. NeuroImage 2017;160:84–96. URL: https://linkinghub.elsevier.com/retrieve/pii/S1053811917302537 . doi: .[24] Preti MG, Bolton TA, Van De Ville D. The dynamic func-tional connectome: State-of-the-art and perspectives. NeuroImage2017;160:41–54. URL: https://linkinghub.elsevier.com/retrieve/pii/S1053811916307881 . doi: .[25] Sporns O, Tononi G, K¨otter R. The human connectome: A structuraldescription of the human brain. PLoS Comput Biol 2005;1(4):e42.[26] Hagmann P, Cammoun L, Gigandet X, Meuli R, Honey CJ, Wedeen VJ,et al. Mapping the structural core of human cerebral cortex. PLoS Biol2008;6:e159.[27] Kaiser M. A tutorial in connectome analysis: topological and spatial fea-tures of brain networks. Neuroimage 2011;57(3):892–907.[28] Dorogovtsev SN, Goltsev AV, Mendes JFF. Critical phenomena in complexnetworks. Rev Mod Phys 2008;80:1275.[29] Boettcher S, Cook JL, Ziff RM. Patchy percolation on a hierarchical net-work with small-world bonds. Phys Rev E 2009;80:041115.[30] Kaiser M, G¨orner M, Hilgetag CC. Criticality of spreading dynamics inhierarchical cluster networks without inhibition. New J Phys 2007;9:110.2231] Beggs JM, Plenz D. Neuronal avalanches in neocortical circuits. J Neurosci2003;23:11167–11177.[32] Vojta T. Rare region effects at classical, quantum and nonequilibriumphase transitions. J Phys A 2006;39:R143.[33] Mu˜noz MA, Juh´asz R, Castellano C, ´Odor G. Griffiths phases on complexnetworks. Phys Rev Lett 2010;105(12):128701.[34] Juh´asz R, ´Odor G, Castellano C, Mu˜noz MA. Rare-region effects in thecontact process on networks. Phys Rev E 2012;85:066125.[35] Villa Mart´ın P, Moretti P, Mu˜noz MA. Rounding of abrupt phase transi-tions in brain networks. J Stat Mech 2015;:P01003.[36] Villegas P, Moretti P, Mu˜noz MA. Frustrated hierarchical synchroniza-tion and emergent complexity in the human connectome network. Sci Rep2014;4:5990.[37] Mill´an AP, Torres J, Bianconi G. Complex network geometry and frus-trated synchronization. Sci Rep 2018;8:9910.[38] Mill´an AP, Torres JJ, Bianconi G. Synchronization in network geometrieswith finite spectral dimension. Physical Review E 2019;99(2):022307.[39] Gallos LK, Makse HA, Sigman M. A small world of weak ties provides opti-mal global integration of self-similar modules in functional brain networks.Proc Natl Acad Sci USA 2012;109:2825–30.[40] Farkas IJ, Der´enyi I. Barab´asi AL, Vicsek T. Spectra of “real-world”graphs: Beyond the semicircle law. Phys Rev E 2001;64:26704.[41] Wang SJ, Hilgetag CC, Zhou C. Sustained activity in hierarchical modularneural networks: Self-organized criticality and oscillations. Front ComputNeurosci 2011;5:30. 2342] Chung FRK. Spectral graph theory. CBMS Regional Conference Series inMathematics 1997;92:212.[43] Pastor-Satorras R, Castellano C, Van Mieghem P, Vespignani A. Epidemicprocesses in complex networks. Rev Mod Phys 2015;87:925–79.[44] H¨utt MT, Kaiser M, Hilgetag CC. Perspective: network-guided patternformation of neural dynamics. Philosophical Transactions of the Royal So-ciety B: Biological Sciences 2014;369(1653):20130522. doi: .[45] Damicelli F, Hilgetag CC, H¨utt MT, Mess´e A. Topological reinforcementas a principle of modularity emergence in brain networks. Network Neuro-science 2019;3(2):589–605. doi: .[46] Clauset A, Shalizi CR, Newman MEJ. Power-law distributions in empiricaldata. SIAM Rev 2009;51(4):661–703.[47] Di Santo S, Villegas P, Burioni R, Mu˜noz MA. Landau–ginzburg theory ofcortex dynamics: Scale-free avalanches emerge at the edge of synchroniza-tion. Proceedings of the National Academy of Sciences 2018;115(7):E1356–65.[48] Buend´ıa V, Villegas P, Burioni R, Mu˜noz MA. Hybrid-type synchronizationtransitions: where marginal coherence, scale-free avalanches, and bistabilitylive together. arXiv preprint arXiv:201103263 2020;.[49] ´Odor G. Critical dynamics on a large human open connectome network.Phys Rev E 2016;94:062411.[50] Goltsev AV, Dorogovtsev SN, Oliveira JG, Mendes JFF. Localization andspreading of diseases in complex networks. Phys Rev Lett 2012;109:128702.[51] Graben PB, Jimenez-Marin A, Diez I, Cortes J, Desroches M, RodriguesS. Metastable resting state brain dynamics. Front Comput Neurosci2019;13:62. 2452] Cook P, Bai Y, Nedjati-Gilani S, Seunarine K, Hall M, Parker G, et al.Camino: Open-Source Diffusion-MRI Reconstruction and Processing. 14thScientific Meeting of the International Society for Magnetic Resonance inMedicine 2006;:2759.
Appendix
Data available to us are from 30 healthy subjects (14 males, 16 females)with age between 22 and 35. Data were provided through the Human Con-nectome Project, WU-Minn Consortium (Principal Investigators: David VanEssen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutesand Centers that support the NIH Blueprint for Neuroscience Research; andby the McDonnell Center for Systems Neuroscience at Washington University.To build the connectivity matrices (for further details see [51]), we processedthe same-subject structure-function triple acquisitions of magnetic resonanceimaging (MRI) – see Appendix for details on acquisition parameters – consist-ing of: 1. High-resolution anatomical MRI (used for the mask segmentation ofgray matter, white matter and cerebrospinal fluid, and for the transformationto common-space of the functional and the diffusion data), 2. Functional MRIat rest (used for extracting region time series of the blood-oxygen-dependentsignal, after removal of movement artifacts and physiological noise, but not theglobal signal regression), and 3. Diffusion MRI (used for building SC matricesafter fitting a diffusion tensor to each voxel, running a deterministic tractogra-phy algorithm using the UCL Camino Diffusion MRI Toolkit [52], and countingthe number of streamlines connecting all pairs within the N = 2514 regions,each one containing on average 66 voxels. Imaging acquisition parameters
Same subject structure-function triple acquisitions were performed using a3T Siemens Connectome Skyra with a 100 mT/m and 32-channel receive coils.25he acquisition consisted of : High-resolution anatomical MRI was acquired using a T1-weighted 3D MPRAGEsequence with the following parameters: TR=2400ms; TE=2.14ms; Flip an-gle=8deg; FOV=224 ×
224 mm ; Voxel size=0.7mm isotropic; Acquisition time= 7 min and 40 sec. Functional data at rest were acquired to obtain the blood-oxygenation-level-dependent (BOLD) signals with a gradient-echo EPI sequence with the followingparameters: TR = 720ms, TE = 33.1ms; Flip angle = 52 deg; FOV = 208 ×
180 mm ; Matrix = 104 ×
90; 72 slices per volume, a total number of 1200volumes; Voxel size = 2mm isotropic; Acquisition time = 14 min and 33 sec.
Diffusion data were acquired with a Spin-echo EPI sequence with the follow-ing parameters: TR = 5520ms; TE = 89.5ms; Flip angle = 78 deg; FOV = 210 ×
180 mm ; Matrix = 168 × Abbreviations
FC: functional connectivity HMN: hierarchical modular network MRI: mag-netic resonance imaging SC: structural connectivity SIS: susceptible-infected-susceptible
Competing interests
The authors declare that they have no competing interests. Here, we only report the parameters which are neeeded for the image preprocessing,but a complete information about the acquisitions can be found at http://protocols.humanconnectome.org/HCP/3T/imaging-protocols.html thics approval and consent to participate The sample included 30 subjects from the MGH-USC Human ConnectomeProject. The research was performed in compliance with the Code of Ethics ofthe World Medical Association (Declaration of Helsinki). All subjects providedwritten informed consent, approved by the ethics committee in accordance withguidelines of HCP WU-Minn.
Funding
We acknowledge the Spanish Ministry and Agencia Estatal de investigaci´on(AEI) through grant FIS2017-84256-P (European Regional Development Fund),as well as the Consejer´ıa de Conocimiento, Investigaci´on Universidad, Juntade Andaluc´ıa and European Regional Development Fund, Ref. A-FQM-175-UGR18 and SOMM17/6105/UGR for financial support. AS and PM acknowl-edge financial support from the Deutsche Forschungsgemeinschaft, under grantsMO 3049/1-1 and MO 3049/3-1.
Availability of data and materials
Numerical simulation data are available on request from the correspondingauthor.
Author contributions
AS: Performed the simulations, Analyzed the data, Wrote the manuscript;PM: Designed the study, Analyzed the data, Wrote the manuscript; ID: Datacuration, Preprocessed the images; JMC: Data curation, Wrote the manuscript.MAM: Designed the study, Analyzed the data, Wrote the manuscript