Resolving inter-regional communication capacity in the human connectome

Abstract Applications of graph theory to the connectome have inspired several models of how neural signaling unfolds atop its structure. Analytic measures derived from these communication models have mainly been used to extract global characteristics of brain networks, obscuring potentially informative inter-regional relationships. Here we develop a simple standardization method to investigate polysynaptic communication pathways between pairs of cortical regions. This procedure allows us to determine which pairs of nodes are topologically closer and which are further than expected on the basis of their degree. We find that communication pathways delineate canonical functional systems. Relating nodal communication capacity to meta-analytic probabilistic patterns of functional specialization, we also show that areas that are most closely integrated within the network are associated with higher order cognitive functions. We find that these regions’ proclivity towards functional integration could naturally arise from the brain’s anatomical configuration through evenly distributed connections among multiple specialized communities. Throughout, we consider two increasingly constrained null models to disentangle the effects of the network’s topology from those passively endowed by spatial embedding. Altogether, the present findings uncover relationships between polysynaptic communication pathways and the brain’s functional organization across multiple topological levels of analysis and demonstrate that network integration facilitates cognitive integration.

Here we develop a simple method to study the capacity for pairs of brain regions (dyads) to communicate with each other. We deconstruct the conventional global approach and estimate communication capacity without averaging over pairs of regions. In addition, we introduce a procedure to standardize the communication capacity between pairs of regions by their communication capacity in a population of rewired null networks, allowing us to identify pairs of regions with greater or less than expected communication capacity. An important advantage of this method is that it can be used in combination with any null network model. Increasingly constrained surrogates constitute increasingly conservative benchmarks that more closely resemble the empirical network under study (Váša & Mišić, 2022). Using them in parallel can simultaneously allow us to control for a covariate or disentangle its effect from those of other unconstrained features in a more liberal null model (Váša & Mišić, 2022). Degree-preserving null models are the most commonly used for network statistic normalization (Humphries, Gurney, & Prescott, 2006;Maslov & Sneppen, 2002;Newman & Girvan, 2004;van den Heuvel & Sporns, 2011;van Wijk, Stam, & Daffertshofer, 2010). They allow us to mitigate the effect of this simple, but influential, graph feature to assess to what extent a network attribute is unexpected in contrast to a random graph that only preserves the empirical network's degree sequence. In doing so, they highlight the effect of more subtle topological properties on a given statistic. Here, the main null model under consideration further constrains the empirical network's weighted degree sequence using a simulated annealing procedure (Cimini et al., 2019;Khachaturyan, Semenovskaya, & Vainstein, 1979;Kirkpatrick, Gelatt, & Vecchi, 1983;Mišić et al., 2015). In parallel, we additionally consider the effect of the human connectome's geometry with a null model that approximately preserves the edge length distribution of the empirical network, in addition to its degree sequence, allowing us to strictly attribute the remaining effects to topology (R. F. . We initially focus on the topological shortest path (hereafter referred to as a "path"), because (a) it is a simple and fundamental method to infer communication, and (b) in many classes of networks, including brain networks, alternative communication mechanisms nevertheless take advantage of shortest paths without any knowledge of the global topology, including diffusion (Goñi et al., 2013Mišić et al., 2015) and navigation (Seguin, Razi, & Zalesky, 2019;Seguin, Mansour, Sporns, Zalesky, & Calamante, 2022;Seguin et al., 2018Seguin et al., , 2020Vézquez-Rodríguez et al., 2020). We then investigate inter-regional communication capacity, mapping it onto large-scale cognitive systems and patterns of functional specialization. Finally, we also consider the relationship between spatial proximity/geometric embedding and communication capacity, as well as alternative communication mechanisms.

