Connectomic Insights into Topologically Centralized Network Edges and Relevant Motifs in the Human Brain

White matter (WM) tracts serve as important material substrates for information transfer across brain regions. However, the topological roles of WM tracts in global brain communications and their underlying microstructural basis remain poorly understood. Here, we employed diffusion magnetic resonance imaging and graph-theoretical approaches to identify the pivotal WM connections in human whole-brain networks and further investigated their wiring substrates (including WM microstructural organization and physical consumption) and topological contributions to the brain's network backbone. We found that the pivotal WM connections with highly topological-edge centrality were primarily distributed in several long-range cortico-cortical connections (including the corpus callosum, cingulum and inferior fronto-occipital fasciculus) and some projection tracts linking subcortical regions. These pivotal WM connections exhibited high levels of microstructural organization indicated by diffusion measures (the fractional anisotropy, the mean diffusivity and the axial diffusivity) and greater physical consumption indicated by streamline lengths, and contributed significantly to the brain's hubs and the rich-club structure. Network motif analysis further revealed their heavy participations in the organization of communication blocks, especially in routes involving inter-hemispheric heterotopic and extremely remote intra-hemispheric systems. Computational simulation models indicated the sharp decrease of global network integrity when attacking these highly centralized edges. Together, our results demonstrated high building-cost consumption and substantial communication capacity contributions for pivotal WM connections, which deepens our understanding of the topological mechanisms that govern the organization of human connectomes.


INTRODUCTION
The white matter (WM) fibers, i.e., bundles of myelinated axons that provide connections between gray matter (GM) regions, serve as the foundation of information communication in the human brain. Their wiring patterns directly determine the topological performance of brain networks. With recent advances in diffusion magnetic resonance imaging (MRI) techniques and computational tractography methods (Mori et al., 1999), the WM fiber bundles in human brain can be reconstructed in a noninvasive way and the whole-brain structural networks can be generated at a macroscopic level (Hagmann et al., 2008;Gong et al., 2009), with nodes representing GM regions and edges representing the characteristics of WM fiber bundles linking GM regions. The mapping and descriptions of topological organization in human brain WM networks [i.e., "human connectomics" (Sporns et al., 2005)] has led to compelling discoveries of topological properties of brain networks, including their small-worldness (Hagmann et al., 2008;Gong et al., 2009), modular structure (Hagmann et al., 2008;He et al., 2009), highly connected hubs (Sporns et al., 2007;Hagmann et al., 2008;Gong et al., 2009;van den Heuvel and Sporns, 2013), and rich-club organization (van den Heuvel and Sporns, 2011;van den Heuvel et al., 2012). Compared to these global and nodal properties, the topological roles of network edges (i.e., WM tracts) that are responsible for information transfer across regions or systems have been less explored. Several long-range cortico-cortical WM connections, such as the corpus callosum, superior longitudinal fasciculus and cingulum, have been observed to be more frequently involved in efficient information transfer in the network than short-range ones (Gong et al., 2009). From an "edge-centric" perspective, de Reus et al. (2014) demonstrated the effects of the disrupted connections on network topologies and within different network communities. However, the underlying principles for the different topological roles of WM tracts are largely unknown.
What, in terms of wiring substrates and contributions, makes certain WM tracts more topologically important than others? Specifically, (i) do the macro-structural communication capacities of WM connections depend upon their microstructural organization and other physical consumptions such as fiber length? Heavier communicational/information loads flowing in these macro-scale brain networks may impose higher demands on biological resources, reflected by high-level microstructural properties (e.g., high axonal density or compact myelin) and long projection distance (Laughlin and Sejnowski, 2003;Kaiser and Hilgetag, 2006); (ii) How do the pivotal WM connections contribute to the topological properties of network hubs and rich-club structure? The WM tracts with high communication capacity could facilitate the transfer of the massive signals generated and processed by the biologically costly hubs (Vaishnavi et al., 2010;Liang et al., 2013;Tomasi et al., 2013;Collin et al., 2014), and therefore, contributing to the formation of the core structure such as rich-club. (iii) What are the pathmotif patterns in the human brain structural networks in terms of pivotal edges? The minimum communication blocks, defined as the shortest paths between regions, indicate the simplest yet complete information routes in a brain network (Milo et al., 2002;van den Heuvel et al., 2012). Their patterns of utilization of pivotal edges, namely their path motifs, could reflect the information transfer strategy of the human brain.
To address these issues, we collected two scanning sessions of diffusion MRI data in 57 healthy adults and built up the human whole-brain structural networks. We first fully chart the pivotal WM connections, characterized by highly topological centralization, which then allow us to study their critical characteristics. The communication capacity of WM edges in the brain network was quantified using the edgebetweenness centrality measurement (Freeman, 1977;Girvan and Newman, 2002), and highly centralized WM connections were identified as pivotal edges, with their wiring substrates (including WM microstructural organization and streamline length) and topological nexus with brain hubs, rich-club and path motif profiles systematically examined. Validations were performed using two sessions of data, with both low-(90 nodes) and high-resolution (1024 nodes) nodal definitions.

