How structure sculpts function: unveiling the contribution of anatomical connectivity to the brain's spontaneous correlation structure

Intrinsic brain activity is characterized by highly structured co-activations between different regions, whose origin is still under debate. In this paper, we address the question whether it is possible to unveil how the underlying anatomical connectivity shape the brain's spontaneous correlation structure. We start from the assumption that in order for two nodes to exhibit large covariation, they must be exposed to similar input patterns from the entire network. We then acknowledge that information rarely spreads only along an unique route, but rather travels along all possible paths. In real networks the strength of local perturbations tends to decay as they propagate away from the sources, leading to a progressive attenuation of the original information content and, thus, of their influence. We use these notions to derive a novel analytical measure, $\mathcal{T}$ , which quantifies the similarity of the whole-network input patterns arriving at any two nodes only due to the underlying topology, in what is a generalization of the matching index. We show that this measure of topological similarity can indeed be used to predict the contribution of network topology to the expected correlation structure, thus unveiling the mechanism behind the tight but elusive relationship between structure and function in complex networks. Finally, we use this measure to investigate brain connectivity, showing that information about the topology defined by the complex fabric of brain axonal pathways specifies to a large extent the time-average functional connectivity observed at rest.


I. INTRODUCTION
In the last two decades a large body of research has demonstrated that spontaneous brain activity forms structured patterns of consistent co-activations across different subsets of brain regions [1][2][3][4][5][6][7][8][9][10] . Contrary to what was somehow implicitly assumed, intrinsic brain activity cannot be considered as a simple sum of rather unpredictable and noisy fluctuations. Spontaneous collective dynamics of different brain regions, measured either with EEG, fMRI or MEG, can be clustered into highly organized and reproducible spatial patterns, referred to as resting-state networks (RSNs), that strikingly resemble those activations observed during the performance of different tasks 7,11,12 . Furthermore, growing evidences suggest that this spontaneous large-scale structure is characterized by a marked temporal organization, mirrored by recurrent alternations of subnetworks that concur in generating a rich dynamical repertoire at different time scales 13,14 . Although the origin and the purpose of spontaneous brain dynamics are still under debate, there is wide agreement that mental states and brain malfunction alter the patterns of spontaneous activity, e.g. the a) Electronic mail: ruggero.bettinardi@upf.edu b) Electronic mail: gorka@Zamora-Lopez.xyz dynamical repertoire of the brain at rest tends to decrease during sleep 15 and under anaesthesia [16][17][18] .
From its discovery, great efforts have been invested in trying to reproduce resting-state brain activity through the use of computational models, in order to obtain a mechanistic explanation of this intriguing but elusive phenomenon. Early models based on the structural connectomes of cats and macaques extensively explored the emerging patterns of correlations in those networks at different spatial and temporal scales [19][20][21] . With the arrival of structural human connectomes obtained through tractography, computational models could finally attempt to fit the empirical correlation structure observed via resting-state fMRI [22][23][24][25][26] . Despite these attempts, we still lack of a precise understanding of the relationship between the shape of the brain's connectome and the emergent patterns of correlations observed during rest. One of the major reasons is that the collective dynamics of a network do not only depend on the shape of the underlying connectivity, but also on the model chosen to simulate the local dynamics of the cortical regions 26,27 .
In this paper we aim at unveiling what is the precise contribution of the anatomical connectivity on the correlation structure observed during rest. To this aim, we first need to derive a theoretical estimate of the expected correlation between nodes due to the network's topology alone; therefore we will try to avoid, as much as possible, the contribution of other factors. To do so, in Section II we introduce a novel graph theoretical quantity measuring the topological similarity of the entire input profiles that two nodes receive from the whole network. This measure, that we named topological similarity, is a generalization of the concept of matching index that explicitly accounts for the fact that in networks information travels along all possible paths, not only along the shortest ones, and that information content tends to decay as it moves away from its source [28][29][30][31] . This measure is based on the concept of network's communicability introduced by Estrada and Hatano 32 , a function that quantifies the strength of the influence that one node exerts over another through all paths of any length assuming an exponential decay of influence with path length. We then demonstrate that the topological similarity function, altough based on pure topological information about the underlying path structure of the network, can indeed be used to approximate the expected correlation structure of networks of time-varying coupled units.
In Section III we systematically investigate the contribution of three topological primitives: the weight of the links, the length of the path and the presence of redundant alternative paths through which information can travel. Topological primitives are fundamental features underlying the network architecture which determine how influence spreads, thus sculpting the similarity of the input profiles of nodes. To do so, we study simple graphs in which these primitive features can be manipulated. We show how these three features are in fact of primary importance in modulating the strength of the influence and the relative topologies into which different nodes are embedded.
Finally, in Section IV we investigate the spontaneous correlation structure of the human brain, demonstrating that taking into account the similarity of the wholenetwork influences shaped by the underlying anatomical topology it is possible to understand and predict the large-scale functional organization observed during rest.