RESULTS
The results are organized as follows. First, we develop a method to standardize communication capacity between pairs of brain regions. We then relate inter-regional communication to the brain's spatial embedding, canonical functional systems, and patterns of functional specialization. All analyses were conducted in a sample of N = 69 healthy participants (source: Dyad: Pair of brain regions. Null models: Models used to produce random surrogate brain networks under specific constraints to benchmark the statistical unexpectedness of an observed brain feature. ▪ Structural connectivity. Structural connectivity was reconstructed from individual participants' diffusion spectrum imaging data using deterministic streamline tractography. A distance-dependent consensus-based thresholding procedure was then used to assemble a group-representative weighted structural connectivity matrix of streamline density (R. F. Betzel, Griffa, Hagmann, & Mišić, 2019;Mišić et al., 2015Mišić et al., , 2018. ▪ Functional connectivity. Functional connectivity was estimated from the same individuals' resting-state functional MRI (rs-fMRI) data using pairwise Pearson correlations among regional time courses. Fisher's r-to-z transformation was applied to individual functional connectivity matrices. A group-average functional connectivity matrix was then computed as the mean across individuals, which was back-transformed to correlation values.
The sample was randomly divided into Discovery (n = 34) and Validation (n = 35) subsets. Analyses were conducted in a high-resolution (1,000 nodes) and low-resolution (219 nodes) parcellation using the Cammoun atlas , a subdivision of the Desikan-Killiany anatomical atlas (Desikan et al., 2006). See Sensitivity Analyses for details.

Benchmarking Dyadic Communication Capacity
To quantify polysynaptic communication capacity between pairs of brain regions, we first compute the topological weighted shortest path lengths on the unthresholded structural connectome ( Figure 1) (Dijkstra, 1959). Shorter weighted path length between a pair of regions Structural connectivity: Map of the anatomical connections among neural elements. Figure 1. Standardization procedure. A population of null structural connectivity matrices that preserve the size, density, and weighted degree sequence of the empirical group-consensus network was generated by randomly rewiring pairs of edges. Weighted shortest path lengths were then computed between every pair of brain regions for the empirical structural brain network and each rewired null. Finally, the path lengths of the empirical network were standardized element-wise against the null population of path lengths from the rewired networks. Lower standardized shortest path length indicates greater communication capacity.
indicates greater communication capacity (Avena-Koenigsberger et al., 2017Rubinov & Sporns, 2010). We simultaneously construct a population of rewired networks that preserve the density and weighted degree sequence of the empirical network (Mišić et al., 2015;Váša & Mišić, 2022). We then compute the path lengths for each rewired network, indicating the communication capacity between pairs of brain regions under the null hypothesis that inter-regional relationships depend only on weighted degree and density ( Figure 1). Finally, we standardize element-wise the empirical path lengths against the population of path lengths in the rewired null networks. The resulting standardized shortest path length matrix quantifies in terms of z-scores how unexpectedly short (<0) or unexpectedly long (>0) communication pathways are between any given pair of brain regions. Figure 2A shows a scatterplot between empirical (abscissa) and rewired (ordinate) path lengths, where each point represents a pair of regions. As expected, the majority of points fall below the identity line (87.15%), suggesting that most path lengths in rewired networks are shorter than in the empirical structural brain network (Watts & Strogatz, 1998). This is in line Figure 2. Benchmarking dyadic communication capacity. (A) Scatterplot between empirical (abscissa) and rewired (ordinate) shortest path lengths, where each point represents a pair of brain regions. Marginal distribution histograms are shown on the top and right axes. Points that appear below the identity line correspond to paths with a shorter length in the rewired networks than in the empirical network, and vice versa for points above the identity line. (B) Distribution of standardized shortest path lengths (z-scores) for all pairs of brain regions. Values less than 0 indicate greater-than-expected communication capacity, and values greater than 0 indicate lower-than-expected communication capacity. (C) Spatial distribution of the top 1% unexpectedly short path lengths. (D) Using the empirical and the standardized path length matrices, closeness centrality (inverse mean path length to the rest of the network) was computed and rank-transformed for every brain region. The two resulting brain maps were then subtracted, resulting in a brain map of the region-wise differences between closeness centrality ranks in the empirical and the standardized networks. Red regions are more integrated in the empirical network, and blue regions are more integrated in the standardized network.
with numerous global accounts of the shortest path length of random networks and their comparison with characteristic path lengths of empirical brain networks (Albert & Barabási, 2002;Bassett & Bullmore, 2006;Hilgetag & Kaiser, 2004;Sporns & Zwi, 2004;Watts & Strogatz, 1998). Interestingly, a number of points reside above identity (12.85%), suggesting that these region pairs enjoy greater-than-expected capacity for communication. Figure 2B further demonstrates this result, showing the distribution of standardized path lengths for all pairs of regions. Negative values indicate dyads with greater-than-expected communication capacity, and positive values indicate dyads with lower-than-expected communication capacity.
To get a sense of how the centrality or "hub-ness" of each brain region changes when path lengths are standardized, we compute the closeness centrality (inverse mean path length to the rest of the network) of each brain region using the empirical and the standardized path length matrices. Figure 2D shows the difference between rank-transformed closeness computed using empirical and standardized path lengths. The figure suggests that the inferred importance of a brain region changes considerably when the procedure is applied. Namely, red regions (e.g., cingulate cortex) are more central in the empirical shortest path length network, and blue regions (e.g., orbitofrontal cortex) are more central in the standardized shortest path length network.

Communication Pathways Delineate Functional Systems
We next consider how communication paths can be contextualized with respect to canonical features of brain networks, including spatial embedding, structure-function coupling, and macroscale intrinsic network organization. Figure 3A shows the relationship between standardized path length and pairwise inter-regional physical distance (left) and pairwise interregional functional connectivity (right). There is a positive association between physical distance and standardized path length, consistent with the notion that areas that are physically further apart have lower communication capacity (Roberts et al., 2016;Seguin et al., 2018;Stiso & Bassett, 2018). There is also a negative association between standardized path length and functional connectivity, consistent with the notion that pairs of areas that are topologically closer have more coherent time courses Honey et al., 2009). Collectively, these results show that standardized path length recapitulates well-known and expected relationships between the topology, geometry, and functional connectivity of the brain (Suárez, Markello, Betzel, & Misic, 2020).
How are communication pathways organized among the canonical macroscale intrinsic networks? Resting-state functional connectivity networks are communities of functionally related areas with coherent time courses that are thought to be putative building blocks of higher cognition (Bellec, Rosa-Neto, Lyttelton, Benali, & Evans, 2010;Damoiseaux et al., 2006;Power et al., 2011;Yeo et al., 2011), but how these networks map onto the underlying communication pathways is not completely understood (Avena-Koenigsberger et al., 2018;Suárez et al., 2020). To address this question, we first stratify brain regions according to their membership in the intrinsic networks derived by Yeo, Krienen and colleagues  ( Figure 3B). Figure 3C (left) shows the mean standardized path length within each intrinsic network (diagonal) and between all pairs of intrinsic networks (off-diagonal). We generally observe shorter path lengths within networks compared with between networks; Figure 3C (right) confirms this intuition, showing that the mean within-network path length is significantly shorter than the mean between-network path length (p spin < 0.001).
Next, we quantify and compare the internal communication capacities of pairs of intrinsic networks by computing the difference between their respective within-network mean Intrinsic networks: Networks of functionally related brain regions with synchronous spontaneous activity. standardized path length ( Figure 3C, middle). We find that the frontoparietal network has consistently greater internal communication capacity compared with other networks, while the somatomotor network has consistently lower internal communication capacity compared with other networks. Interestingly, Figure 3C (left) also indicates that communication pathways internal to the frontoparietal network are the only ones to exhibit a greater-than- path length (y) between two brain regions grows as a function of the Euclidean distance (x) between them (left). The black line corresponds to the fitted exponential y = −9.31e −0.03x + 4.33. Functional connectivity (y) between two brain regions decays as a function of the standardized shortest path length (x) between them (right). The black line corresponds to the fitted exponential y = 0.12e −0.18x . (B) Standardized shortest path length matrix with brain regions ordered based on their affiliations to the Yeo intrinsic networks. (C) Left: Heatmap of the mean standardized path lengths across node pairs belonging to the same intrinsic network (diagonal) and to different intrinsic networks (off-diagonal). A blue square identifies a negative mean standardized path length, indicative of shorter-than-expected communication pathways with greater-than-expected communication capacity. Middle: Heatmap of the pairwise differences of the means among Yeo intrinsic networks, calculated as the mean value of the network on the y-axis minus the mean value of the network on the x-axis, with the mean value corresponding to the mean standardized path length across node pairs belonging to the same network (diagonal elements of the left heatmap). A purple square indicates significant difference of the means based on network label permutation using spatial autocorrelation-preserving null models (Bonferroni corrected, α = 0.05), whereas "n.s." denotes not significant differences. The frontoparietal network displays a consistently shorter mean standardized path length (i.e., higher internal communication capacity) compared with other networks, whereas the somatomotor network exhibits a consistently greater mean standardized path length (i.e., lower internal communication capacity) in comparison to other networks. Right: The mean within-network standardized path length is significantly shorter than the mean between-network standardized path length (p spin < 0.001). (D) Same as panel C but for shortest path lengths standardized using a geometry-preserving null model. Intrinsic networks: vis = visual, sm = somatomotor, da = dorsal attention, va = ventral attention, lim = limbic, fp = frontoparietal, dm = default mode. expected communication capacity, characterized by a negative mean standardized path length.
Moreover, we reproduce this analysis using a geometry-preserving null model ( Figure 3D). Once again, we find significantly shorter path lengths within Yeo networks (p spin < 0.001; Figure 3D, right) and identify the frontoparietal network and the somatomotor networks as consistently exhibiting the highest and lowest internal communication capacity, respectively ( Figure 3D, middle). However, most pairwise comparisons between intrinsic networks are now not significant as assessed by spatial autocorrelation-preserving nulls (Bonferroni corrected, α = 0.05). Furthermore, the frontoparietal network does not present a negative mean standardized path length anymore, indicating that its greater-than-expected communication capacity can be partly attributed to its spatial embedding ( Figure 3D, left).