Participants
The data employed in this study were a subset of the Connectivity-based Brain Imaging Research Database (C-BIRD) at Beijing Normal University. This subset includes data from 57 participants (male/female: 30/27; age: 19-30 years) who completed two MRI scan sessions at an interval of approximately 6-weeks (40.94 ± 4.51 days). All participants were right-handed and had no history of neurological or psychiatric disorders. Written informed consent was obtained from each participant, and this study was approved by the Institutional Review Board of the State Key Laboratory of Cognitive Neuroscience and Learning at Beijing Normal University. These data have been released in the Consortium for Reliability and Reproducibility (CoRR) dataset (http://fcon_1000.projects.nitrc.org/indi/CoRR/ html/bnu_1.html; Lin et al., 2015).

Data Acquisition and Preprocessing
All MRI data were obtained using a Siemens Trio Tim 3.0 T scanner (Siemens Medical Systems, Erlangen, Germany) with a 12-channel phased-array head coil in the Imaging Center for Brain Research, Beijing Normal University. Diffusion weighted imaging data were acquired using a single-shot twice-refocused spin-echo diffusion echo-planar imaging sequence. The sequence parameters were repetition time (TR) = 8000 ms, echo time (TE) = 89 ms, 30 non-collinear diffusion directions with b = 1000 s/mm 2 and an additional volume with b = 0 s/mm 2 , data matrix = 128 × 128, field of view (FOV) = 282 mm × 282 mm, 2.2 mm slice thickness, isotropic voxel size (2.2 mm 3 ), bandwidth = 1562 Hz/pixel, and 62 transverse slices without gap covering the whole brain, number of excitation = 2. Three-dimensional high-resolution brain structural images were acquired by using a T1-weighted, sagittal 3D magnetization prepared rapid gradient echo (MP-RAGE) sequences. The sequence parameters were TR/TE/inversion time = 2530 ms/3.39 ms/1100 ms, flip angle = 7 • , FOV = 256 mm × 256 mm, in-plane resolution = 256 × 256, slice thickness = 1.33 mm, and 144 sagittal slices covering the whole brain. The data of session 1 were used for the main analyses and the data of session 2 were used for validation analysis.
For each participant, the diffusion MRI data were preprocessed with FMRIB's Diffusion Toolbox of FSL (Version 5.0; http://www.fmrib.ox.ac.uk/fsl) to correct artifacts induced by head motion and eddy currents by applying an affine alignment of each diffusion-weighted image to the b0 image. Figure 1 illustrates the flowchart of the construction of wholebrain structural networks, which was outlined as follows.

Network Node Definition
In this study, we employed the Automated Anatomical Labeling Atlas (AAL) (Tzourio-Mazoyer et al., 2002) to parcel the brain into 90 cortical and subcortical regions. The detailed procedure for brain WM networks construction has been described previously (Gong et al., 2009). In brief, for each individual, the T1-weighted images were firstly coregistered to the averaged b0 image (two b0 images were obtained for each subject in this study) in native diffusion space using a linear transformation. The transformed T1-weighted images were then nonlinearly transformed to the ICBM152 T1 template in the Montreal Neurological Institute (MNI) space and the transformation matrices were estimated. Finally, the inversed transformation was used to warp the AAL atlases from the MNI space to the native diffusion space, therefore, obtaining the parcellation for each individual. Notably, discrete labeling values in the atlas were preserved by the use of a nearest-neighbor interpolation method. All of these linear and nonlinear mappings were implemented by using the SPM8 package (www.fil.ion.ucl.ac.uk/spm/software/spm8/).

White Matter Tractography
Reconstruction of the whole-brain WM tracts was performed using DTIstudio (version 3.0.3) based on the Fiber Assignment by Continuous Tracking (FACT) algorithm (Mori et al., 1999). Fiber tracking was computed by seeding each voxel with fractional anisotropy (FA) value greater than 0.2. Fiber tracking was stopped at voxels where FA < 0.2 or if the turning angle between adjacent steps was greater than 45 • . All the fiber pathways in the brain were reconstructed using the deterministic tractography method.

Network Edge Definition
To perform the analyses on both the group and individual level, we defined network edges in two ways. For the group-level network, from the set of 57 individual connectivity matrices, a group-level connectivity matrix was computed by selecting all connections that were present in at least 50% of the group of individuals (de Reus and van den Heuvel, 2013a). For the individual networks, we selected a threshold value for the number of streamlines. Two regions were considered to be structurally connected if there are at least 3 streamlines with end points located in these two regions. All of these networks were unweighted. Given the high consistency between the results of group-level and individual level analyses, we mainly reported the results of group-level WM network and treated individual network analyses as a validation.

Network Analysis
Network analyses were carried out by using the Matlab BGL package (www.stanford.edu/~dgleich/programs/matlab_ bgl/) and the GRETNA toolkits (Wang et al., 2015). The 3D visualizations for brain networks were generated using the BrainNet Viewer (Xia et al., 2013).

Identification of Pivotal Edges
We calculated the edge betweenness centrality (EBC) for each edge in the WM network. The betweenness of an edge, B edge i , is a global centrality measures that captures the influence of an edge over information flow between other nodes in the network. B edge i is defined as the number of shortest paths between any pairs of other nodes that pass through the edge (Freeman, 1977;Girvan and Newman, 2002): where, j and k are any two nodes that are not linked directly with edge i in the network, σ j,k is the total number of shortest paths between nodes j and k and σ j,k(i) is the number of shortest paths between nodes j and k that pass through edge i. The normalized EBC was then calculated as follows: > 1, i.e., the EBC > 1 SD above mean) were identified as the pivotal edges in the networks. Additionally, to examine the EBC distribution of the brain networks, we used three possible forms: a powerlaw, P(x) ∼ αx β ; an exponential model P(x) ∼ αexp(βx); and an exponentially truncated power-law, P(x) ∼ αx β exp(x/γ). The cumulative distribution was used to reduce the effects of noise (Strogatz, 2001), and the goodness of fitting was tested using R 2 values.