II. HOW TOPOLOGY SCULPTS THE CORRELATION STRUCTURE OF NETWORKS
In order to properly address the contribution of network's topology on the emerging spontaneous correlation structure, we first need to acknowledge that the collective behaviour of a set of coupled dynamical units depends on three principal ingredients: (i) the structure of the network, (ii) the local dynamics of the nodes and (iii) the coupling function determining how information is passed from one node to another. In fact, for a fixed network, changing the local dynamical model of the nodes and the coupling function usually leads to different collective dynamics 26,27,33 . Therefore, in order to estimate the contribution of structure alone we need to set apart the role of the other two factors.
Typically, the activity of two nodes exhibits statistical dependence either if they are connected by means of a direct link, or if the aggregate inputs they receive from the entire network are similar, independently of whether there is a link between the two or not. Because information in a network rarely travels exclusively along the shortest paths 34,35 but instead diffuses along the whole network, we realise that the total influence of one node over another mainly depends on three topological features: (i) the strength of the coupling between them, usually represented by the weights of the links, (ii) the graph distance between the two nodes and (iii) the propagation of the influence along multiple alternative paths 36,37 . We refer to these three topological features as the topological primitives because, in combination, they encode the shape of the network's path structure. As such the analysis of the building blocks underlying larger complex networks is a necessary step to understand the contribution of network structure to the emergent activity.
In general, the influence of a direct link is greater than the influence exerted over longer paths, as the latter is mediated through third nodes. In fact, in real systems the "power" of the signals or their information content naturally decays along the path 28-31 , unless there exists an active mechanism which amplifies the incoming signal at the cost of energy. Following the same rationale, it is unlikely that influence or information from one node to another propagates only along a single, selected path, unless there are specific gating mechanisms controlling the detailed routing of information over all existing paths. The total number of paths (of all lengths) between any two nodes in a network is in fact infinite. It is well known that the total number of paths m l ij of length l between nodes i and j in a graph grows with l. This number is given exactly by the l th power of the adjacency matrix A, m ij = (A l ) ij 38,39 . Thus, the total number of paths leaving from node i and arriving at node j is given by the sum: This number typically diverges and thus, for the dynamics within a network to remain bounded, the amount of influence needs to decay faster with the length than the growth in the number of paths. Mathematically, the problem consists in finding a set of coefficients {k l } for which the series ∞ l=0 k l A l converges for any adjacency matrix, A. While the solution to this problem is not unique, Estrada and Hatano 32 proposed an exponential decay of the influence with path length and introduced the communicability measure C. The communicability function thus corresponds to the matrix exponential of A, which can be expanded into a series of powers with coefficients k l = 1/ l!: From a physical perspective, the communicability is analogous to the Green's function of the network 32,40 and expresses how local perturbations propagate along the system. Perturbation of a given node's dynamics in fact first propagates to its direct neighbours, affecting their activity; the activity of the neighbours will in turn propagate to their neighbours (as well as back to the node that was perturbed in the first place), and interact with their intrinsic local dynamics. This simple propagation mechanism implies that the effect of a local perturbation could be perceived also by distant nodes but attenuated and modulated by the dynamics of each node along the route. Due to its correspondence to the Green's function, the communicability can be tuned using a constant global coupling parameter that uniformly scales the weights of all links in A 40,41 , allowing to search through the emerging collective dynamics over multiple scales. The generalised, tunable, communicability is then: When g is weak, perturbations quickly decay, producing local correlations only around the node's neighbourhood. As g grows perturbations propagate deeper into the network, giving raise to stronger correlations over more distant nodes. Zamora-López et al. 41 have shown that considering the communicability as the propagator kernel for the diffusion of Gaussian noise along the network, it reveals an equivalent correlation structure as networks of generic and widely used models, e.g. Kuramoto oscillators and neural masses (see Supplementary Information in Ref. 41). Under these assumptions (Gaussian white noise sources and exponential decay), in Ref 41 it was shown that the covariance matrix of the system, Σ, can be analytically estimated as Σ = C · C T where C T is the transposed communicability. The cross-correlation matrix R of the system is then calculated, as usual, normalising the rows of the covariances Σ ij by the autocovariances Σ ii as, R ij = Σ ij /Σ ii . Despite the merit of being analytical, this estimate of the network's crosscorrelation matrix still relied on assuming very simple dynamical model, the diffusion of Gaussian white noise, and did not fully disentangle the unique contribution of the topology. In the present work we want to address the question whether it is possible to estimate the most likely correlation structure that a network gives rise, based only on its topological properties. As mentioned before, it is legitimate to assume that the activities of two nodes will exhibit statistical dependence either if they are connected by means of a direct link, or if the aggregate input they receive from the entire network is similar. With this in mind, we realised that the column vectors c j of the communicability matrix C indeed represent the input profile of the influences node j receives from all nodes in the network along all possible paths, including the influence a node has on itself due to recurrent paths. Therefore, it should be possible to predict analytically the magnitude of the expected correlation between any two nodes in a network by comparing the similarity of their input profiles. The resemblance between two input profiles can in principle be calculated using any measure of similarity between multidimensional vectors. Here we choose the cosine similarity because it returns results bounded between −1 and 1, equivalent to the cross-correlation measure we aim at comparing. Consequently, we define the topological similarity, T ij , between two nodes as the cosine of the angle between the corresponding columns c i and c j in the communicability matrix: being , the inner product and the vector norm. The definition of T depends uniquely on the topological constraints of the network encoded in the adjacency matrix A, plus the realistic assumption (embedded in C) that the influence or the information content decays with the length of the path. Despite being a purely topological measure, we recognized that, when the variance of the nodes is the same for all nodes, T corresponds exactly to the correlation matrix R of the Gaussian diffusion dynamics studied in Ref 41. Therefore, T formally closes the cycle for the search of a direct relation between the structure of a network and the expected correlation pattern that this structure will tend to generate. Furthermore, this analytical measure unveils the fundamental mechanism behind the contribution of network structure to the emerging function, which can be understood in terms of the similarity of influences that two nodes receive from the whole network via its complete path structure.
Finally, we shall notice that T can be estimated for both directed and undirected graphs, as well as for weighted networks. In the case of undirected graphs C is symmetric but if the links are directed, then the columns of C determine the input profiles and the rows represent the profile of output influences of the nodes. It must be noted that despite the measure can be computed for any weighted adjacency matrix, it does not always make sense to do so. Because communicability is a measure of influence along the paths, it only has a direct physical meaning when the weights of the links represent the coupling strength between the nodes, the flow capacity of the link or a compatible physical sense.
Summarising, here we have introduced a topological estimator of a network's correlation structure. Although it ignores any specific dynamics of the nodes, it accounts for the fact that information or influence within a net-work propagates along all possible alternative paths and that naturally decays for longer paths. The measure actually quantifies the similarity between the input profiles of nodes. In the following section we systematically investigate how three fundamental features of a network (the weights of the links, the path length and the redundancy of paths) influence both the communicability and the topological similarity between nodes that, as shown above, can in fact be used as a proxy of their expected correlation.