Communication Capacity and Functional Specialization
Given that communication capacity is regionally heterogeneous and maps onto intrinsic networks, we ask whether regional communication capacity is related to functional specialization. Figure 4A shows the mean standardized path length from each region to the rest of the network, with red indicating greater integration with the network and yellow indicating lower integration.
We statistically compare this map with a library of meta-analytic task-based fMRI activation maps from the Neurosynth repository (Hansen et al., 2021;Yarkoni et al., 2011). Each of the Neurosynth brain maps consists of region-wise measures of the probability that a particular term is reported in a study if an activation was observed in a given region. In this analysis, we only consider the intersection of terms from the Neurosynth database and the Cognitive atlas , comprising a total of 123 cognitive and behavioral terms. Figure 4B shows statistically significant spatial correlations between the node-wise standardized path length map and each of the Neurosynth term maps, as assessed using spatial autocorrelation-preserving null models (Alexander-Bloch et al., 2018;; p spin < 0.05 in gray; Bonferroni corrected, α = 0.05 in green). We find anticorrelations with higher order cognitive terms (e.g., "monitoring," "strategy"). This suggests that areas that communicate closely with many other areas in the connectome are associated with higher order cognitive function. In other words, cognitive integration appears to be supported by network integration.
To better understand why lower standardized path length is associated with higher order cognitive function, we compare the regional map of standardized path length with maps of weighted degree (sum of edge weights incident on a node; strength), betweenness (proportion of shortest paths that traverse a node), and participation (distribution of node links among functional network communities) ( Figure 4C). All three topological features were computed on the empirical weighted structural network. As expected due to the standardization procedure, there is no significant correlation with weighted degree (r s = −0.19, p spin = 0.10). However, standardized path length is significantly negatively correlated with betweenness (r s = −0.45, p spin < 0.001) and participation (r s = −0.37, p spin < 0.001), suggesting that regions that are closely integrated into the connectome can better sample information from multiple specialized communities.
Furthermore, we replicate this analysis using a geometry-preserving null model ( Figure 4D-E). Once again, we find significant anticorrelations between node-wise standardized path lengths and higher order cognitive terms (e.g., "monitoring", "expectancy"). This suggests that the observed association between cognitive and network integration is not passively endowed by the brain's physical embedding, but rather driven by its topological organization.

Extending to Multiple Communication Models
So far, we have only considered path length as a proxy for communication. However, there exist numerous other models of communication in the connectome (Avena-Koenigsberger et al., 2018). Here we extend the framework developed in the previous sections to additional measures of communication proximity that have been proposed for brain networks, including search information Rosvall, Grönlund, Minnhagen, & Sneppen, 2005), path transitivity , communicability (Crofts & Higham, 2009;Estrada & Hatano, Figure 4. Communication capacity and functional specialization. (A) Brain map of mean standardized path length from each node to the rest of the network, with red denoting a greater integration of the node within the network and yellow denoting a lower integration. (B) Statistically significant spatial correlations based on spatial autocorrelation-preserving nulls (p spin < 0.05 in gray; Bonferroni corrected, α = 0.05 in green) between node-wise mean standardized path lengths and meta-analytic probabilistic functional activation maps from the Neurosynth platform, associated with 123 terms overlapping the Neurosynth database and the Cognitive atlas. An important set of anticorrelations suggests that highly integrated nodes are associated with higher order cognitive functions. (C) Relationships between node-wise mean standardized path length and topological features of the empirical weighted structural network. As expected due to the standardization procedure, node-wise mean standardized path length is not significantly correlated with weighted degree (r s = −0.19, p spin = 0.10; left), while it is significantly negatively correlated with betweenness (r s = −0.45, p spin < 0.001; middle) and participation (r s = −0.37, p spin < 0.001; right), suggesting that regions that are more topologically integrated also have a more diverse connection profile among functionally specialized intrinsic networks. (D-E) Same as panel C but for shortest path lengths standardized using a geometry-preserving null model. 2008), and mean first-passage time (Goñi et al., 2013;Noh & Rieger, 2004). As before, we first standardize each communication matrix using a population of rewired networks ( Figure 5A). We then extract all dyadic (i, j ) elements from each communication matrix and assemble them into the columns of a dyads × communication models matrix.
Applying principal component analysis (PCA) to this matrix identifies a single dominant component that accounts for 62.79% of the dyad-level variance in communication. The resulting PC1 scores are then used to constitute an "aggregated" communication matrix that shows the capacity for communication among all pairs of brain areas across multiple communication models. The greatest contribution to the aggregated communication measure is from search information, path length, and mean first-passage time, with modest contribution from communicability. Note that, unlike the original communication measures, this aggregated communication measure is not necessarily transitive.
Overall, we observe similar results using the aggregated multicommunication measure ( Figure 5). Namely, we find a positive relationship between standardized communication distance and physical distance (r s = 0.53, p ≈ 0; Figure 5B), a negative relationship between standardized communication distance and functional connectivity (r s = −0.21, p ≈ 0; Figure 5B), and significantly shorter communication distance within canonical intrinsic networks than between networks ( p spin < 0.001; Figure 5B). From left to right: Positive Spearman correlation between PC1 score and Euclidean distance (r s = 0.53, p ≈ 0), negative Spearman correlation between functional connectivity and PC1 score (r s = −0.21, p ≈ 0), and significantly lower within-network than between-network PC1 score (p spin < 0.001). Figure S1A in the Supporting Information shows the same procedure applied at the node level. As before, we find consistent results, showing that areas that are closer in communication distance to other areas in the connectome tend to be associated with higher order cognitive function ( Figure S1C) and greater participation in network communities (r s = −0.40, p spin < 0.001; Figure S1D; right).