The Wiring Substrates of the Pivotal Edges
To assess whether the betweenness of an edge depends on its wiring basis, including the WM microstructural properties and streamline length, we calculated the correlations between EBC and each of the WM tract properties and further compared the differences in these measurements between pivotal and non-pivotal edges. Previous research has demonstrated that the microstructural properties of WM such as axonal membrane or myelin can influence the degree of diffusion anisotropy indicated by the indices of diffusion tensor imaging (Beaulieu, 2002). Here, four WM diffusion indices were measured to assess the WM microstructural properties, including the FA, the mean diffusivity (MD), the axial diffusivity (AD) and the radial diffusivity (RD) (Basser, 1995;Song et al., 2002). These four metrics estimate different aspects of diffusion properties of the WM tissues: i) FA expresses the degree of anisotropy of a diffusion process in a given voxel: FA = 3 2 , wherẽ λ = (λ 1 + λ 2 + λ 3 ) 3; ii) MD is often used to estimate the total level of diffusion: MD = (λ 1 + λ 2 + λ 3 )/3; iii) FIGURE 1 | A flowchart for the construction of the whole-brain WM network. (A) The T1-weighted image was firstly rigidly coregistered to the averaged b0 image in native diffusion space; (B) The transformed T1 image was then nonlinearly transformed to the ICBM152 T1 template in the MNI space, and the transformation matrix T was estimated; (C) the inversed transformation T-1 was used to warp the AAL atlas from the MNI space to the native diffusion space, obtaining the parcellation for each individual; (D) the whole-brain WM tracts were reconstructed by using deterministic tractography; (E) the WM fibers connecting each pair of regions were determined for each subject, thus the individual WM networks were constructed; (F) the group-level connectivity matrix was computed by selecting all connections that were present in at least 50% of the group of individuals; and (G) both the individual and group-level networks were further analyzed by using graph theoretical methods.
AD is a metric of the level of diffusion in the direction of the first eigenvector and of local fiber orientation: AD = λ 1 ; and iv) RD is an estimation of the amount of diffusion in the perpendicular to the first eigenvector and of the level of myelination of WM: RD = (λ 2 + λ 3 )/2. For each edge, these four metrics were estimated by averaging the values of the voxels that the streamlines passed through, respectively. The streamline length of each edge, which represents its wiring costs (Kaiser and Hilgetag, 2006;Bullmore and Sporns, 2012), was estimated from the average length of the interconnecting streamlines in each individual network. For the group-level network, all the above metrics were calculated by averaging the values over existing edges across all individuals. Subsequently, Spearman's correlations were used to analyze the correlations between EBC and each of these fiber properties across all edges in the networks. Furthermore, the differences in these fiber properties between pivotal and non-pivotal edges were determined by permutation tests. To evaluate the total performance and consumption of the pivotal edges, we further calculated the proportion of EBC and streamline length of the pivotal edges. This was done by summing up the EBC and streamline length of all the pivotal edges and then dividing them by the total EBC and streamline length of the whole brain, respectively. The proportion curves of EBC and streamline length were plotted to demonstrate whether the pivotal edges have over-average values, in which the x-axis represents the proportion of edges sequenced in EBC and the y-axis represents the accumulate proportion of EBC or streamline length. Furthermore, the cost-performance for each edge was estimated by dividing the EBC by the streamline length and was compared between pivotal and non-pivotal edges.

Contributions to the Nodal Centralities
Considering that the topological properties of nodes and edges are highly interdependent in the brain network, we examined the relationship between the betweenness of the edge and the average nodal centralities of the two nodes it links. These nodal properties included nodal degree, efficiency, and betweenness. The nodal degree K i is a basic topological property, which is defined as the number of links connected to a node: where a ij = 1 if a connection exists between node i and node j in the unweighted network. The nodal efficiency reflects the averaged communication capability of the given node to others, which is defined as the averaged reciprocal of the shortest path length from the given node to other nodes (Latora and Marchiori, 2001): where L −1 ij is the reciprocal of the shortest path length between nodes i and j. The definition for nodal betweenness B node i is similar to edge betweenness, which is defined as the number of shortest paths between pairs of other nodes that pass through the node (Freeman, 1977): where ρ j,k(i) is the number of shortest paths between nodes j and k that pass through node i. Spearman's correlations were adopted to analyze the correlations between EBC and nodal degree, efficiency and betweenness, respectively, across all edges in the networks. The differences on these nodal properties between pivotal and non-pivotal edges were further determined by using permutation tests.

Contributions of the Pivotal Edges to the Rich-Club Structure
The rich-club structure in networks is present when the highdegree nodes of a network tend to be more densely connected among themselves than expected by chance, thus forming a core architecture in the brain network. Such a structure has recently been revealed not only in animal brains (Harriger et al., 2012; Reus and van den Heuvel, 2013b) but also in the human brain (van den Heuvel and Sporns, 2011;Collin et al., 2014). To assess whether the edges belonging to the rich-club had higher EBC and whether the rich-club had more pivotal edges, we examined the rich-club architecture of the WM network. The rich-club coefficient (k) was first calculated as follows: where k is the nodal degree to define hubs, E >k is the number of links among these hubs, and N >k is the number of hub nodes. Then, the (k) was normalized to norm(k) by comparing to the mean (k) of 1000 random networks of equal size and similar connectivity distribution: where random (k) is the averaged rich-club coefficient over the 1000 random networks. An increasing normalized norm (k) > 1 over a range of k reflects the existence of rich-club architecture in a network. In this study, we selected the k where the norm (k) reached the peak value (k = 14) as the threshold for hub definition, which represented the most significant rich-club architecture (we also validated the results of other thresholds, e.g., k > 9, see validation results). Once the rich hub nodes were determined, the edges in the network can be divided into three categories according to the nodes they linked: (i) "rich-club connections" linking rich-club nodes, (ii) "feeder connections" linking rich-club nodes to non-rich-club nodes, and (iii) "local connections" linking non-rich-club nodes to each other (van den Heuvel et al., 2012). Finally, the differences in EBC and proportion of number/EBC of the pivotal edges among these three categories of connections were determined by using ANOVA with post-hoc permutation tests and Chi square tests, respectively.