III. EFFECT OF TOPOLOGICAL PRIMITIVES
In this section, we will investigate how three topological primitives we defined before, namely links' weights, graph distance and the presence of multiple alternative paths, sculpts the influence that one node exerts over another, setting aside the role of local nodes' dynamics. To this aim, we will focus on three simple classes of graphs: (i) chains, (ii) cycles and (iii) path-redundant networks. We will evaluate how manipulating critical parameters of these graphs leads to changes in the influence between given pairs of nodes a and b (measured by their communicability, C ab ) and in their topological similarity T ab .
The contribution of link weight in the simplest case can be understood when analyzing how changes in the weight modulate the mutual influence of two nodes directly connected by a single link (see Supplemental Figure S1). As expected, incresing link weight is associated with increase in both the influence that one node exert over the other, as well as in their topological similarity.
A simple manner to obtain a better intuition of the behavior of both communicability and topological similarity is by studying how increasing the graph distance between two reference nodes affects both C ab and T ab in simple chain topologies. See top-left panels of Figure 1.
In general, it is possible to see how increasing the length of the path separating the two nodes correspond to a decrease in both their communicability and topological similarity, behavior directly determined by the decay formalized in the definition of the communicability. From the example, it is indeed evident how, for increasing lenghts of the chain, the input profiles of the two nodes at the ends of the chain (highlighted with green rectangles in the communicability matrices, top-left panels of Figure 1) become more and more antithetic, reflecting opposed whole-network influences. This behavior is captured by the corresponding decrease in the nodes topological similarity. The effect of uniformly varying the weights of all links in a chain is mainly quantitative: in fact, augmenting the weights in chains of fixed length increases the strength of the influence of each node over all other nodes in a way that is inversely proportional to the distance separating them, top-right panel of  Figure S2 to appreciate the interaction of chain length and links' weight on C and T between the nodes at the two ends of a chain).
Due to the definition of communicability, the length of the path separating any two nodes as well as the weights of the links are the most important parameters defining the strength of their mutual influence and therefore their topological similarity as well. As such, chain topologies can be thought as a baseline to compare how more complex motifs such as cycles and path-redundant topologies, which indeed incorporate chain topologies, modulate both the resulting communicability and the topological similarity between given nodes.
With this in mind, we compared C ab and T ab of the three model graphs (chains, cycles and path-redundant motifs) having corresponding graph diameter, i.e. having same longest path L. Bottom-left panel of Figure 1 provides a schematic representation of different motifs having equivalent longest paths. For the case of chains and path-redundant topologies, we computed C ab and T ab for those two nodes that were more distant, in other words, those at the two extremes of (each) path, whereas for cyclic topologies, we always selected two adjacent nodes: this choice allowed us to clearly disentangle the contribution that the direct link has on both the commmunicability and on the topological similarity, above and beyond the modulation produced by chains of different length.
As expected, in chains both C ab and T ab decay as the distance between the extremal nodes increases, see bottom-right panels of Figure 1, blue lines. On the other hand the effect of the direct link is well illustrated for the case of cyclic architectures (bottom-right panels in Figure 1, red lines). In fact, the presence of a direct link importantly enhances and poses a lower bound for both C ab and T ab , that however exhibit a chain-like decay as the indirect path between them increases.
This approach makes the difference between chains and cyclic primitives straightforward; however, in order to being able to properly understand how increasing path redundancy affects the communicability and hence topological similarity, it can be of great help to directly quantify the difference between C ab and T ab calculated from chains and from redundant topologies for comparable path length. The two lowest-right panels in Figure 1 illustrate three examples of these differences, namely the cases with two (dark green lines), three (light green lines) and four (orange lines) redundant paths. From this analysis, we appreciate that increasing the number of alternative paths does increase both the total influence and the topological similarity between the two reference nodes at both ends of the paths, but that the magnitude of this increase decays with the length of the path, vanishing for paths longer than 5 links.
Studying these simple networks (single link, chain, cycles and redundant paths) gives the opportunity to understand the central relevance of the three topological primitives in shaping how the influences of nodes unfould through the graph. This information in turn determines how strong will be the expected correlation between any pair of nodes, approximated by their topological similarity. These simplified network models can thus be used to summarize some of the simplest forms of interactions sustained by these topological features.