Disentangling the Contributions of Topology and Geometry
Physical proximity is an important predictor of connection probability and connection weight (R. F. Roberts et al., 2016), and therefore by extension, communication capacity between areas. We therefore seek to disentangle the contribution of geometry from the contribution of topology to the results reported so far. Initially, we utilized degreepreserving rewiring to construct populations of null networks and standardize path length (Maslov & Sneppen, 2002;Mišić et al., 2015). To additionally control for the role of geometry, we repeat the experiments using a more conservative null model that approximately preserves the edge length distribution and the edge length-to-weight mapping, in addition to degree sequence (R. F. . Supporting Information Figure S2 shows the results when applying the geometry-preserving null model. As for the degree-preserving null model, Figure S2A shows that the majority of path lengths are shorter in the rewired networks than in the empirical network (86.04%). However, the majority of points (region pairs) now lie closer to the identity line, suggesting that the physical embedding of the structural brain network contributes in pulling its regions apart topologically. Furthermore, controlling for geometry in the standardization procedure accentuates the integration and segregation of certain regions, providing a clearer, more homogeneous picture of changes in closeness centrality ( Figure S2B). The orbitofrontal cortex becomes more integrated and the cingulate cortex and precuneus become more segregated when strictly considering the empirical network's topology. Interestingly, the paracentral lobule and the motor cortex now move in directions opposite to those observed when using the null model that did not preserve the connectome's geometry. Indeed, the geometry-preserving standardization procedure moved the paracentral lobule topologically closer to the rest of the network and the motor cortex further. This suggests that the physical distance of these regions to the rest of the brain interacts with its topology to centralize the motor cortex and segregate the paracentral lobule.
As expected, controlling for edge length considerably attenuates the relationship between standardized path length and Euclidean distance (Supporting Information Figure S2C, left). Importantly, most of the other results are largely preserved, including shorter standardized path lengths within compared with between networks (p spin < 0.001; Figure 3D, right) and a relationship between node-wise communication and integrative function ( Figure 4E). Collectively, this control experiment suggests that most of the results reported previously-except the relationship with Euclidean distance-are mainly driven by the topological organization of the connectome rather than spatial embedding and geometric relationships. An important exception is the higher-than-expected internal communication capacity of the frontoparietal network, which seems to have been partly driven by the brain's geometry.

Sensitivity Analyses
We test the replicability of the findings in the Validation sample (Supporting Information Figure  S3) and seek to assess the sensitivity of the results to a variety of processing choices. We repeat Rewiring: Randomly swapping pairs of connections under specific constraints.
We find consistent results across all sensitivity analyses. This includes significantly shorter communication pathways between cortical regions belonging to the same intrinsic functional network than between regions belonging to different intrinsic networks, as well as significant relationships between a region's topological integration and its association with higher order executive functions (Supporting Information Figures S3, S4, S5). The frontoparietal network is consistently identified as exhibiting the highest internal communication capacity. Communication pathways internal to the frontoparietal network also have a negative mean standardized path length, indicative of a greater-than-expected communication capacity (Figures S3,  S4, S5).
Next, we seek to test the extent to which the results are influenced by the inclusion of direct monosynaptic pathways (i.e., paths between directly anatomically connected nodes; Supporting Information Figure S6). We therefore repeat the analyses in strictly polysynaptic communication pathways (i.e., paths between pairs of nodes that are separated by two or more anatomical connections; Supporting Information Figures S6C, S7, S8, S9). Again, we find consistent results with the notable exception of the communication capacity of the frontoparietal network ( Figure S7; original result in Figure 3C). This indicates that the greater-than-expected communication capacity of this network is partly driven by monosynaptic connections.
Recently, there has been a growing interest to relate patterns of brain communication with the higher order connectivity of brain networks (Battiston et al., 2021;Crofts, Forrester, & O'Dea, 2016;Griffa et al., 2022;Sizemore et al., 2018). In alignment with this new line of inquiry, and to demonstrate the flexibility of our standardization procedure, we apply it to a recently introduced higher order multimodal communication measure: structure-function bandwidth (Parsons et al., 2022). Based on a multilayer framework in which structural and functional connectivity are considered simultaneously, bandwidth between two brain regions in the functional connectivity layer is defined according to the minimum edge weight of a path connecting these nodes in the structural connectivity layer. The maximum bandwidth is then selected across all paths of a given length. Here, we consider triangles composed of two-hop structural paths closed by a functional edge. We weigh structure-function bandwidth by the functional edge weights to provide a generalized communication measure taking into account all functional connections enclosing triangles. This yields a structure-function bandwidth matrix that we concurrently standardize using the same procedure previously described, with structure-function bandwidth being recomputed in a population of rewired structural networks while the functional network is left intact. We then consider how our generalized structurefunction bandwidth measure is organized among intrinsic functional networks. We observe distinct patterns of intra-and inter-network structure-function bandwidth when comparing the empirical and the standardized measures (Supporting Information Figure S10). This analysis demonstrates how the standardization procedure can be readily extended to accommodate questions about higher order brain network architecture.

DISCUSSION
In the present report, we introduce a simple method to standardize communication path lengths in brain networks. These results showcase how dyadic relationships can be resolved and studied while accounting for more basic topological and geometric features of the network. Building on hierarchically constrained null models, this rigorous standardization procedure enables specific quantification and localization of the degree of unexpectedness of communication measures. In contrast to classical approaches, considering communication capacity at a finer granularity allows us to revisit previous investigations of brain hubs, recapitulate well-known geometric and functional attributes of inter-regional communication, and uncover new relationships between the human connectome's topological and functional architectures across multiple topological levels of analysis in a nuanced and principled way.
Classical studies focused on global path length or efficiency of brain networks. These reports found evidence of near-minimal path length characteristic of small-world architecture across multiple species and reconstruction methods (Bassett & Bullmore, 2017;Hilgetag & Kaiser, 2004;Kaiser & Hilgetag, 2006;Sporns & Zwi, 2004;Watts & Strogatz, 1998). In addition, empirical studies found that low characteristic path length or high global efficiency is associated with greater cognitive performance (Li et al., 2009;van den Heuvel, Stam, et al., 2009) and is concomitant with healthy neurodevelopment (Baum et al., 2017;Fan et al., 2011;Hagmann et al., 2010;Khundrakpam, Lewis, Zhao, Chouinard-Decorte, & Evans, 2016;Khundrakpam et al., 2013). Altogether, these findings speak to the behavioral and biological relevance of global accounts of communication capacity.
However, global communication measures such as characteristic path length or global efficiency are effectively summary statistics of a myriad of complex inter-regional communication relationships. In the present study, we focus specifically on dyadic communication while statistically controlling for fundamental topological and geometric features (i.e., degree and spatial position). We find that most communication pathways between brain regions are longer than expected on the basis of their degree and/or spatial position. Despite the standardization procedure, we still recapitulate fundamental features of inter-regional communication, such as positive relationships with spatial proximity (Roberts et al., 2016;Seguin et al., 2018;Stiso & Bassett, 2018) and negative relationships with resting-state functional connectivity Honey et al., 2009;Seguin et al., 2020;Vázquez-Rodríguez et al., 2019;Zamani Esfahlani et al., 2022).
We note that the procedure used here is similar to what is typically done when considering global communication statistics, such as characteristic path length. Namely, characteristic path length is often normalized by the mean characteristic path length across a population of null networks, such as when estimating the small-world coefficient (Humphries et al., 2006;Watts & Strogatz, 1998). We build on this general approach by finely resolving dyadic communication relationships. In addition, by standardizing each dyadic path length as opposed to normalizing by the mean, we implicitly take into account variance across the null population (Bassett et al., 2008).
The standardization procedure alters the centrality rankings of brain regions, suggesting that taking the constraints of degree into account can lead to different inferences about the functional importance of brain regions. For instance, we find that the orbitofrontal cortex is more integrated in the standardized network than the empirical network, and that parts of the medial prefrontal cortex, cingulate cortex, and precuneus are less integrated in the standardized network. These results run counter to numerous classical investigations of brain hubs that did not explicitly control for degree when estimating different centrality measures based on path length, such as betweenness and closeness (Hagmann et al., 2008;van den Heuvel & Sporns, 2011, 2013b. A large body of neuroimaging studies has highlighted the existence of a number of macroscale communities of functionally related brain regions with correlated resting-state fMRI signals (Damoiseaux et al., 2006;Power et al., 2011;Uddin, Yeo, & Spreng, 2019;Yeo et al., 2011). Previous efforts were made to relate these patterns of synchronized spontaneous activity to the underlying anatomical scaffolding of the brain (van den Heuvel & Hulshoff Pol, 2010). Specifically, some studies have investigated the structural underpinnings of the functional links within and between resting-state networks by identifying specific white matter tracts that could mediate these relationships (Greicius, Supekar, Menon, & Dougherty, 2009;van den Heuvel, Mandl, Kahn, & Hulshoff Pol, 2009). More recently, a growing interest in the functional predictive utility of communication measures based on structural connectivity has led to investigations into the relationships between functional brain networks and underlying patterns of polysynaptic communication. It has been shown that the intrinsic functional hierarchy of the brain guides communication trajectories and allows for signal diversification (Vézquez-Rodríguez et al., 2020). Furthermore, it was found that the modular boundaries of resting-state functional connectivity were approximated by modules of polysynaptic communication distance . In line with these studies, we map communication pathways among canonical intrinsic functional networks . We find that standardized path lengths are significantly shorter within than between intrinsic networks, suggesting that the topological organization of the human connectome contributes in giving rise to macroscale intrinsic patterns of functional interactions.
Moreover, by organizing communication pathways within individual intrinsic networks and between pairs of networks, we identify the frontoparietal network as exhibiting the highest internal communication capacity. Previous findings had associated a greater structural global efficiency of a frontoparietal network (mean closeness centrality across nodes of the network to all brain regions) with a higher working memory capacity (Pineda-Pardo, Martínez, Román, & Colom, 2016). The present result suggests that, in addition to global integration of frontoparietal nodes, the higher order executive control functions that have been widely attributed to frontoparietal networks (Dosenbach et al., 2007;Fedorenko, Duncan, & Kanwisher, 2013;Laird et al., 2011;Niendam et al., 2012;Power & Petersen, 2013;Vincent, Kahn, Snyder, Raichle, & Buckner, 2008;Yeo et al., 2015) might also benefit from an unexpectedly high level of internal communication capacity.
What are the functional consequences of lower or greater communication capacity? Comparison with meta-analytic maps of functional specialization suggests that regions that are topologically closer to others tend to be associated with higher order cognitive functions such as monitoring and strategy. In other words, we find that greater network integration is associated with cognitive integration (Cole, Yarkoni, Repovš, Anticevic, & Braver, 2012;Mišić et al., 2015;Sporns, 2013;Worrell, Rumschlag, Betzel, Sporns, & Mišić, 2017). This is consistent with numerous theories that posit that patterns of regional specialization arise from connectivity profiles (Bressler & Menon, 2010;Mars et al., 2016Mars et al., , 2018McIntosh, 2000;Passingham, Stephan, & Kötter, 2002) and topological embedding of brain regions (Petersen & Sporns, 2015;Sporns, 2011Sporns, , 2013Wig, 2017). Moreover, we find that nodal topological integration is positively associated with the number of shortest paths traversing it and the diversity of a region's connections among intrinsic functional communities. This is in line with previous accounts of hub characteristics (Hagmann et al., 2008;van den Heuvel, Kahn, Goñi, & Sporns, 2012;van den Heuvel & Sporns, 2011, 2013a, 2013b. Here, we further show these features to be related even when controlling for the effects of degree and strength. Previous reports have also shown positive associations between a node's involvement in complex tasks and the diversity and flexibility of its functional links to the rest of the brain, especially in frontoparietal regions (Bertolero, Yeo, & D'Esposito, 2015;Cole et al., 2013;Power et al., 2011). Altogether, the present results complement these findings, suggesting that brain regions that subserve higher order cognition also benefit from a structural substrate for the diversification and integration of information.
This work is part of a wider trend in the field to infer and quantify the potential for communication among brain regions based on their wiring patterns (Avena-Koenigsberger et al., 2018;Graham, Avena-Koenigsberger, & Mišić, 2020;Graham & Rockmore, 2011;Srivastava et al., 2020;Zhou et al., 2022). Although we mainly focus on shortest paths, multiple alternative communication protocols have been proposed (Avena-Koenigsberger et al., 2017;R. F. Betzel, Faskowitz, Mišić, Sporns, & Seguin, 2022;Crofts & Higham, 2009;Goñi et al., 2014;Mišić et al., 2015;Seguin et al., 2018;Zhou et al., 2022). The present standardization procedure can be readily applied to any dyadic communication measure. Combining additional measures of decentralized communication such as search information and path transitivity, we find results consistent with those derived using path length. As the field moves towards more biologically realistic and validated communication protocols, future studies could adapt this standardization procedure to accommodate emerging measures of communication capacity. As an example, we show that it can be readily adapted to a higher order multimodal communication measure defined in a multiplex brain network combining structural and functional connectivity.
In addition, the present procedure standardizes communication measures using two common types of rewiring null models. Here we focus on disentangling the contribution of the structural brain network's topology from the background effect of spatial embedding. We therefore apply one null model that preserves the (weighted) degree sequence and another that additionally preserves wiring length (R. F. . Interestingly, most effects are preserved when applying the geometric null model, suggesting that they are driven by topology rather than by spatial embedding. More generally, this highlights the fact that any null model-embodying a specific null hypothesis-could be used for the standardization procedure to selectively tease apart features that shape communication in brain networks (Váša & Mišić, 2022;Zamora-López & Brasselet, 2019).
The present findings should be viewed in light of multiple methodological limitations. First, all networks were reconstructed using diffusion-weighted MRI, a method known to be subject to false positives and false negatives (de Reus & van den Heuvel, 2013;Maier-Hein et al., 2017;Thomas et al., 2014). Although we attempt to mitigate this limitation by splitting the sample and repeating all analyses, future work is necessary to improve the quality of connectome reconstruction. Second, networks reconstructed from diffusion-weighted MRI are undirected, limiting the biological plausibility of these networks and our capacity to fit communication models to them. Third, we focus on several well-studied and mathematically fundamental communication protocols, including centralized (e.g., shortest paths) and decentralized (e.g., search information) measures, but this selection is non-exhaustive and alternative communication measures could be considered in future work.
In summary, we present a simple method to resolve communication capacity in brain networks. The method is based on conventional procedures already commonplace in connectomics, but allows researchers to focus on dyadic communication. This procedure is inherently flexible, being able to accommodate emerging communication measures and null models. As the field develops increasingly sophisticated and biologically realistic generative models of inter-regional communication, this procedure will allow greater insight into features that shape signaling patterns in brain networks.

Data Acquisition
A total of n = 69 healthy participants (25 females, age 28.8 ± 8.9 years old) were scanned at the Lausanne University Hospital in a 3-Tesla MRI Scanner (Trio, Siemens Medical, Germany) using a 32-channel head coil . The protocol included (a) a magnetizationprepared rapid acquisition gradient echo (MPRAGE) sequence sensitive to white/gray matter contrast (1 mm in-plane resolution, 1.2 mm slice thickness), (b) a diffusion spectrum imaging (DSI) sequence (128 diffusion-weighted volumes and a single b0 volume, maximum b-value 8,000 s/mm 2 , 2.2 × 2.2 × 3.0 mm voxel size), and (c) a gradient echo-planar imaging (EPI) sequence sensitive to blood oxygen level-dependent (BOLD) contrast (3.3 mm in-plane resolution and slice thickness with a 0.3 mm gap, TR 1,920 ms, resulting in 280 images per participant). The last sequence was used as part of an eyes-open resting-state fMRI (rs-fMRI) scan in which the participants were not overtly engaged in a task. Informed written consent was obtained for all participants in accordance with institutional guidelines and the protocol was approved by the Ethics Committee of Clinical Research of the Faculty of Biology and Medicine, University of Lausanne, Switzerland.

Network Reconstruction
Structural connectomes were reconstructed for individual participants using deterministic streamline tractography and divided according to a gray matter parcellation of 1,000 cortical nodes . The analyses were also repeated at a coarser 219 corticalregion resolution. White matter and gray matter were segmented from the MPRAGE volumes using the FreeSurfer version 5.0.0 open-source package, whereas DSI data preprocessing was implemented with tools from the Connectome Mapper open-source software (Daducci et al., 2012), initiating 32 streamline propagations per diffusion direction for each white matter voxel . Structural connectivity was defined as streamline density between node pairs, that is, the number of streamlines between two regions normalized by the mean length of the streamlines and the mean surface area of the regions (Hagmann et al., 2008). The fMRI data underwent regression of physiological variables, including white matter, cerebrospinal fluid, and motion (estimated via rigid body coregistration). BOLD time series were subsequently subjected to a lowpass temporal Gaussian filter with 1.92 s full width half maximum, and motion "scrubbing" (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) was performed after excluding the first four time points for the time series to stabilize. Functional connectivity was then computed as the Pearson correlation coefficient between the fMRI BOLD time series of each node pair.
The data were randomly split into Discovery (n = 34) and Validation (n = 35) subsets. We then generated group-representative brain networks for each subset to amplify signal-to-noise ratio using functions from the netneurotools open-source package (https://netneurotools .readthedocs.io/en/latest/index.html). A consensus approach that preserves (a) the mean density across participants and (b) the participant-level edge length distributions was adopted for the structural connectomes (R. F. Betzel et al., 2019). First, the cumulative edge length distribution across individual participants' structural connectivity matrices is divided into M bins, M corresponding to the average number of edges across participants. The edge occurring most frequently across participants is then selected within each bin, breaking ties by selecting the higher weighted edge on average. This procedure was applied separately for intra-and interhemispheric edges to ensure that the latter are not underrepresented. The selected edges constitute the distance-dependent group-consensus structural brain network. Finally, the weight of each edge is computed as the mean across participants. The grouprepresentative functional connectivity matrix was defined as the group average following Fisher's r-to-z transformation. The final group-consensus matrix was back-transformed to correlation values.

Communication Models
In this section, we define the analytic communication measures associated with the network communication models considered in the present study and provide their implementation details. All the communication measures, with the exception of path transitivity, were computed using functions from the open-access Python version of the Brain Connectivity Toolbox (https://github.com/aestrivex/bctpy; Rubinov & Sporns, 2010). Path transitivity was implemented in Python based on a MATLAB script openly available in the Brain Connectivity Toolbox.
Shortest paths. Let A denote a weighted adjacency matrix. To identify the sequence of unique edges π u→v = {A ui , …, A jv } spanning the minimum length between nodes u and v (i.e., shortest path), we first defined a monotonic transformation from edge weights, namely streamline density in the present case, to edge lengths L, which can be more intuitively considered as the cost of traversing this edge. We used the negative natural logarithm: L ui = −log(A ui ), mapping greater streamline density to lower cost of signal propagation. Dijkstra's algorithm (Dijkstra, 1959) was used to identify shortest paths, and their length was computed as the sum L uv = L ui + … + L jv of traversed edges' lengths.
Search information. Search information quantifies the amount of information required by a naïve random walker to travel along a specific path in a network Rosvall et al., 2005). More specifically, here we consider shortest paths π u→v , capturing the accessibility of these optimal routes in network topology. This analytic measure is derived from the probability that a random walker starting at u follows Ω u→v = {u, …, v}, the sequence of nodes visited along π u→v to reach v. This probability depends on the strength of the nodes comprised in Ω u→v and can be expressed as follows: Search information can then be defined as S π u→v ð Þ¼− log 2 P π u→v ð Þ ð Þ : Note that this definition assumes that the random walker has no memory of its previous step, which would reduce the information needed to determine the next. Importantly, this measure is not symmetric. That is, S(π u→v ) ≠ S(π v→u ). Search information is contingent on the assignment of the path's source and target nodes.
Path transitivity. Path transitivity is a measure of transitivity based on the shortest path between a pair of nodes. It quantifies the density of one-step detours (triangles) available along the path . Intuitively, it can be considered to be the accessibility of the shortest path or its robustness to a deviation of the signal traversing it. Path transitivity can be computed as the average pairwise matching index of the nodes comprising the path.
The matching index between two nodes is a measure of the similarity of their connectivity profiles, excluding edges incident on both nodes (Hilgetag, Burns, O'Neill, Scannell, & Young, 2000). The matching index M between nodes s and t is defined as where Θ is the Heaviside step function.
Building on this definition, the path transitivity P of the shortest path π u→v can be defined as Communicability. Communicability between two nodes is defined as the length-weighted sum of all walks between them, with longer walks being more penalized (Estrada & Hatano, 2008).
The communicability matrix C of pairwise communicability estimates between all nodes in the network is calculated as the matrix exponential of the adjacency matrix: C = e A . Following Crofts and Higham (2009), in the case of a weighted structural connectivity matrix, we first normalize the adjacency matrix as D 1/2 AD 1/2 , where D is the diagonal weighted degree matrix. The objective of this procedure is to reduce the undue influence of nodes with large strength.
Mean first-passage time. The mean first-passage time from node u to node v is the expected number of hops in a random walk evolved by a random walker starting at node u before arriving for the first time at node v (Noh & Rieger, 2004). Considering nodes of an undirected connected network as the states of an ergodic Markov chain allowed Goñi et al. (2013) to derive the mean first-passage time T between nodes u and v as where w is the steady-state vector of the underlying Markov process and Z is the fundamental matrix calculated as (I − P + W ) −1 . P is the transition probability matrix computed as D −1 A, where D is the diagonal weighted degree matrix. W is a square matrix whose columns correspond to w. In the present study, following Zamani Esfahlani et al. (2022), we standardize as z-scores the columns of the matrix of pairwise mean first-passage time among all nodes in the network to remove nodal bias.
Structure-function bandwidth. Structure-function bandwidth quantifies how well structural connectivity throughput mediates functional connectivity in a multiplex network composed of a structural connectivity layer and a functional connectivity layer. Bandwidth between two nodes in the functional connectivity layer is defined according to the minimum edge weight of a path connecting these nodes in the structural connectivity layer. The maximum bandwidth is then selected across all paths of a given length. Here, we consider triangles composed of two-hop structural paths closed by a functional edge. We weigh structure-function bandwidths by the functional edge weights to provide a generalized communication measure taking into account all node pairs in the functional connectivity layer. The functional connectivityweighted bandwidth B between nodes u and v is computed as follows: where F is the weighted functional connectivity matrix, and S is the binary structural connectivity matrix. The (1 − S uv ) term excludes all cases where a direct structural edge also connects nodes u and v. Every possible intermediary node i in the structural connectivity layer is considered to complete the triangle.

Null Models
Standardization. Network features are interrelated, with many complex network properties depending on basic features (Váša & Mišić, 2022). To mitigate the effect of differences in simple features on the topological relationships under study, we provide a frame of reference to communication measures. First, we built a population of 100 null structural connectivity matrices by randomly rewiring pairs of edges of the empirical group-consensus networks, systematically disrupting their topology while maintaining basic network features. Next, communication measures were computed for every node pair of the empirical structural brain networks and all 100 surrogate graphs. Finally, the empirical communication matrices were demeaned and standardized to unit variance elementwise using the null population. This procedure yielded, in units of standard deviation, an approximate communication measure in relation to what would be expected by chance in a similar network.
Two null models were considered as part of this procedure, resulting in a hierarchy of preserved network features. Both surrogate models guarantee the connectedness of the produced rewired network, that is, no node is disconnected from the rest of the network.
The first method randomly swaps pairs of edges (approximately 10 swaps per edge) while maintaining network size (i.e., number of nodes), density (i.e., proportion of possible edges expressed), and degree sequence (i.e., number of edges incident to each node) (Maslov & Sneppen, 2002). We applied an implementation of this technique openly available in the Python version of the Brain Connectivity Toolbox (https://github.com/aestrivex/bctpy; Rubinov & Sporns, 2010). To extend this procedure to weighted structural connectivity matrices, we then used simulated annealing to further preserve the empirical network's strength sequence (sum of edge weights incident to each node) (Mišić et al., 2015). Simulated annealing is a stochastic search method used to approximate the global optimum of a given function (Kirkpatrick et al., 1983). This is achieved through the Metropolis procedure (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953), controlled by the temperature parameter T. The simulated annealing process is initiated at a high temperature that allows the exploration of costly system configurations, preventing the algorithm from getting stuck in local minima. Throughout the process, the system is slowly cooled down while descending along the optimization landscape, with increasingly limited uphill rearrangements and smaller, fine-tuned changes in the system cost. Here, we minimize the cost function E defined as the sum of squares between the strength sequence vectors of the empirical and the randomized networks: where s i and s 0 i are the strengths of node i in the empirical and the null networks, respectively. To optimize this function, weights were randomly permuted among edges. Rearrangements were only accepted if they lowered the cost of the system or followed the probabilistic Metropolis acceptance criterion: r < exp(−(E 0 − E )/T ), where r ∼ U(0, 1). We used an annealing schedule composed of 100 stages of 10,000 permutations with an initial temperature of 1,000, halved after each stage.
A limitation of this rewiring procedure is that, on average, randomly flipping pairs of edges will result in unrealistically large distances between nodes (Váša & Mišić, 2022). To address this issue, we use a Python implementation of an approach proposed in R. F. Betzel and Bassett (2018), openly available as a function in the netneurotools package (https://netneurotools .readthedocs.io/en/latest/index.html). This method applies the same iterative rewiring procedure as the first, but with additional constraints. Edges are binned according to the Euclidean distance between the centroids of their associated parcel pair. The number of bins was determined heuristically as the square root of the number of edges in the empirical group-consensus network. Swapping is then performed within each bin to approximately preserve the edge length distribution, in addition to exactly reproducing the network's size, density, and degree sequence like the first method. The maximum total number of edge swaps to perform was set to the network size times 20. In maintaining the geometry of the empirical network, this rewiring model provides a more representative surrogate, resulting in a more conservative null model.
Collectively, the two null models used for the standardization of the structural connectivity matrices express a hierarchy of constraints. When considered in parallel, they allow us to distinguish the contribution of the structural brain network's topology from effects passively endowed by its spatial embedding when studying the connectome's architecture (Roberts et al., 2016).
Significance testing. To assess the significance of statistics based on node-level communication measures, we relied on spatial permutation null models . We generated null distributions of statistical estimates derived from permuted brain maps, while preserving the spatial autocorrelation of the original data to respect the assumption of exchangeability.
First, we parcellated the FreeSurfer fsaverage surface according to the Cammoun atlas  using tools from the Connectome Mapper (Daducci et al., 2012). A spherical projection of the fsaverage surface was then used to assign spherical coordinates to each parcel; centroids were defined as the vertices closest (in terms of Euclidean distance) to the center-of-mass (i.e., arithmetic mean across all the vertices' coordinates) of each parcel. Random rotations were then applied to one hemisphere of this spherical representation of the atlas to disrupt its topography. Rotations were then mirrored to the other hemisphere. Finally, each parcel in the original brain map was reassigned a rotated parcel using the Hungarian algorithm (Kuhn, 1955). In contrast to other reassignment heuristics, this method attempts to reassign each rotated parcel to a unique original parcel, that is, to retain the exact original data distribution. This is particularly important when testing network-based statistics as in the present study . Overall, this procedure was repeated 1,000 times using an openly accessible function from the netneurotools package (https://netneurotools.readthedocs .io/en/latest/index.html) to generate null statistical distributions.
For the correlational analyses, brain maps of communication measures were subjected to this spatial autocorrelation-preserving permutation procedure. For each 1,000 rotated nulls, a correlation coefficient was computed between the surrogate brain map and the statistical map under study (i.e., maps of topological features or Neurosynth functional activation maps), yielding a null distribution of 1,000 coefficients. The original rho was then compared against this null distribution to assess its significance by computing a two-sided p value (p spin ) as the proportion of more extreme null coefficients.
For the partition of the communication measures within and between Yeo intrinsic networks , network labels were first permuted according to the spatially constrained Hungarian method . Dyadic communication measures associated with a node pair belonging to the same intrinsic network, as identified by the permuted labels, were considered within-network, whereas measures associated with a node pair belonging to different intrinsic networks were considered between-networks. The difference between the mean values of the distributions associated with these two categories was then computed for each 1,000 permutation, once again generating a null distribution of this statistical estimate against which the empirical difference was compared to produce a two-sided p value.
A similar method was employed to assess the significance of the differences between mean internal standardized path lengths for pairs of Yeo intrinsic networks.

Yeo Intrinsic Networks
When stratifying brain regions according to their membership in canonical macroscale functional systems, we used the seven intrinsic networks derived by Yeo et al. (2011) via clustering of resting-state fMRI data from 1,000 subjects. A parcellation of the seven resting-state networks in the FreeSurfer fsaverage5 surface space was downloaded from https://github.com /ThomasYeoLab/CBIG/. Nodes of the Cammoun parcellations were then labeled using a winner-take-all approach in which each parcel was attributed the most common intrinsic network assignment across its vertices.

Neurosynth
Functional activation maps synthesizing more than 15,000 published fMRI studies into probabilistic measures of the association between individual voxels and cognitive terms of interest were obtained from the Neurosynth platform (https://github.com/neurosynth/neurosynth; Yarkoni et al., 2011). This association measure quantifies the probability that a particular term is reported in a study if an activation was observed in a given region. The probabilistic maps were extracted from Neurosynth using 123 cognitive terms overlapping the set of keywords from the Neurosynth database and the Cognitive atlas , a public knowledge base of cognitive science. The list of cognitive and behavioral terms ranges from generic concepts (e.g., "attention", "emotion") to specific cognitive processes ("visual attention", "episodic memory"), behaviors ("eating", "sleeping"), and emotions ("fear", "anxiety").
For each of the 123 terms, volumetric reverse inference maps were generated from Neurosynth and projected to the FreeSurfer fsaverage5 mid-gray surface with nearest neighbor interpolation using FreeSurfer's mri_vol2surf function (v6.0.0; https://surfer.nmr.mgh.harvard.edu/). The resulting surface maps were then parcellated according the both the 219 and 1,000 cortical node resolutions of the Cammoun atlas  to obtain node-wise mean probabilistic measures.

Topological Features
In this section, we define the topological features examined in this study and provide their implementation details. All the graph measures, with the exception of closeness centrality, were computed using functions from the open-access Python version of the Brain Connectivity Toolbox (https://github.com/aestrivex/bctpy; Rubinov & Sporns, 2010).
▪ Degree. Binary degree corresponds to the number of edges incident on a node, whereas weighted degree (strength) corresponds to the sum of edge weights incident on a node. ▪ Closeness centrality. The closeness centrality of a node is the inverse of the mean shortest path length between this node and all the other nodes in the network (Rubinov & Sporns, 2010). The closeness centrality C of a node u is defined as where n is the number of nodes in the network and L uv corresponds to the shortest path length between nodes u and v. Note that this measure can be problematic when standardizing shortest path lengths, because of negative contributions to the mean. Nevertheless, we maintain this conventional definition, because all node-wise mean standardized shortest path lengths obtained were positive. ▪ Betweenness centrality. The betweenness centrality of a node is the fraction of shortest paths between all pairs of nodes in the network that contain this node (Freeman, 1977). The betweenness centrality B of a node i is defined as where λ uv is the total number of shortest paths from node u to node v and λ uv (i) is the number of these paths that include node i. Brandes' algorithm (Brandes, 2001) was used to compute the betweenness centrality of each individual node of the structural connectivity matrices under study. ▪ Participation coefficient. The participation coefficient is a measure of the diversity of a node's connection profile among the communities of the network (Guimera & Nunes Amaral, 2005). The participation coefficient P of a node i is defined as where s i,c corresponds to the sum of the weights of all edges incident on node i and nodes in community c, s i is the strength of node i, and C represents the number of communities. A participation coefficient of 0 indicates that the totality of a node's edges are connected to other nodes within its community. The closer the participation coefficient is to 1, the more evenly distributed are its edges among the communities. Note that this measure presumes an established community structure. Here, we used the functional networks derived by Yeo et al. (2011) as communities to compute participation coefficients for each individual node of the structural connectomes under study.

Principal Component Analysis
To take into account a range of models of communication in the brain in addition to shortest path routing, we generated an aggregate measure of the capacity for communication using principal component analysis. First, the standardized communication matrices derived following the procedures detailed above were z-scored across all elements. For each communication matrix, we then extracted all pairwise measures of communication distance, with the exception of undefined diagonal entries, and assembled them into the columns of a matrix X of size d × c, where d is the number of dyads and c is the number of communication models. PCA was then applied to this matrix using the PCA function from the scikit-learn package for machine learning in Python (Pedregosa et al., 2011). First, X was mean centered, that is, demeaned column-wise to obtain X c . Then, a full singular value decomposition (SVD) was applied to X c such that where U and V are orthonormal matrices of sizes d × c and c × c, consisting of the left and right singular vectors, respectively, and S is a diagonal matrix of singular values of size c × c (Eckart & Young, 1936).
In this analysis, we kept only the first component, corresponding to U 1 , the first column of U, which accounted for 62.79% of the total variance in dyadic communication capacity across models. The matrix multiplication U 1 S yielded the principal component scores used in this study.
The same analysis was then reproduced at the node level with standardized node-wise mean vectors of communication distance directly constituting the columns of X. Note that before computing the row averages of the communication matrices, we enforced symmetry of the search information and mean first-passage time matrices by replacing them with (C + C T )/2, where C is the communication matrix and C T is its transpose. In doing so, we consider both the incoming and the outgoing nodal communication capacity. The same procedure was also applied to the matrix of principle component scores prior to the partition of its elements within and between Yeo intrinsic networks.