Communication Length and Path Motifs
To further assess the importance of pivotal edges in brain communication, we analyzed the communication length and path motifs of the WM network. Both metrics were based on the path length between any pair of nodes in the network. Firstly, all 4005 (n × (n − 1) 2, n = 90) unique shortest paths between all 90 nodes in the WM network were traced. Second, the communication length of each shortest path was calculated as the total streamline length of the edges that were used in the path. Subsequently, for each of the shortest paths, the streamline length spent on pivotal or non-pivotal edges was calculated. Finally, once aggregated across all shortest paths, several indices were calculated, including (i) the percentage of paths through pivotal edges, which was calculated by dividing the number of shortest paths passing through pivotal edges by the total number of shortest paths (i.e., 4005); (ii) the percentage of communication lengths of paths through pivotal edges, which was computed by dividing the sum of total streamline length of those shortest paths through pivotal edges by the sum of the total streamline length of all the shortest paths; and (iii) the percentage of communication length spent on the pivotal edges in paths through pivotal edges, which was defined as the ratio between the sum of streamline length pivotal edges of every shortest paths and the total streamline length of all the shortest paths walking through pivotal edges.
The ordered sequence of the pivotal or the non-pivotal edges on the routes of each shortest path were referred to as the "path motifs." Therefore, six patterns of path motifs were identified, including the "N" (non-pivotal), "P" (pivotal), "N-P" (non-pivotal-pivotal), "N-P-N" (non-pivotal-pivotalnon-pivotal), "P-N-P" (pivotal-non-pivotal-pivotal), and "N-P-N-P" (non-pivotal-pivotal-non-pivotal-pivotal). Notably, these path motifs only represented the changes on edge types, but not the exact edge numbers on the path. For example, the "N-P-N" included paths with three edges (NPN), four edges (NNPN, NPPN, and NPNN) or more edges (e.g., NNNPN). Almost all paths could be classified into the six path motifs according to their sequences of edges along the path traveled. Subsequently, the distribution of path motifs was obtained by counting the numbers of the shortest paths in each motif pattern. To assess whether the frequency of each path motif in the WM network was at chance, we generated 1000 matched random networks, identified their pivotal edges, and calculated their path motif distributions. The Z score for each path motif was then computed by subtracting the mean value from the value of proportion in the real WM network and dividing by the standard deviation of those in random networks. Furthermore, we examined the proportions of various path motifs across brain systems to investigate if different brain systems communicated via different path motif patterns. Specifically, all brain regions were first classified into ten different brain systems (i.e., the frontal, parietal, temporal, occipital, and subcortical systems in each hemisphere). Thus, all paths could be allocated into 55 different groups, according to the positions of the two brain regions they connected. Then the proportion of the path motifs for each brain region pair could be measured. The 55 groups could be further classified into four categories, including within system, intra-hemispheric between systems, inter-hemispheric homotopic and heterotopic paths. Finally, a hierarchical clustering analysis was performed to determine whether the four categories had different patterns of path motifs. Briefly, a dissimilarity matrix was calculated by estimating the Euclidean distance between each pair of path groups, and agglomerative hierarchical cluster trees were generated based on the dissimilarity matrix with the weighted linkage agglomerative algorithm.

Vulnerability and Network Robustness
The vulnerability is widely used to quantitatively measure the damage on the network performance caused by the simulated failure of its elements (Costa et al., 2007). To calculate the vulnerability of an individual edge in the WM network, we removed the edges one by one from the network and calculated the changes in global efficiency of the resulting networks. To test the effects of pivotal edges and non-pivotal edges on network performance, we compared the vulnerability values of these two groups using a permutation test.
Network robustness, characterized by the degree of tolerance against random failures and targeted attacks, is usually associated with the stability of a complex network. Here, we investigated the robustness of the networks by the removals of edges (Kaiser et al., 2007;He et al., 2008). To address the random edge failure tolerance of the WM network, we first randomly removed one edge from the networks and then measured the changes in the global efficiency and the size of the largest connected component of the networks. Then, we continued to select and remove additional edges from the networks randomly and recomputed the two measures. To evaluate the targeted attack tolerance, we repeated the above processes but removed the edges with the greatest betweenness. Considering the dynamical compensation in the human brain, we recalculated the edge betweenness after each removal.

Validation Analysis
To evaluate the reproducibility of our findings, we validated our main findings via the following four procedures. Firstly, we parceled the whole brain using a randomly generated high-resolution template with 1024 nodes (Zalesky et al., 2010), and reconstructed the WM network. We repeated our analysis to determine whether our main findings are independent from node definitions. Secondly, we performed the analysis on each of the individual WM networks to assess the reproducibility of those main findings on the individual level. Thirdly, we analyzed the imaging data of the same individuals scanned after an interval of approximately 41 days to determine whether there were consistent topological organizations of the WM network across two scans. Finally, to assess whether our group-level findings are sensitive to network construction thresholds, two additional thresholds were applied. We reconstructed the group-level networks by selecting all connections that were present in at least 40 or 60% of the group of individuals.
All of the comparisons between pivotal and non-pivotal edges were performed using permutation tests. Briefly, for each metric, we initially calculated the difference of the mean values between pivotal and non-pivotal edges. An empirical distribution of the difference was then obtained by randomly reallocating all of the values into two groups and re-computing the mean differences between the two randomized groups (10,000 permutations). The original difference between pivotal and non-pivotal edges was assigned a p value as the proportion of random values in the obtained empirical distribution. The 95th percentile points of the empirical distribution were used as critical values in a one-tailed test of whether the observed group differences could occur by chance.