IV. UNDERSTANDING THE BRAIN'S SPONTANEOUS CORRELATION STRUCTURE
In the previous section we have analyzed how simple network architectures could affect both communicability and topological similarity between a given pair of nodes. However, real networks are made of interwined assemblies of those topological motifs that form intricate architectures and, together with the particular dynamical properties characterizing the system at hand, determine the emergence of complex patterns of interactions. In this section, we will try to find out) how much of the complex pattern of spontaneous correlations empirically observed in the resting brain can be explained just by the topology of the underlying anatomical structure.
The main assumption we make is the following: if two brain regions receive similar influences from the entire network, then the probability that they will exhibit consistent co-activations is high; on the other hand, if the inputs they receive from the whole network are very different, it is more likely that their covariance is weak. In Sections II and III, we showed how the the structure of a network can be used to estimate the strength of the influence that one node exerts over another trough the communicability C, and that C can be further used to quantify the topological similarity, T , of the input profiles between any pair of nodes. We then demonstrated that T can in fact be interpreted as the expected correlation between the nodes in the network. As such, we will use T to estimate the unique contribution of topology to the empirical time-average correlation structure of the brain's spontaneous activity measured using resting-state fMRI.
To this aim, we will compute the topological similarity matrix from the group-average structural connectivity matrix (SC), and use it to estimate the correlation structure mirrored by the group-average empirical functional connectivity matrix (FC) obtained from restingstate fMRI. See Section VI for details about the empirical SC and FC matrices. The SC matrix stores the information about axonal pathways reconstructed using whole-brain diffusion tensor imaging and tractography, thus defines the whole-brain anatomical wiring diagram of the brain. It should be noted that, despite their reproducibility, current methods used to reconstruct fiber bundles from diffusion imaging are characterized by intrinsic limitations constraining their accuracy: in fact, it is well-known that these reconstruction algorithms tend to favor the shortest, straightest and simplest path between any two reference voxels 42 , which in turn impair their ability to accurately detect crossing fibers and long inter-hemispheric axons 43,44 . With this in mind, we focused only on intra-hemipsheric structural and functional connectivity, in order to avoid the confounding effects of accumulating errors due to the above mentioned limitations in sampling inter-hemispheric pathways.
As mentioned in Section II, the communicability can be scaled by a factor g controlling the emerging collective dynamics, thus leading to the possibility of obtaining as many topological similarity matrices as the values of g used to scale C (see Figure 2 for three instances of T obtained for increasing values of g).
For each hemisphere, we thus searched for the topological similarity matrix T that best explained the observed intra-hemispheric correlation structure, by optimizing the communicability matrix according to a global scaling factor, g (Panel G and H in Figure 3).
As a mesure of model fitting, we used the mean absolute error (MAE), an outlier-robust alternative of the mean squared error (MSE), a classic statistic to quantify the goodness of an estimator. From the scatter plots (Panels I,L of Figure 3) and the best-fitting T matrices (Panels M,N of Figure 3), it is possible to appreciate how properly accounting for the overall input pattern sustained by a given topology can indeed reduce the error in estimating the emergent correlations from the raw structural connectivity alone; in fact, in the example we pass from a E(SC,FC)≈ 0.42 to E(T ,FC)≈ 0.15).
These results demonstrate that knowledge of the topology of whole-network input patterns of different brain regions, sustained by direct and indirect routes of multiple interweaved axonal bundles, can very much predict the time-average correlation structure observed from spontaneous BOLD fluctuations, above and beyond the information about direct anatomical connections stored in the SC matrices.
In the following, we will investigate how explicitly introducing local nodes' dynamics could affect the prediction of the average correlations. To this aim, we will make use of the Hopf normal model, that have successfully used to predict mesoscopic brain activity 45,46 , and evaluate their ability to explain brain functional connectivity. By this comparison, we will be able to gain some insight about the contribution of adding local dynamics in determining the observed correlation structure, as well as to have a better understanding of the role of the underlying topology, shared either by the two models.

A. Introducing local dynamics
We will now explicitly introduce local dynamics to simulate large-scale brain activity using the connectional architecture defined by the empirical SC matrix. By doing this, we will be able to estimate, through numerical simulations, the correlation between different brain regions, and then compared these results with the average empirical correlation matrix (Empirical FC) obtained from resting-state fMRI. Numerical simulations will be achieved through the use of the Hopf normal model, a dy-namical model able to reconcile noise-based approaches with models based on oscillators. This formalism is based on the normal form of a Hopf bifurcation 45,47 , a type of bifurcation that occurs when a system characterized by a stable fixed point loses its stability by exhibiting oscillations. As such, this model allows transitions between asynchronous noise activity and oscillations, thus making it a good candidate to reproduce empirical data as observed either with EEG, MEG or fMRI [45][46][47] . This model rely on the choice of two parameters, namely g, the global scaling of the strength of all the links in the network, and α, the parameter controlling the dynamical regime of each node. The bifurcation parameter will be set to 0, meaning that all nodes lie at the bifurca- tion working point, region that has been demonstrated to give good approximate of the empirical resting-state FC 46 , and we will thus optimize for g. Results from numerical simulations are depicted in Figure 4. The introduction of complex local dynamics does indeed increase the capacity of predicting the time-average correlation structure (MAE(T )≈ 0.15, whereas MAE(Hopf)≈ 0.11)), even though the magnitude of this increase is relatively small, suggesting that most of the average structure observed in correlated spontaneous BOLD fluctuations are in fact largely determined by the topology of the underlying network's architecture.