Results
We used diffusion MRI tractography approaches to construct the brain networks at both the group and the individual levels (Figure 1). The group-based brain network with 90 nodes was fully connected, with 431 WM connections and a connection density of 10.8%. For the individuals, the densities of the brain networks ranged from 8 to 12% (mean ± SD: 9.7% ± 0.9%). We further identified the pivotal WM connections (i.e., highly centralized edges) from the brain networks and systematically examined their physical characteristics and topological contributions to network communications. Given compatible results between 90-node and 1024-node networks and between the group-and individual-level networks, we mainly report the findings from the group-based network analyses with 90 nodes, unless specifically mentioned.

Edge Betweenness Distribution
We used the EBC (Freeman, 1977;Girvan and Newman, 2002) to quantify the topological centrality or communication capacity of the WM edges in the structural brain network. The EBC of an edge measures the frequency with which the shortest path between any region pair passes through this edge. The EBC distribution of the network was best fitted by the exponentially truncated power-law form [P(x) ∼ αx β exp(x/γ)] rather than the power-law [P(x) ∼ αx β ] and or exponential [P(x) ∼ αexp(βx)] models (Figure 2A). The estimated parameters were: α = 0.98, β = 0.20, γ = −15.42 and R 2 = 0.998, respectively. Such a model indicates that i) the WM edges play heterogeneous roles in information communication across the network, and ii) the brain network includes some highly centralized edges but prevents the existence of extremely centralized WM connections with overly heavy loads.

Identifying Pivotal Edges
The WM edges with higher EBC (>1 SD above the whole-brain mean) are referred to as pivotal edges in the present study.
Forty-eight of the 431 edges (11.1%) were identified as pivotal edges in the whole-brain WM network ( Figure 2B and Table 1), including 16 inter-hemispheric and 32 intra-hemispheric (22 inter-and 10 intra-lobe) connections. The spatial distribution pattern of the group-based pivotal edges was largely consistent with the probability map of pivotal WM edges across individuals ( Figure 2B). Specifically, these pivotal edges were primarily located in human major WM tracts, involving the corpus callosum (16/48), the cingulum (7/48), the uncinate fasciculus (4/48), and several projection tracts linking the subcortical with cortical regions (5/48) ( Figure 2C and Table 1).

The Wiring Substrates of the Pivotal WM Edges
We further explored whether topologically centralized network edges involve high-level microstructural organization and expensive physical consumption. The WM diffusion properties, including the FA, MD, and AD, exhibited significantly higher values in the pivotal edges compared to the non-pivotal ones (all ps < 0.0004, permutation tests) ( Figure 2D). Notably, the RD did not show significant difference in EBC between the pivotal and non-pivotal edges (p = 0.41, permutation test). The pivotal WM edges exhibited significantly greater streamline lengths than the non-pivotal ones (p < 0.0001, permutation test) (Figure 2E), and the EBC was positively correlated with streamline length across edges (Spearman ρ = 0.29, p < 0.0001) ( Table S1), suggesting that these centrally embedded WM connections tend to span longer physical distances. The ratio of EBC to streamline length, quantifying the communication capacity per unit of streamline length or cost-performance, was significantly higher in the pivotal edges than the non-pivotal edges (p < 0.0001, permutation test) ( Figure 2E). Lastly, we charted the curves for the proportion of edges vs. proportions of EBC and streamline length, and found that the pivotal edges (10.8% in number) consumed 16.9% of the streamline length but contributed 31.3% to the total communication capacity (in terms of EBC) of the whole brain ( Figure 2F). These results together suggest the costly but highly cost-efficient signature of the pivotal edges.

Contribution to the Nodal Properties
We explored the relationship between the EBC values of WM edges and their linked nodes' properties (nodal degree, efficiency, and betweenness). The pivotal edges had significantly greater contributions to all three nodal properties than the non-pivotal ones (all ps < 0.0001, permutation tests) (Figure 3). Moreover, significant positive correlations were found over all three nodal centralities and across all edges, with the Spearman's correlation coefficients of 0.35, 0.38, and 0.59, respectively (all ps < 0.0001) ( Table S1). These results indicated a strong topological nexus between pivotal edges and pivotal nodes in the WM network.