V. SUMMARY AND DISCUSSION
The brain is a complex system and as such its overall dynamics cannot be fully understood without taking into account the rich patterns of interactions into which its components are inherently embedded into. The collective dynamics in a network results from the complex interplay between its underlying structure, the nature of the local dynamics of individual nodes, their specific working point and the manner in which they are coupled 26,27,33 . Here, we aimed at understand to what extent does the structure of a network drive the emerging patterns of interactions. Is it possible, knowing the complete wiring diagram of a network, to estimate its most likely correlation structure?
Setting aside the large influence from the local dynamics, we propose that the topological features which play major roles in shaping the interaction between two nodes are: (i) the strength of the links, (ii) the length of the path between the nodes, and (iii) the routing of information along multiple and redundant paths. A strong direct connection between two nodes is usually a reliable indicator of the strength of their functional interaction. In fact, in general direct links lead to more effective communications since the content of information (e.g., the amplitude of the perturbations) tends to decay over longer paths 28 . However, the flow of information from one node to another does not follow a unique path, but rather spreads along several. Therefore, the total influence that one node exerts over another is accumulated over all possible paths, of all lengths [29][30][31] .
The goal of this paper is that of answering the question of how the structure of a network contributes to sculpt the expected pattern of correlations that emerge when the network hosts a dynamical process. A direct application is that of understanding the extent to which the structural connectivity of the brain determines the largescale correlation structure observed during rest. Accordingly, we have developed a graph measure, the topological similarity T , which estimates the expected correlation between two nodes based on the similarity of the "influences" both receive from the whole network. If two nodes receive the same sets of inputs, then they should be strongly correlated. On the contrary, if they share no common inputs, their correlation tends to be weaker. Hence, this measure can be considered as a generalization of the matching index. In the analysis of graphs the matching index is typically used to quantify the number of common neighbours shared by two nodes. However, the matching index only accounts for the direct neighbours and ignores the rest of the complex fabric of interactions sustained by the entire network. To perform a complete comparison of the inputs accumulated over all paths, while accounting for the natural decay of signal power or information content in any real system, we considered the columns of the communicabiilty measure 32 as they directly correspond to the whole-network input profiles of the nodes. Although the actual decay rate of the signals may differ across real systems, our choice of the communicability guarantees (due to its exponential decay) that the accumulation of perturbations over all possible paths, Eq. (1), converges for all adjacency matrices A. However, any other formalism implementing a decay of the influence between nodes as a function of their graph distance can be used to compute the topological simmilarity T , with the constrain that such formalism must be based on a set of coefficients assuring series convergence for any adjacency matrix.
Recently, Zamora-López et al. 41 found that considering the communicability as the propagator kernel for the diffusion of Gaussian noise along the network allows to analytically estimate the time-average cross-correlation matrix R of the system. We realise that both approaches, the one starting from a dynamical system describing the propagation of perturbations 41 , and the other based uniquely on topological constrains and the assumption of exponential decay of influences, are indeed equivalent. In fact, we find that T ≡ R when the variance ξ i of the Gaussian noise is the same for all nodes in the Gaussian diffusion system. Therefore we can conclude, with a high confidence, that T represents the most likely expected correlation structure of a network only due to its underlying topology.
As argued before, the other crucial ingredients for the collective dynamics on a network are the local dynamics of the nodes and the coupling functions, both of which can be either linear or nonlinear. The topological similarity T thus captures the tendency of the nodes of the network to correlate or to synchronise with each other. The contribution of the local dynamics and of the coupling function is to modulate this original background tendency set by the connection topology and either enhance or disrupt the underlying patterns of correlation.
In order to gain understanding on how each of the three topological features, namely, link weight, path length and path redundancy, precisely affect how influence propagates through the network, we have first investigated the behaviour of C and T between selected nodes in simple networks: chains and cycles of varying length, and pathredundant graphs. Not surprisingly, we have found that the strength of a direct link between two nodes is a major contributor to the intensity of their expected correlation. However, the presence of common inputs or redundant paths between them enhances the intensity of their interaction beyond the baseline determined by the strength of the link. If there is no direct link between two nodes, common inputs and redundant paths can trigger strong correlations between them. However, this influence tends to decay with the length of the paths. Although the precise decay rate of the influence depends on the characteristics of the real system, it is worth noting that real networks with diameter larger than five are rather uncommon. Specially in brain and neural networks which are dense. From a more general perspective, we shall notice that large efforts have been devoted in the literature to investigate which of the common network properties, e.g., degree distribution, clustering coefficient, communities or motifs, determine more prominently the collective dynamics on a network, particularly its capacity to synchronise [48][49][50][51][52][53] . After our observations we can conclude that what truly matters for the collective dynamics is the path structure of the network which determines how far do perturbations reach and how their effect accumulates over redundant and recurrent paths. Thus, we foresee that, to precisely understand how the typical network properties determine the collective dynamics, future work needs to identify how those properties alter the underlying network's path structure.

Limitations and outlook
The results we have here presented come with some limitations which shall be underlined. (A) The topological similarity T we have introduced represents an estimation of the time-averaged cross-correlation. Thus, it only estimates the spatial correlations between brain regions and does not capture temporal correlations. (B) Communicability measure assumes an exponential decay of the influence with pathlength. This is not the only choice possible and other sets of coefficients k l exists which will lead the weighted version of the sum in Eq. (1) converge. The precise decay rate of influence, perturbations or information in real systems will vary according to the system's nature. Thus, in some real applications it might be possible to identify the decay rate with pathlength and identify the "right" set of coefficients k l . As such, T is not restricted by the set of coefficients used in the communicability measure proposed by 32 .(C) We found here more convenient to quantify topological similarity as the cosine similarity of the input profiles, Eq. (4), as we found it more accurate in estimating the similarity between input profiles of small size, as in the case of short chains. However, other measures are possible, e.g., correlation or euclidean distance, without conceptually modifying the meaning of topological similarity. (D) for the analysis on the brain's connectivity and spontaneous correlations, we considered both hemispheres as if they were independent. The reason was to avoid biases due to the unreliability of tractography to identify inter-hemispheric fibers. Still, our comparisons are biased to some extent because the SC-based T and simulations consider both hemispheres as independent while the empirical resting-state measurements reflect the activity of the brain regions, which are, certainly, embedded on the whole network. Finally, we want to highlight the usefulness of both topological similarity T and the correlation matrix R of the Gaussian diffusion system 41 to explore the functional connectivity of synthetic and empirical networks. Because they estimate the expected correlation analytically, they are very fast to compute. For example, they allow to compare the effects of network perturbations, e.g., node or link lessions, without the need to run computationally expensive simulations. The only difference between T and R is that the formalism of R allows also to explore the effects of simulated inputs by increasing of decreasing the variance of the Gaussian noise at selected nodes.

VI. MATERIALS AND METHODS
To obtain the population average structural and functional connectomes twenty-one healthy volunteers (mean age 21.56 years; standard deviation 1.84 years; all males; all right handed) participated in five (5) resting-state and two (2) DTI scanning sessions. All participants signed informed consent to participate. The study was conducted at the School of Psychology, Birmingham and was approved by the University of Birmingham Ethics Committee.

Data acquisition
The scanning sessions were conducted at the Birmingham University Imaging Centre using a 3T Philips Achieva MRI scanner with a 32-channel SENSE head coil. T1-weighted anatomical data (175 slices; 1 × 1 × 1 mm 3 resolution) were collected during the first session only. DTI data were collected in two sessions (23.3 2.5 days apart). The DTI acquisition consisted of 60 isotropically-distributed diffusion weighted directions (b=1500 smm-2; TR=9.5s; TE=78ms; 75 slices; 2 × 2 × 2 mm 3 resolution) plus a single volume without diffusion weighting (b = 0 smm-2, denoted as b0). The DTI sequence was repeated twice during each session, once following the Anterior-to-Posterior phase-encoding direction and once the Posterior-to-Anterior direction, to correct for susceptibility-induced geometric distortions 54 . Resting-state data were collected in five sessions (the first and the last collected in the same scanning session as the DTI data) using whole brain echo-planar imaging (EPI) (TR = 2s; TE = 35ms; 32 slices; 2.5 × 2.5 × 4 mm 3 resolution). Participants were instructed to have their eyes open and maintain fixation to a white dot presented at the centre of the screen.