Contribution to the Rich-Club Architecture
We examined the contribution of pivotal WM edges to the richclub architecture of the WM network. Figure 4A illustrates the curve of the normalized rich-club coefficient, norm (k), over a range of nodal degree, k. The norm (k) were larger than 1 at k > 8, indicating the existence of a rich-club structure in the WM network (van den Heuvel and Sporns, 2011). Here, we chose the peak norm (k), where the nodes with k > 14 were considered the brain hubs (see Table S2 for results of other thresholds). We identified 11 network hubs that were primarily located in the bilateral precuneus, the bilateral orbital part of superior frontal gyrus, the right dorsolateral superior frontal gyrus, the bilateral calcarine sulcus, the left middle occipital gyrus, the right superior occipital gyrus and the bilateral putamen ( Figure 4B). Most of these hubs (72.7%) were structurally connected with the pivotal edges. Further, we divided the whole-brain WM connections into three categories according to the types of nodes they linked (van den Heuvel et al., 2012): rich-club connections between rich-club nodes (n = 21), feeder connections between rich-club and nonrich-club nodes (n = 142) and local connections between two non-rich-club nodes (n = 268). Significant differences in EBC were observed among these three edge categories [F (2, 428) = 44.9, p = 2.0 × 10 −18 ], with a descending order of the rich-club, feeder and local connections (post hoc comparisons, permutation tests, all ps < 0.001) ( Figure 4C). Importantly, the proportions of the number of pivotal WM edges among the three categories, which represents the network building contribution, were significantly different [χ 2 (2) = 73.5, p = 1.1 × 10 −16 ]: 57.1% (12/21) within the rich-club connections, 19.7% (28/142) within the feeder connections and 3.0% (8/268) within the local connections ( Figure 4D). The proportion of EBC of the pivotal edges among the three categories, which represents the network communication contribution, was also significantly different [χ 2 (2) = 2490.3, p < 1.0 × 10 −64 ]: 83.8% within the rich-club connections, 41.4% within the feeder connections and 11.3% within the local connections ( Figure 4E). These results suggest that the rich-club architecture of the brain networks was topologically supported by the pivotal WM edges.
FIGURE 3 | Contribution of the pivotal edges toward the centrality of the nodes they linked. The pivotal edges had significantly greater contributions to all three nodal properties (nodal degree, efficiency, and betweenness) than the non-pivotal ones. The contribution toward nodal centralities of an edge was estimated by averaging nodal properties of its two linking nodes.

FIGURE 4 | Pivotal edges and the rich-club structure. (A)
The normalized rich-club coefficient norm (k) of the group-level WM network was above 1 for a range of k from 9 to 16. The peak at k = 14 was selected as the hub threshold for further analysis. (B) The network hubs were mainly located in the medial line of the brain and the connections of the brain network can be further classified into three categories: rich-club (red), feeder (yellow) and local (blue) connections.

Communication Length of Pivotal Edges
When investigating the shortest paths (i.e., the minimum communication block) between any two nodes in the brain network, 58% of the paths (4005 in total) were found to pass through the pivotal edges. The pivotal-edge related paths accounted for 66% of the total communication length (communication length of one shortest path was defined as the total streamline length of the edges along the path). Specifically, when considering only the pivotal edge related paths, the total streamline length of these pivotal edges accounted for 67% of the total communication length, suggesting their high utilizations in brain communication.

Path Motifs of the WM Network
Every shortest path between nodes walks through a series of edges, of which the ordered sequence of the pivotal or the ordinary edges on the routes was referred to as the "path motifs." Six types of path motif were identified in the wholebrain WM network (Figure 5A), and their appearing frequencies were statistically compared to those of 1000 equivalent random networks. The comparison revealed that the paths in the brain network exhibited significantly greater percentages of several pivotal-edge related connection types (e.g., "P, " "N-P, " "N-P-N, " and "N-P-N-P, " all Zs > 4.4, Figure 5B), especially the nonpivotal to pivotal to non-pivotal (N-P-N) path motif (Z = 23.1). Such a path-motif distribution profile indicates the central role of pivotal edges in the communicational organization of brain circuits.

Path Motifs across Different Brain Systems
Hierarchical clustering analysis of the path motifs across different brain systems revealed four distinct patterns, with gradual changes from the topologically simplest path motif "N" to complex "P" related patterns ( Figure 5C): (i) the paths within the frontal and occipital systems in each hemisphere, and several occipital related intra-hemispheric paths mostly travel through the non-pivotal edges (mean percentage of "N": 89.2%), indicating a plainest communication pattern between these regions; (ii) the pattern with a moderate percentage of motif "N" (61.0%) and a small percentage of motif "N-P" (26.3%) was observed in the paths within the parietal, temporal and subcortical systems, and in most of the intra-hemispheric paths; (iii) the pattern with fewer motif "N" (35.8%) but more "P" related motifs [57.3%, including "N-P" (34.8%) and "N-P-N" (22.5%)] were primarily distributed in the homotopic paths between bilateral frontal systems and heterotopic paths related to frontal and occipital systems; and iv) the last pattern with mostly the "P" related motifs [88.4%, including "P" (11.6%), "N-P" (47.5%), and "N-P-N" (29.3%)] primarily existed in the homotopic and heterotopic paths among bilateral parietal, temporal and subcortical regions. These findings suggest that the path motifs are distinctively differentiated in their communication across different brain systems: the heterotopic systems between hemispheres and extremely distant intrahemispheric systems (e.g., between frontal and occipital systems) tend to adopt the pivotal related paths, while other intrahemispheric systems tend to use plain communication patterns.

Lesion Simulations
Two simulation strategies were used to evaluate how a "lesion" of the pivotal WM edges influenced brain network performance. First, we calculated the vulnerability of each edge in the WM network, which was defined as the change in the global efficiency of the network after eliminating this edge from the whole-brain network (Costa et al., 2007). Not surprisingly, the pivotal edges showed significantly greater vulnerability than the non-pivotal ones (p < 0.0001, permutation tests) (Figure 6A), especially the precuneus-related fiber tracts (Table 1). Second, we evaluated the topological robustness of the brain networks against random failure and targeted attacks of WM edges (Kaiser et al., 2007;He et al., 2008). Both the global efficiency and the size of the largestconnected component slowly declined in response to the random failure; In contrast, these network properties decreased rapidly in response to a targeted attack, with an over 40% reduction when 20% of the most centralized edges were attacked (Figures 6B,C). These simulation analyses highlight the topological significance of the pivotal WM edges in the global integrity of brain networks.

Validation Results
We evaluated the reproducibility of our main findings in several different ways, involving thresholds for rich-club, high-resolution brain parcellations, analyzing individual WM networks, data from another session and different connectivity thresholds during network constructions. We found the main results remained unchanged (Supplemental Results, Figures S1-S3, Tables S2-S6), indicating a robust reliability of our findings.