Whole-brain DTI tractography
We processed the DTI data in FSL version 5.0.8 (FMRIB Software Library, http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) on a Red Hat Linux operating system. We corrected the data for susceptibility distortions, eddy currents and motion artifacts 55 . We subsequently rotated the gradient directions (bvecs) to correct them for motion rotation 42,56,57 . We then used the corrected gradients in the Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques (BedpostX) tool to generate a distribution model in each voxel 58 . We used the default parameters in BedpostX for the diffusion modelling: 2 fibers per voxel, weight of 1 for the secondary fibers, discard of the first 1000 iterations before sampling.
We parcellated the brain into 116 areas using the Automated Anatomical Labeling (AAL) atlas 59 . We followed a 4-step registration procedure to align the AAL atlas from MNI to native space: (a) align the non-weighted diffusion volume (b0) of each session to their midspace and create a midspace-template, (b) align the midspacetemplate to the anatomical (T1) scan, (c) align the T1 to the MNI template of FSL, and (d) invert and combine all the transformation matrices of the previous steps to obtain the MNI-to-native registration. The results of each step were visually inspected to ensure that the alignment was successful.
Step (a) controls for potential bias towards the first session, when the T1 was acquired (similar methodology to Ref 60. The first two registrations are 6dof linear transformations (rigid-body) since we aligned images of the same subject, while the third is 12-dof nonlinear to warp the participants brain around the MNI template. The final matrix of step (d) was applied to the AAL atlas using nearest-neighbour interpolation to preserve the labels of the areas.
We calculated the number of probabilistic streamlines starting from each AAL area and reaching any other AAL area by feeding the BedpostX model to the Probabilistic Tracking algorithm (ProbtrackX) 61 . The parameters we used in ProbtrackX are: 5000 samples per voxel, 2000 steps per sample until conversion, 0.5mm step length, 0.2 curvature threshold, 0.01 volume fraction threshold and loopcheck enabled to prevent streamlines from forming loops. We normalised the number of streamlines by the size of the seed area and thresholded streamlines lower than 1% of maximum (i.e. setting them to zero). We subsequently computed the undirected structural connectivity matrix by averaging the normalised streamlines from area i to area j and from area j to area i. The results for each subject (in each DTI session) were organised into 42 weighted adjacency matrices A of size 116 × 116.

Population-average structural connectome
To estimate the population average structural connectivity (SC) we pooled the 42 SC matrices together (2 per subject) and considered only the reduced parcellation into 90 brain areas (45 per hemisphere) by excluding the 26 regions of the cerebellum and the vermis. The 42 SC matrices contained a variable number L of undirected links ranging from L = 895 for the sparsest case (density ρ = 0.22) to L = 1279 for the densest (ρ = 0.32). We noticed that the simple average of the matrices into a single SC matrix by averaging the 42 values each link takes along the pool leads to an average connectome with strongly biased network properties. For example, this plain average SC matrix contained L = 1967 links, which is almost twice the number of links as in the individual matrices. In order to avoid this problem we have devised a method which automatically removes outlier links before performing the average. For each link (i, j) we have initially a set of 42 weigths {w s ij } where s = 1, 2, . . . , 42. The method searches for outlier weights (data-points falling out of 1.5 times the interquartile range) and removes them from the data pool. The search is iteratively repeated until no further outliers are detected and then the population-average SC weight for the link (i, j) is calculated as the average weight of the surviving values. In practice, the method converges very rapidly and it rarely performs more than 2 iterations per link. This method allows to clean the data without having to set an arbitrary hard threshold 62 for the minimally accepted prevalence of the link. Full details of the method are currently in preparation and will be presented somewhere else. The resulting population-average SC matrix out of our iterative pruning method contains L = 1189 links (ρ = 0.30), which lies within the range of connectivity for the individual 42 matrices. For the simulations we treated the left and the right hemispheres independently, as two matrices of N = 45 ROIs because tractography is known to largely miss interhemispheric connections.

Resting-state time-courses and functional connectome
We pre-processed the EPI resting-state data in FSL version 5.0.8 (FMRIB Software Library, http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) on a Red Hat Linux operating system using MELODIC (Multivariate Exploratory Linear Optimized Decomposition into Independent Components). We corrected the data for motion and slice scan timing, removed the non-brain tissue, applied 5mm FWHM spatial smoothing and removed spike motion artifacts using WaveletDespike 63 . We subsequently applied high-pass temporal filtering and then extracted the average time-course from each AAL area. To estimate the population-average functional connectivity (FC) matrix we considered again only N = 90 brain areas excluding the cerebellum and the vermis. We concatenated the 105 sequences of resting-state signals (21 subjects, 5 sessions per subject) into a single long multivariate time-series and computed the Pearson correlation (z-Fisher corrected) for every pair of signals. The opposite procedure, to compute an FC matrix per session and averaging over the 105 FC matrices leads to almost identical results.

Hopf normal model
Within this model, the temporal evolution of the activity z of node j is given in the complex domain as: z j = ρ j e iθj = x j + iy j Where ω is the node's intrinsic frequency of oscillation, α is the local bifurcation parameter (local because the model allows the possibility to assign a different value of α for each node in the network) and η is additive Gaussian noise with standard deviation σ. This system above shows a supercritical bifurcation at α = 0. Specifically, if α j is smaller than 0, then the local dynamic has a stable fixed point at z j = 0, while for α j values larger than 0 there exists a stable limit-cycle oscillation of frequency f = ω/2π . Whole-brain dynamics are described by the following coupled equations: Where C ij is the anatomical connectivity between nodes i and j, g is the global coupling factor and the standard deviation of gaussian noise is σ = 0.02. In this model the simulated activity corresponds to the BOLD signal of each node. The intrinsic frequency of each node was estimated as the peak frequency in the associated narrowband (i.e., 0.04 -0.07 Hz 64 ) of the empirical BOLD signals of each brain region. We simulated, for ech of the two hemispherese (45 ROIs each), 330000 points using Euler's method for integration (dt = 0.001). The connectivity between all the regions of interest was defined using the empirical structural connectivity matrix (SC), and obtained timeseries were then used to compute the simulated correlation matrix (Simulated FC).

ACKNOWLEDGMENTS
We would like to thank Rui Wang and Caroline di Bernardi Luft for their help in collecting the data used in this study.