DISCUSSION
We mapped the small proportion of WM connections in human whole-brain networks that were highly topologically centralized in terms of global brain communication. These pivotal WM connections exhibited higher levels of WM microstructural organization and consumed longer streamline lengths, while they topologically contributed significantly to the brain's hub and rich-club architecture as well as the communication blocks, especially in the routes between inter-hemispheric and extremely remote intra-hemispheric systems. Simulation models showed that the integrity of brain networks would decrease sharply when pivotal edges were under attack. These signatures of pivotal WM connections were highly reproducible across individuals, scan sessions, spatial scales and network construction strategies.

Spatial Distribution of Pivotal WM Connections
We found that different WM tracts play heterogeneous roles in global information communication in brain networks. The most centralized WM connections were primarily located in several major long-range WM pathways, such as the corpus callosum, cingulum and inferior fronto-occipital fasciculus, which are the vital communication lines across regions to support human cognition. For instance, one third of the pivotal edges belonging to the inter-hemispheric connections located in the corpus callosum, which coordinates numerous functional integrations processes including perception, attention, memory, language and reasoning (Gazzaniga, 2000). The cingulum that connects medial frontal, parietal and temporal systems, is the principal route of the default mode network (Greicius et al., 2009), damage to which is associated with various cognitive impairments (Metzler-Baddeley et al., 2012). Several projection tracts are also topologically centralized, involving fibers linking the subcortical nuclei such as the putamen, the thalamus and the caudate nuclei, which are hubs of the human WM networks FIGURE 5 | Path motifs of the whole-brain WM network and of different brain systems. (A) Examples of the six types of path motifs in the whole-brain WM network. Path motifs were defined as the ordered sequence of the pivotal or the non-pivotal edges on the routes of each shortest path. N, non-pivotal edge; P, pivotal edge. (B) The frequency percentage and normalized distribution of path motifs derived by comparing the actual frequency of each path motif to that of 1000 equivalent random networks. The "non-pivotal to pivotal to non-pivotal" (N-P-N) path motif was the most frequent path motif in the brain network (Z = 23.1). (C) The bottom matrix shows the proportions of path motifs between each pair of the brain systems. The following hierarchical clustering analysis revealed four path-motif distribution patterns, of which the path motif "N" decreased gradually while the "P" related path motifs constantly accumulated. Notably, most within-system and intra hemispheric paths communicated with the style "N," while the inter-hemispheric paths, especially paths between heterotopic systems, utilized the pivotal edges more often. Fro, frontal; L, left; Occ, occipital; R, right; Par, parietal; Tem, temporal; Sub, subcortical; IntraHemi, intra hemispheric; InterHemi, inter hemispheric. (van den Heuvel and Sporns, 2011;Crossley et al., 2014) and are involved in various brain functions, (Burgess et al., 2002;Packard and Knowlton, 2002;Kunimatsu and Tanaka, 2010). From an anatomic embedding perspective, these WM tracts, especially the posterior cerebral WM, have also been demonstrated to be more involved in the construction of the brain network as the voxels within these regions had more network edges passing through them (Owen et al., 2015). Together, these pivotal WM connections provide critical high-throughput communication channels or shortcuts between cortical and subcortical regions.

The WM Microstructural Organization and Wiring Cost of the Pivotal WM Connections
The pivotal WM connections showed higher levels of microstructural organization, as indicated by greater FA, MD, and AD values than non-pivotal ones, which was mainly driven by the higher AD value in pivotal WM tracts. Such results might reflect more orderly fiber organization, greater axonal diameter, larger packing densities, higher proportion of myelinated axons of these pivotal edges (Basser, 1995;Beaulieu, 2002). The pivotal edges also showed greater physical consumption, indicated by streamline length, suggesting an extraordinary wiring cost for building these topologically centralized shortcuts that facilitate the link across distant brain regions. Intriguingly, the unit communication capacity (dividing EBC by streamline length) is higher for pivotal edges, indicating their rewarding building cost. Together, the high-level microstructural organization and longer axonal fibers empower these pivotal edges to have faster communication routes with shorter transmission delays, which consequently facilitate synchronous information processing, increase signal transfer robustness and reduce noise during communication (Laughlin and Sejnowski, 2003;Kaiser and Hilgetag, 2006;Collin et al., 2014). This phenomenon, wherein a small number of high-quality and long-range pivotal edges consume substantial resources to ensure high levels of both local and global information integration of the brain, provides additional support for the concept of cost-efficiency balance of neural circuitry formation (Bullmore and Sporns, 2012).

Contributions of the Pivotal Edges toward the Network Hubs/Rich-Club Architectures
The pivotal edges had a significantly greater contribution to nodal topological properties than the non-pivotal ones. Biologically costly hubs are supposed to generate and process massive signals to achieve the functional integration (Sporns et al., 2007;Buckner et al., 2009;Liang et al., 2013;Tomasi et al., 2013;van den Heuvel and Sporns, 2013). Nearly 80% of pivotal WM connections were directly linked with hubs. These pivotal WM connections with high-throughout transfer capability could fulfill the communicational requirement of the hubs and further nourish their developments. In fact, a significant positive correlation between nodal centrality and the anatomical distance of edges has been reported (Alexander-Bloch et al., 2013;Crossley et al., 2014). The associations between the edge-betweenness centrality and nodal centralities observed here provide more direct evidences for the interactive nexus between pivotal edges and brain hubs. Furthermore, approximately 20% of the pivotal WM edges connect non-hub regions. Such organization profiles avoid the over-aggregation of pivotal WM connections and formation of dominating clusters that shoulders extremely large communication loads. While clusters with extremely huge loads could be consumptive and vulnerable to pathogenic processes, the distributed localization of pivotal edges may help launch a compensation mechanism in cases of disease.
Recent studies have identified significantly denser connections between hub regions compared with non-hub regions, forming a rich-club core architecture in the structural brain networks in humans (van den Heuvel and Sporns, 2011) and other species (Harriger et al., 2012;de Reus and van den Heuvel, 2013b;Towlson et al., 2013). The rich-club organization is biologically costly in terms of metabolism and wiring cost but provides functional benefits by enhancing both the global information flow and the resilience of the network to hub attacks (van den Heuvel and Sporns, 2011;van den Heuvel et al., 2012). Here, from a perspective of network edges, we showed significant contributions of the pivotal WM connections to either the building number or the communication load within the richclub architecture, suggesting that the pivotal WM connections facilitate the communication in the rich-club. Notably, some pivotal connections were those feeder connections linking hubs and non-hubs, reflecting their key roles in guiding and shunting of information flow around the communication cores.

Utilization of the Pivotal WM Connections in Communication Strategy
One particularly intriguing finding here is that the communication blocks (i.e., the shortest paths) in the brain network are mainly associated with pivotal edges: the popular communication motif follows a route sequence of edges, where signals pass through increasingly centralized edges and then decreasing centrality (following an "N-P-N" path motif), indicating the pattern of first traveling on a side road, then turning onto highways, and finally leaving the arterial traffic. Interestingly, van den Heuvel et al. (2012) examined the relationship between the shortest communication paths and the rich-club architecture, and revealed a "zoomingout/zooming-in" structure of shortest paths by which the signals fed into, traversed, and exited the rich-club. They contended that this pattern is an expression of the degree-based "greedy routing" (Kleinberg, 2000) navigation strategy, an efficient communication scheme in technology and transportation networks in which the travel paths are selected only on the basis of local information without the knowledge of the global topography of the network (Simsek and Jensen, 2008;Boguna et al., 2009). Information about the global brain structure is probably absent from the local neurons or brain regions; therefore, sending information to distant targets by traveling through the topologically centralized WM connections for their spreading possibility could be an optimal path. Notably, we also observed diverse communication path motifs across different brain systems: inter-hemispheric paths communicating brain systems are more dependent on pivotal WM connections than intra-hemispheric within-or between-system paths; path motifs of the parietal or subcortical regions are far more complex than those of the systems with homogeneous and/or primary functions (e.g., occipital). These patterns inspire new questions about exactly how specific path motifs are related to the types of functions of a particular system.
In summary, our findings provide crucial evidence for the mechanisms underlying efficient communication strategies in human WM networks: investing in a "highway" system and utilizing it in the support of hubs, rich-club structures, and most-prevalent path-motifs. Such findings demonstrate a specific manifestation of the cost-efficiency principle in brain WM networks and deepen the understanding of signal exchange patterns across regions.

Limitations and Further Considerations
Several methodological limitations warrant further consideration. First, diffusion tensor imaging in associated with inaccuracies in resolving crossing fibers and sharp angulations of tracts , which could lead to false-negatives in tracing for long-range fibers and false-positive connections between nearby regions. We applied a higher FA threshold (0.2) in fiber tracing to minimize possible artifacts due to the acquisition and tractography noise. However, the unbalanced numbers of long-and short-range connections might overestimate the importance of long-range connections, and result in the increase of FA of pivotal edges due to the fact that longer edges pass through high myelinated central white matter tracts. Future studies could be conducted by collecting diffusion spectrum imaging or high angular resolution diffusion imaging and reconstructing WM pathways with complex geometries, despite of the noise due to the high b-values (i.e., strong gradient) and long acquisition times. Second, the diffusion properties, including the FA, MD, AD, and RD, were estimated to represent the microstructural organization of the pivotal edges. Although these metrics provide approximate reflections of specific microstructural attributes such as fiber organization, axon density and myelination, the accurate biological interpretation for these indices remains ambiguous and controversial (Jones et al., 2013). Future studies using advanced diffusion MR models such as CHARMED (Assaf et al., 2004;Assaf and Basser, 2005), AxCaliber (Assaf et al., 2008), and ActiveAx (Alexander et al., 2010), and other MR contrast mechanisms such as quantitative magnetization transfer imaging (Sled and Pike, 2001;Cercignani and Alexander, 2006) might provide more accurate microstructural information. Third, relating the functional characteristics (e.g., dynamics and causality) (Hiltunen et al., 2014;Zalesky et al., 2014) of pivotal edges in functional networks to the pivotal WM connections is important for understanding the possible mechanisms of structural and functional coupling. Fourth, incorporating computational modeling approaches (Honey et al., 2009;Raj et al., 2012;Vertes et al., 2012;Deco et al., 2014) into the analysis of pivotal WM connections can further emulate their topological properties in an empirical situation. Conversely, optimizing the model parameters by taking into account the physical and topological characters of pivotal WM connections might help design artificial networks. Finally, a simulation "lesion" analysis was performed here to evaluate the topological influences of the pivotal edges on global brain communication. Real disease data are desirable to ascertain how brain lesions located in pivotal WM connections affect the topological performance of brain networks and subsequently cognitive dysfunctions.

AUTHOR CONTRIBUTIONS
MX, YB, and YH: designed research; MX, QL, YB, and YH: performed research; MX, QL, and YH: contributed analytic tools; MX, QL, and YH: analyzed data; and MX, YB, and YH wrote the paper.