Abstract
The human brain is a complex network of interlinked regions. Recent studies have demonstrated the existence of a number of highly connected and highly central neocortical hub regions, regions that play a key role in global information integration between different parts of the network. The potential functional importance of these “brain hubs” is underscored by recent studies showing that disturbances of their structural and functional connectivity profile are linked to neuropathology. This study aims to map out both the subcortical and neocortical hubs of the brain and examine their mutual relationship, particularly their structural linkages. Here, we demonstrate that brain hubs form a so-called “rich club,” characterized by a tendency for high-degree nodes to be more densely connected among themselves than nodes of a lower degree, providing important information on the higher-level topology of the brain network. Whole-brain structural networks of 21 subjects were reconstructed using diffusion tensor imaging data. Examining the connectivity profile of these networks revealed a group of 12 strongly interconnected bihemispheric hub regions, comprising the precuneus, superior frontal and superior parietal cortex, as well as the subcortical hippocampus, putamen, and thalamus. Importantly, these hub regions were found to be more densely interconnected than would be expected based solely on their degree, together forming a rich club. We discuss the potential functional implications of the rich-club organization of the human connectome, particularly in light of its role in information integration and in conferring robustness to its structural core.
Introduction
The human brain is a network. At the large scale, its complex wiring diagram forms a network of hundreds of brain regions and thousands of interconnecting white matter axonal pathways (Sporns, 2011). Network studies of the human brain are motivated by the idea that brain function is not solely attributable to individual regions and connections, but rather emerges from the topology of the network as a whole, the connectome of the brain (Sporns et al., 2005; Bullmore and Sporns, 2009).
Looking at the brain as an integrative complex system, one principal marker for its structural ordering is its degree distribution, expressing the probability P(k) of a node to share a link with k other nodes in the network (Barabasi and Bonabeau, 2003). Studies of brain networks have shown a heavy-tailed degree distribution (Eguíluzet al., 2005; Hagmann et al., 2007; van den Heuvel et al., 2008b), indicating the existence of a number of highly connected regions (Hagmann et al., 2008). It has been noted that some of these brain regions play a central role in the overall network organization, as indexed by a high degree, low clustering, short path length, high centrality and participation in multiple communities across the network, identifying them as “brain hubs” (Sporns et al., 2007; van den Heuvel et al., 2010). Examining the function and role of these hubs is of special interest as they play a central role in establishing and maintaining efficient global brain communication, a crucial feature for healthy brain functioning (Bassett et al., 2009; van den Heuvel et al., 2009a). First studies have identified a number of key cortical hubs (Hagmann et al., 2008; Gong et al., 2009; van den Heuvel et al., 2010), but many organizational properties of brain hubs—particularly their structural linkages—have yet to be revealed.
The so-called “rich-club” phenomenon in networks is said to be present when the hubs of a network tend to be more densely connected among themselves than nodes of a lower degree (Colizza et al., 2006). The name arises from the analogy with social systems, where highly central individuals—being “rich” in connections—often form a highly interconnected club (Zhou and Mondragon, 2004; McAuley et al., 2007). The presence, or absence, of rich-club organization can provide important information on the higher-order structure of a network, particularly on the level of resilience, hierarchal ordering, and specialization (Colizza et al., 2006; McAuley et al., 2007). The strong rich-club tendency of power grids, for example, is related to the necessity of the network to easily distribute the load of one station to the other stations, reducing the possibility of critical failure, while the absence of rich-club organization in protein interaction networks has been suggested to reflect a high level of functional specialization (Colizza et al., 2006; McAuley et al., 2007). This study was aimed at examining the rich-club organization of the human connectome. Diffusion tensor imaging (DTI) was used to reconstruct the connections of the structural brain network, and rich-club organization was explored by taking into account both the density and efficacy of the white matter connections of the brain.
Materials and Methods
Subjects
A group of 21 healthy participants with no history of neurological or psychiatric disorders were included in this study [mean age (SD), 29.95 (8.3) years; 15 males, 6 females], after obtaining written informed consent, in line with the Declaration of Helsinki approved by the medical ethics committee for research in humans of the University Medical Center Utrecht (Utrecht, The Netherlands).
MR acquisition
MR data was acquired on a 3 tesla Philips Achieva clinical scanner at the University Medical Center Utrecht using an 8 channel SENSE receiver head coil. Data acquisition included high-angular DTI and an anatomical T1 image for anatomical reference. Within each 30 min scanning session, two DTI sets were acquired, each consisting of 30 weighted diffusion scans and 5 unweighted B0 scans. Acquisition parameters for the DTI-MR included the following: parallel imaging SENSE with a p reduction of 3; high angular gradient set of 30 different weighted directions (Jones et al., 1999; Jones, 2004); TR/TE, 7035/68 ms; 2 × 2 × 2 mm; 75 slices covering the whole brain; b weighting of 1000 s/mm2 (van den Heuvel et al., 2008a, 2009b). The second set was acquired with similar acquisition parameters, except for a reversed k-space readout. Second, directly after the acquisition of the diffusion-weighted scans, a T1-weighted image was acquired for anatomical reference and definition of the different structural nodes of the network, using the following scanning parameters: 3D fast field echo using parallel imaging; TR/TE, 10/4.6 ms; FOV, 240 × 240 mm; 200 slices covering the whole brain; 0.75 mm isotropic voxel size.
DTI preprocessing and deterministic fiber tracking
DTI preprocessing included the following steps. First, DTI images were corrected for possible susceptibility distortions, often reported in single-shot EPI DTI images, by estimating a field distortion map based on the two b = 0 images. The distortion map was computed using the information of the two b = 0 images, which were acquired with a reversed k-space readout (Andersson et al., 2003). This distortion map was applied to the two sets of 30 weighted images (Andersson et al., 2003), resulting in a corrected DTI set, consisting of a single unweighted corrected b = 0 image and 30 corrected weighted images. Second, diffusion images were corrected for eddy-current distortions (Andersson and Skare, 2002) and realigned to the b = 0 image. Third, a tensor was fitted to the diffusion profile within each voxel using a robust tensor fitting method (Chang et al., 2005), and the preferred diffusion direction within each voxel was computed as the principal eigenvector of the eigenvalue decomposition of the fitted tensor. Fourth, based on computed eigenvalues, the level of fractional anisotropy (FA) of each voxel was computed, providing information on the level of preferred diffusion direction within a voxel (Beaulieu and Allen, 1994). Fifth, using the extracted information on the preferred diffusion direction with each voxel in the brain mask, the white matter tracts of the brain network were reconstructed using the deterministic fiber tracking, based on the FACT (fiber assignment by continuous tracking) algorithm (Mori et al., 1999; Mori and van Zijl, 2002). Within each voxel in the brain mask (i.e., the total of white and gray matter), eight seeds were started, evenly distributed over the volume of the voxel. A streamline was started from each seed following the main diffusion direction from voxel to voxel, thus reconstructing white matter fibers. This procedure resulted in a large sample of all possible (reconstructable) fiber tracts of the brain network. A fiber streamline was stopped when the fiber track reached a voxel with a FA value <0.1 (indicating a low level of preferred diffusion within that particular voxel), when the trajectory of the traced fiber left the brain mask or when the fiber tract made a sharp turn of >45°.
T1 preprocessing
Within each individual dataset, the T1 images were used for anatomical reference and for the selection of the nodes of the brain network. Gray/white matter tissue classification was performed on the basis of the T1-weighted MR image using the well validated Freesurfer suite (V5; http://surfer.nmr.mgh.harvard.edu/) (Fischl et al., 2004). In addition, the Freesurfer suite was used for automatic segmentation of subcortical structures—thalamus, pallidum, caudate, putamen, accumbens, hippocampus, amygdala—and automatic parcellation of the reconstructed cortical surface into 68 distinct brain regions. In total, 82 brain regions were selected, representing the nodes of the individual brain networks.
Reconstruction of structural brain networks: unweighted, streamline-weighted, FA-weighted
Reconstruction of individual networks.
Individual brain networks were modeled based on the set of reconstructed fiber tracts F combined with the collection of segmented brain regions (Hagmann et al., 2008; van den Heuvel et al., 2010). The different steps are illustrated in Figure 1. A network can be mathematically described as a graph G = (V,E), consisting of a collection of nodes V and a set of connections E. V was taken as the collection of 82 distinct segmented brain regions (Fig. 1II). Next, within each individual dataset, for each pair of brain regions i and j it was determined whether there were fibers present in the total collection of reconstructed tracts F that interconnected region i and region j. This was done by evaluating for all fibers in F whether a fiber touched both region i and region j. In total, this procedure resulted in (on average) the use of 53% of the total collection of reconstructed fiber streamlines [group mean (SD), 53 (2.3)]. Forty-seven percent of the total collection of fiber streamlines failed to connect a pair of brain regions, resulting from exclusively passing through white matter or touching only a single brain region. Next, when a connection was present between regions i and j, information on this connection was included in the connectivity matrix M. Four variations of M were computed as follows.
(1) Unweighted Munw.
In the unweighted brain networks connections were represented in a binary fashion. When a streamline between region i and region j was present, the connection was included in the connectivity matrix Munw, by setting the value of cell Munw(i,j) to 1 and zero otherwise. This procedure was then repeated for all nodes i and j in V generating an unweighted undirected connectivity matrix M, representing the presence or absence of white matter connections of the individual structural brain networks.
(2) Weighted by streamline count Mw-nos.
In addition to an unweighted graph approach, in which only information on the existence of a connection is present, one can include information on the strength of the existing connections, including additional information on the quantity and/or efficacy of the connections. To this end, in addition to Munw, information on the quantity of white matter connectivity between region i and region j was computed, setting cell (i,j) of a weighted connectivity matrix Mw-nos to the sum of the existing streamlines that resulted from the fiber selection procedure.
(3) Weighted by streamline count, corrected for ROI volume Mw-nosROI.
However, as across the set of regions, the volumes of the segmented regions are not uniform, the size of a ROI may influence the fiber selection procedure: bigger regions may have a higher probability of being touched by one of the fiber streamlines. Therefore, to control for this effect, a second NOS-weighted Mw-nosROI was computed, in which, during the formation of the connectivity matrix M, the number of streamlines between region i and j was normalized by the sum of the volumes of ROI i and j. However, it has been suggested that volume correction may potentially overcompensate volume-driven effects on the streamline count. In addition, some studies have pointed out the effect of a distal bias in streamline tractography, as the number of fiber streamlines between two regions of the network decreases as a function of the distance that separates them (Zalesky and Fornito, 2009). Longer streamlines include a greater number of propagation steps, each corresponding to a new decision point in the streamline procedure. The degree to which different normalization strategies may interact, and to which extent they may bias the data to one end or the other, is unknown. It is for these reasons that we choose to focus mainly on the unweighted and NOS-weighted results.
(4) Weighted by FA values Mw-FA.
Information on the integrity of the interregional white matter connections was incorporated by estimating FA values for each of the interconnecting tracts. FA values express the level of anisotropic diffusion of white matter in a brain voxel and are a commonly used metric to examine the microstructural aspects of brain connectivity. It is believed that a large contribution to the direction-dependent diffusion signal comes from axonal membranes hindering the diffusion of water molecules, with higher FA values expressing a higher level of microstructural organization of white matter connections (Beaulieu, 2002). Lower levels of FA are commonly used as a marker for white matter damage in patient studies (Sun et al., 2003; Kubicki et al., 2005; Verstraete et al., 2010). Furthermore, higher FA values of white matter tracts have been linked to faster task performance (Ewing-Cobbs et al., 2006; Gold et al., 2007), making FA values a possible marker for the efficacy of brain connections. To this end, the points along each fiber were colored with their corresponding voxelwise FA values, and for each of the existing connections between region i and region j the average FA was computed as the mean of the FA values of all included streamlines that formed the connection between region i and j. The mean FA value was then entered in the FA-weighted connectivity matrix Mw-fa.
Group connectomes: weighted and unweighted
For both the unweighted and weighted networks, a group-averaged network was computed. (1) For unweighted, from the set of 21 individual connectivity matrices M, a group average connectivity matrix Mgroup was computed by selecting all connections that were present in at least 75% of the group of subjects. (2, 3) Next, across the individual weighted connectivity matrices Mw-nos/Mw-nosROI, a group-averaged weighted Mgroup was computed by averaging the cell values of the individual matrices for those connections present in the unweighted matrix M in the subgroup of subjects (i.e., >75%, leaving out zero entries). (4) Using a similar approach, a group-averaged weighted Mgroupw-fa was computed, by averaging the FA values of the connections over the group of subjects.
Rich-club organization
Graph metrics.
Graph theory was used to examine the topology of the reconstructed brain networks. Characteristic measures of network organization were computed, including the (node-specific) degree k, clustering coefficient, characteristic path length, betweenness centrality, normalized clustering coefficient and normalized path length (both normalized relative to a set of 100 comparable random graphs), global efficiency, assortativity, and modularity. Graph metrics were computed using the Brain Connectivity Toolbox as described previously (Rubinov and Sporns 2010). The centerpiece of the present paper is the investigation of rich clubs of hub nodes, which we examine by computing the rich-club coefficients Φ(k) across a range of degree k of the networks.
Unweighted network.
For a given M, the degree of each node i in the network was determined, counting the number of links that node i shared with k other nodes in the network. All nodes that showed a number of connections of ≤k were removed from the network. For the remaining network, the rich-club coefficient Φ(k) was computed as the ratio of connections present between the remaining nodes and the total number of possible connections that would be present when the set would be fully connected. Formally, the rich-club coefficient Φ(k) is given by the following (Zhou and Mondragon, 2004; Colizza et al., 2006; McAuley et al., 2007):
Weighted networks.
Following a similar principle, for weighted networks, a weighted rich-club parameter Φw(k) was computed (Opsahl et al., 2008). First, all connections of the examined network were ranked in respect by weight, resulting in a vector Wranked. Next, within M, for each value of k, the group of nodes with a degree larger than k was selected. Next, the number of links E>k between the members of the subset was counted, together with their collective weight W>k, computed as the sum of the weights of the resulting E>k connections. The weighted rich-club parameter Φw(k) was then computed as the ratio between W>k and the sum of the weights of the strongest E>k connections of the whole network, given by the top E>k number of connections of the collection of ranked connections in Wranked. Formally, Φw(k) is given by the following (Opsahl et al., 2008):
Normalization.
The rich-club effect can be quantified by the rich-club coefficient Φ(k). Since random networks—like the Erdos–Renyi model—also show an increasing function of Φ(k), due to the fact that nodes with a higher degree also have a higher probability of being interconnected by chance alone, Φ(k) is typically normalized relative to a (set of) comparable random network(s) of equal size and similar connectivity distribution, giving a normalized rich-club coefficient Φnorm (Colizza et al., 2006; McAuley et al., 2007). An increasing normalized coefficient Φnorm of >1 over a range of k reflects the existence of rich-club organization in a network. Therefore, in this study, for each examined network the rich-club curve was compared with the rich-club curve of a set of random networks, created by randomizing the connections of the network, keeping the degree distribution and sequence of the matrix intact (Rubinov and Sporns, 2010). For each network, m = 1000 random networks were computed and from each of the randomized networks, for each level of k, the rich-club coefficient Φrandom (or Φrandomw) was computed. Next, the overall Φrandom(k) was computed as the average rich-club coefficient over the m random networks. The normalized rich-club coefficient Φnorm(k) was computed as follows: For simplicity, in the remaining text, Φnorm(k) and Φnormw(k) will be referred to as Φ(k) and Φw(k), respectively.
Statistics.
To assess statistical significance of rich-club organization, permutation testing was used (Bassett and Bullmore, 2009; van den Heuvel et al., 2010). Obtained from the population of 1000 random networks (see above), the distribution of Φrandom(k) yielded a null distribution of rich-club coefficients obtained from random topologies. Next, for the range of k expressing rich-club organization, it was tested whether Φ significantly exceeded Φrandom (averaged over the examined range of k) and a (one-sided) p value was assigned to Φ as the percentage of Φrandom that exceeded Φ.
s-core decomposition
To further examine the organization of the connectome, the s-core structure of the group brain network was computed. The s-core is as a weighted equivalent of the more commonly known k-core. The k-core of an unweighted graph G is defined as the maximal connected subgraph of nodes in G in which all nodes have at least k connections and in a similar fashion, the s-core of a weighted graph Gw is the subgraph of nodes of Gw in which all connections show a summed weight of s or higher. s-core decomposition can provide insight into the hierarchical organization of a network (Chatterjee and Sinha, 2008; Hagmann et al., 2008). s-core decomposition proceeds by successively pruning the connections of the network, along the following steps.
For a particular sum of weights s:
Step 1. Remove all nodes of whose sum of weights <s, resulting in a pruned connectivity matrix M′.
Step 2. From the remaining set of nodes, compute the connectivity strength s′ for each node. If nodes are found that have a lower level of connectivity than s′, step 1 is repeated to obtain a new M′; otherwise, proceed to step 3.
Step 3. The remaining subset of nodes forms the s-core of the network.
Finally, for each node i in the network its core-level can be determined, as the maximal s-core node i is participating in. As such, a higher core-level of a node expresses a more central role of the node in the overall network.
Modularity, provincial and connector hubs
To examine the community (module) structure of the brain network, and the role of hubs in interconnecting distinct modules, module partitioning was performed (Rubinov and Sporns, 2010). Nodes identified as hubs were further classified into “provincial” and “connector” hubs, based on their level of participation in their local module and their level of connectedness to other modules. The level of “intramodule” connectivity versus “intermodule” connectivity of a node can be expressed by the “participation index” of a node (Sporns et al., 2007; Rubinov and Sporns, 2010). The participation index is formally given by the following: with Nm, the number of modules; ki, the degree of node i; and kim, the number of connections from node i to module m. From the selected hubs, connector hubs—interconnecting modules—were defined as nodes with a Pindexi > 0.5, and provincial hubs—connecting nodes within a module—were selected as nodes with a Pindexi ≤0.5, with the threshold 0.5 marking the top 15% of the node-specific Pindex values across the 82 nodes of the network (Sporns et al., 2007).
Rich-club centrality
To examine the role of the rich club in the global network structure the centrality of the rich club was measured. To this end, the percentage of shortest paths between any two non-rich-club nodes that passed through the rich club was computed (Mw-nos). Two heuristics were examined: (1) the percentage of the number of shortest paths that passed through at least one of the rich-club nodes, normalized to the number of all shortest paths in the network, and (2) the number of shortest paths between all nodes in the network that passed through at least one rich-club edge.
Rich club in targeted and random attack
Network rich clubs have been noted to play a central role in overall network structure (Xu et al., 2010) having a strong positive impact on the global efficiency of the network. The role of a node (or a set of nodes) in the level of global efficiency of a network can be evaluated by examining the damage inflicted by attack on that node, simulated as a decrease in the weights of its connections. We further examined the role of the rich club in global brain network structure by probing for the effects of network damage on global efficiency due to “attack” (illustrated in Fig. 7a). Two forms of attack can be distinguished: “targeted attack” and “random attack”. (1) Targeted attack refers to examining the effects of damage to a specific set of nodes/edges in the network. (2) In contrast, random attack refers to examining the effects of damage resulting from attacking a set of random nodes/edges in the network. In this study, we compared the effect on global efficiency of the network following targeted attack on rich-club connections versus random attack.
Targeted attack to rich-club connections.
In this condition, the weights (NOS-weighted) of the connections of the rich club were attacked at two levels, inflicting 50 and 100% damage to all weights of connections that interconnected the rich-club nodes, by setting the weight of each of the rich-club connections to the following: with τ, the level of inflicted damage, thus obtaining a damaged network Mτ. The total damage to the network was computed as follows: To assess the potential impact on global information flow we computed, for each Mτ, the level of global efficiency (Rubinov and Sporns, 2010) normalized to the level of global efficiency of M.
Random attack.
The effect on global efficiency of targeted attack on rich-club connections was compared with random damage to the network. To simulate random attack, we randomly selected a set of O non-rich-club connections of the brain network Mw-nos and damaged this set at two levels of τ (i.e., 50 and 100%).
Random attack to hub (non-rich-club) connections.
In addition to the random attack condition, in which any connection in the network could be attacked, a more restricted random condition was examined, in which damage was limited to the connections of the 12 rich-club nodes with the rest of the network (i.e., excluding connections between rich-club members). This condition serves to examine whether the effects of targeted attack could be attributed to damaging the rich-club connections and not to just to lowering the nodal degree of the hubs of the brain network, regardless of whether they formed a club or not. Similar to the target condition, damage was inflicted by reducing the weights of the connections of the selected random set by 50 or 100%, respectively, and the level of global efficiency of the resulting damaged network Mτ was computed.
For both conditions of random attack, the total weight (i.e., the sum of the weights of the connections) of the damaged connections was equal to the total damage to the rich club at level τ. For each level of damage (i.e., 50 and 100%), a sample of 10,000 random cases was examined.
High-resolution analysis
The main focus of this study is the examination of rich-club organization of the human connectome. For this, we used a regional approach, parcellating the brain network in 82 brain regions. However, it has been noted that graph metrics can be quite dependent on the resolution of the network (van den Heuvel et al., 2008b; Wang et al., 2009; Fornito et al., 2010; Zalesky et al.). In the absence of objective whole-brain parcellation strategies that identify coherent brain regions based on anatomical or functional criteria (Wig et al., 2011), some studies have advocated the use of higher-resolution networks, going up to >1000 smaller parcels, instead of using a coarse and to some extent arbitrary parcellation scheme of 82 brain regions (Hagmann et al., 2008; van den Heuvel et al., 2008b). We, therefore, performed an initial high-resolution network analysis, examining rich-club organization at high resolution. This analysis included the following steps. First, instead of parcellating the cortical surface into 68 large brain regions, for each individual dataset, the left and right cortical mantle was divided into a total of 1170 parcels (585 in each hemisphere), following a similar procedure as described by Hagmann et al. (2008) (using the Freesurfer software suite). Note that, in this parcellation scheme, only cortical regions are included. Second, for each individual dataset, this parcellation was used to create a high-resolution connectivity matrix Mhw-nos, following the exact same steps as described for the regional analysis (see above). Third, from the individual matrices, a group consensus connectivity matrix was constructed, by averaging the values of the individual high-resolution connectivity matrices. For the group matrix, pathways that were present in 50% or more of the group of subjects were taken into account. Fourth, similar to the regional approach, Φw(k), Φrandomw(k) and Φnormw(k) were computed from the resulting group network and the rich-club nodes were examined. Rich-club nodes (i.e., nodes with degree >k) were labeled with the anatomical name from the regional parcellation scheme. To determine the extent to which brain regions participated in the rich club, for each region i a participation ratio Pr was computed, calculated as the number of rich-club parcels that overlapped with region i, normalized to the number of participating parcels based on chance considering the volume of a region (i.e., the number of parcels in region i) and the total size of the rich club (i.e., number of rich-club parcels). As such, a participation ratio Pri > 1 is reflecting an above average participation of region i in the rich club.
Results
Graph metrics
Table 1 summarizes the characteristic graph measures of the group connectomes (unweighted, NOS-weighted, FA-weighted networks, high-resolution network) and of the individual networks, describing clustering, path length, global efficiency, and assortativity. As expected, results show a high level of (normalized) local clustering (γ > 1) and a high level of (normalized) global efficiency (λ ∼ 1), confirming a small-world organization of structural brain networks (Hagmann et al., 2007, 2008; van den Heuvel et al., 2010). Figure 2 summarizes the node-specific degree k, betweenness centrality, path length, and clustering coefficient values of the nodes of Mgroupw-nos. We note a high level of consistency of the normalized clustering coefficient and path length across the four different network types (unweighted and NOS/NOSroi- and FA-weighted), suggesting substantial agreement in the topological organization between the four networks.
Rich-club organization of the brain network
Figure 3 shows the normalized rich-club coefficient curves Φ(k) (Fig. 3a) and Φw(k) (Fig. 3b–d) for the group-averaged networks. Rich-club coefficients Φ(k) of the connectome are given in dark gray, rich-club coefficients of a set of comparable random networks Φrandom(k) are given in light gray, their ratio expressing the normalized rich-club coefficient Φnorm(k) is given in red.
Unweighted network
Figure 3a shows Φ(k), Φrandom(k), and Φnorm(k). For a range of k = 11 to k = 17, Φnorm(k) clearly shows characteristic rich-club organization, with Φnorm(k) > 1 and increasing over k [p = 0.028, permutation testing (see Materials and Methods)]. Examining the individual structural networks, all subjects showed comparable Φnorm(k) functions.
Weighted networks
Figure 3b shows the rich-club coefficient curves Φw-nos(k) of the NOS-weighted group averaged network, the (average) rich-club coefficient values of the set of comparable random networks Φrandw-nos(k), and their ratio, giving Φnormw(k) values. The figure shows rich-club organization of Φw-nos(k) > Φrandw-nos(k) for k = 11 to k = 21 (p < 0.001, permutation testing). Figure 3c shows the rich-club coefficients Φw-nosROI depicting rich club organization for k = 10 to k = 17 (p = 0.0050, permutation testing). In addition, Figure 3d shows the rich-club coefficients for the FA-weighted group averaged network, being Φw(k), Φrandw(k), and Φnormw(k) values, illustrating rich-club organization for k = 9 to k = 23 (p < 0.001, permutation testing). Overlapping the unweighted data, both the FA-weighted, NOS-weighted, and NOS-ROI-weighted networks displayed rich-club organization, with Φnormw > 1 for k ≥ 9.
Rich-club formation
Figure 4, c and e, shows the anatomical arrangement of the rich club in the group-averaged structural brain network, with the size of the nodes marking node degree. Figure 4c illustrates all connections of the group-averaged brain network with the rich-club connections (i.e., connections between rich-club regions) marked by thickened dark blue lines; Figure 4e shows only the connections of the rich club. In addition, Figure 4d shows all connections made by rich-club regions to other regions of the brain. Within-rich-club connections are marked in dark blue, and all other connections by light blue lines. Figure 4d illustrates that almost all brain regions are directly linked to the structural rich-club core of the brain network.
Rich-club regions
The illustrated rich club (of k = 17) comprised 12 regions, including the superior parietal, precuneus, superior frontal cortex, putamen, hippocampus, and thalamus, all in both left and right hemispheres. A high level of consistency of rich-club participation among regions identified in the group-averaged matrix was found across the group of individual networks.
s-core decomposition
Figure 5 displays the s-core decomposition of the structural brain network (NOS-weighted). Figure 5a depicts the s-core decomposition of the total brain network, showing a high core level of subcortical and superior frontal regions (red) and medial parietal regions (orange). In addition to s-core decomposition of the whole-brain network (including both subcortical and cortical regions), s-core decomposition was also performed on only cortical regions (i.e., excluding subcortical regions of the network) to compare s-core decomposition to the results of previous studies (Hagmann et al., 2008). Figure 5b depicts s-core decomposition of the neocortex, with a high core level of bilateral superior frontal, precuneus, and superior parietal regions, largely consistent with previous results on s-core structure of the neocortex (Hagmann et al., 2008). In Figure 5, a and b, the color of the nodes depicts their core level. The figure clearly indicates a high core level of the previously identified rich-club regions (i.e., bilateral precuneus, superior frontal cortex, superior parietal cortex, subcortical regions thalamus, hippocampus, and putamen and bilateral precentral regions).
Modules, provincial and connector hubs, centrality
Module partition using Newman's modularity metric (Newman, 2006) of the group averaged connectome resulted in four distinct modules. Figure 6a illustrates the modules obtained for the group averaged NOS-weighted network, showing the formation of a left anterior (red module), right anterior (yellow), left posterior (green), and right posterior (orange) subnetwork. In addition, Figure 6b illustrates the weighted backbone of the NOS-weighted brain network (k average of 5) (Rubinov and Sporns, 2010) with the provincial hubs (connecting nodes within a module) in yellow and the connector hubs (interconnecting different modules) marked in red. Figure 6b shows that all midline cortical rich-club nodes (i.e., bilateral precuneus, superior frontal, superior parietal) are connector hubs, playing an important role in between-module connectivity, while subcortical rich-club regions (bilateral thalamus, putamen) play an important role in module structure.
Betweenness centrality
Figure 6c–e depicts the level of centrality of the nodes of the network. The betweenness centrality of a node/edge in a network expresses the number of shortest paths between any two nodes in the network that pass through that particular node/edge. Figure 6c–e illustrates the level of node and edge betweenness centrality for the unweighted, NOS-weighted, and FA-weighted networks, respectively. The size of the nodes illustrates the level of betweenness centrality of the nodes, and the thickness of the edges depicts the level of edge centrality. As expected, we find a high level of centrality for rich-club nodes (in all three conditions). Interestingly, only for the NOS-weighted network, rich-club connections showed a higher level of edge centrality.
Rich-club centrality
Examining the centrality of the observed rich club, 89% of all shortest paths between all non-rich-club regions in the network passed through one or more of the rich-club nodes. For comparison, we randomly selected 1000 sets of nodes corresponding in size to the number of nodes contained in the rich club and again traced all short paths between nodes outside this set. We found that, on average, only 32% of all paths passed through the randomly selected node sets, significantly lower than the percentage of paths through the rich club (p < 0.001, permutation testing). Moreover, almost 66% of the shortest paths between any two non-rich-club nodes in the network were found to pass through one or more of the rich-club edges (and thus pass at least two rich-club regions), compared with an average of only 4% in the random condition (p < 0.001, permutation testing).
Rich-club attack versus random attack
Figure 7b shows the effects of targeted rich-club attack and random attack on the structural brain network for 50 and 100% attack, for the Mw-NOS weighted network with the rich club defined at k = 17 (overlapping the rich-club regions as illustrated in Fig. 4). The figure shows that targeted attack on rich-club connections resulted in a strong reduction in global efficiency with a maximum of 18% when all rich-club connections are destroyed (black, targeted attack; 100% damage), compared with a maximum of ∼5% in the most severe of the two random conditions (light gray, random attack; dark gray, random attack to non-rich-club hub connections) (value of p < 0.0001; 10,000 permutations). Figure 7b depicts the total level of damage to the network (i.e., the percentage of weight removed from the network) at each level of damage for the specialized rich-club attack and random attack condition.
High-resolution network
In addition to the regional analysis, an exploratory high-resolution analysis was performed. Figure 8a shows the high-resolution parcellation scheme (1170 parcels). Figure 8b shows the Φw(k), Φrandomw(k), and Φnormw(k) curves (high-resolution NOS-weighted). For a range of k = 7 to k = 15, Φnormw(k) was found to be significantly greater than 1 (p < 0.001, permutation testing of 1000 random networks), suggesting a rich-club organization of the human connectome at high spatial resolution. Figure 8c shows the group connectome. Rich-club nodes and the connections interconnecting them (i.e., the rich club) at each level of k are colored from green to yellow to red, with the orange and red nodes marking the highest level rich-club nodes and their rich-club connections. Figure 8d shows the same network, but displaying only the rich-club connections (k > 9). To define the anatomical relevance of these high-resolution nodes, rich-club parcels were overlaid with the regional parcellation scheme. Next, for each cortical region i, a rich-club participation ratio Pri was computed (see Materials and Methods), with regions showing a Pri > 1 showing an above average level of rich-club participation. In accordance with the rich-club results of the regional analysis, regions showing Pr > 1 included the left and right precuneus (1.39 and 1.78, respectively), as well as the left and right superior frontal (1.15 and 1.60, respectively). At the rich-club threshold of k > 9, the superior parietal cortex, however, revealed a Pr of ∼1, not suggesting a level of rich-club participation as shown in the regional analysis. Lowering the threshold to k > 8 did yield a Pri of 1.18 for the left superior parietal cortex, which suggests that only restricted portions of the superior parietal cortex are involved in rich-club formation. In addition to the superior frontal cortex, precuneus and, to some extent, parietal cortex, the high-resolution analysis revealed additional rich-club involvement of other frontal, parietal, and temporal brain regions (k > 9). These regions comprised posterior midline regions, including left and right posterior cingulate (Pri = 1.59 and 1.91, respectively), the isthmus of the cingulate (1.39; 1.91), the right cuneus (3.00), the lingual gyrus (1.57; 1.27), the pericalcarine cortex (3.82; 3.82), and the right paracentral cortex (1.21); anterior midline regions, including medial orbitofronal cortex (1.53; 1.47), the left and right caudal anterior cingulate (1.70; 1.91) and left rostral anterior cingulate cortex (2.29), right caudal middle frontal cortex (1.57); temporal regions, including left and right entorhinal cortex (2.29; 1.27), left fusiform gyrus (2.11), left and right parahippocampal gyrus (2.55; 3.35), and superior temporal cortex (1.74; 2.04); and, finally, frontal-temporal regions, including the left and right insular cortex (2.29; 2.12).
Discussion
The main finding of this study is the existence of rich-club organization of the human connectome. Bilateral frontoparietal regions, including the precuneus and superior frontal and parietal cortex, as well as important subcortical regions including the hippocampus, thalamus, and putamen were found to be not only individually central but also densely interconnected, together forming a rich club. Significant rich-club organization was observed in all four regional based versions of structural brain networks we examined, a regional-based unweighted network representing the existence of interregional connections, NOS-weighted/NOS-ROI-weighted networks incorporating information on the quantity of these connections, and a FA-weighted network incorporating information on the efficacy of brain connectivity (Fig. 2a–d). Furthermore, rich-club organization was also found in a high-resolution network (Fig. 8), strongly suggesting the existence of a robust densely interconnected rich club in the human connectome.
Our findings are in line with previous reports on a central hub role of neocortical precuneus and superior frontal regions (Sporns et al., 2007; Hagmann et al., 2008; van den Heuvel et al., 2010), and on the existence of a densely connected structural core in posterior medial cortex including the precuneus and superior parietal cortex (Hagmann et al., 2008). The observed rich-club organization of the human network extends earlier observations in nonhuman species reporting on the tendency of hubs in the macaque and cat cortex to form “hub complexes” of high interconnectivity density (Sporns et al., 2007; Zamora-López et al., 2009) and reinforces the results of a recent study suggesting a hierarchical assortative organization in human brain networks, a network topology in which high-degree nodes exhibit a tendency to be interconnected (Bassett et al., 2011). Our study now shows that prominent brain hubs are not only interconnected, but form a dense rich club, extending from the structural core along the cortical midline and into specific subcortical regions (Figs. 4, 5). The aggregation of hubs into a rich club suggests that the communication hubs of the brain do not operate as individual entities, but instead act as a strongly interlinked collective. One plausible explanation for the above random level of rich-club connectivity might be the tendency of the brain to provide a certain level of resilience to its core, in case of malfunction of one of its key hubs, regardless of a higher cost of the needed white matter wiring (Kaiser et al., 2007). Failure of a hub can have a severe effect on the level of global communication efficiency of a network, due to its central role in the network (Albert et al., 2000; Callaway et al., 2000; Zhou and Mondragon, 2004).
What is the relationship of the rich club to known functional systems and resting-state networks (RSNs)? One possible hypothesis is that the rich club corresponds to, or overlaps, one or more RSNs. We note that there is indeed some overlap between the rich club and the default mode network (DMN) of the brain (Greicius et al., 2003; Mason et al., 2007; Raichle and Snyder, 2007). However, the overlap is only partial, with several DMN regions (e.g., the inferior parietal lobule, medial temporal, and medial prefrontal cortex) not found among rich-club regions (Fig. 4). An alternative hypothesis suggests that the rich club may serve to link different functional modules in the brain, through partial overlap with several RSNs possibly including the default mode network, the salience network, the executive network, and primary motor, visual, and auditory networks. Indeed, our initial high-resolution analysis revealed the rich club to be a distributed set of nodes participating in the default mode network (e.g., posterior cingulate cortex/precuneus, medial orbitofrontal cortex), the salience network (e.g., insula, caudal anterior cingulate cortex), the visual (cuneus, lingual gyrus) and auditory (superior temporal cortex) networks, and to some extent the executive control network (e.g., superior frontal and left superior parietal). Further studies of rich-club organization in highly resolved structural networks are needed to determine its relation to key components of resting-state functional connectivity.
A central role of the rich club of the brain may be further illustrated by an analysis of the effects of targeted and random attack on global efficiency. Attacks that specifically targeted rich-club connections impaired global efficiency approximately three times more than attacks that inflicted randomly distributed damage to the network (Fig. 7b), suggesting a central role of rich-club regions in overall brain communication. A hypothesized central role of rich-club regions in long-distance brain communication is in agreement with recent functional connectivity studies, showing an important role of functional hubs in optimizing global brain communication efficiency for healthy cognitive brain functioning (van den Heuvel et al., 2009a). Together, our structural analysis tends to suggest that the topological central embedding of the rich club of the brain makes it a focal point for whole-brain communication.
The finding of a rich-club organization of the connectome may provide new insight into how disorders affect brain topology and functioning. Studies have suggested connectome abnormalities in a wide range of neurological and psychiatric brain disorders, for example, Alzheimer's disease, amyotrophic lateral sclerosis, Parkinson's disease, schizophrenia, and autism, with each disorder affecting brain communication efficiency and capacity in a different manner (Bosboom et al., 2009; Stam et al., 2009; Lynall et al., 2010; van den Heuvel et al., 2010; Verstraete et al., 2010; Zalesky et al., 2011). In the context of our performed study, it is tempting to speculate that diseases that affect the rich-club core of the brain, as suggested for Alzheimer's disease and schizophrenia, may have more global effects on brain communication and thereby affect multiple cognitive domains. Indeed, computational neural mass modeling of structural damage to the brain network have provided first insight into the dynamical effects of virtual brain lesions, indicating that damage to topologically central brain regions have greater effects on cortical synchronization (Honey and Sporns, 2008; Alstott et al., 2009). The present study predicts that such dynamic effects of brain lesions are particularly pronounced if lesions disturb the rich-club organization of the brain.
Some issues have to be taken into account when interpreting the results. First, the rich-club curves as shown Figure 3 tend to suggest that the range of nodal degree k over which the brain network shows rich-club organization may be to some extent constrained by a lower and upper bound. Although our analysis is not specifically designed for this, there may be room for speculation about possible biological origins of these bounds. The lower bound may perhaps reflect a crossover point at which low degree nodes are eliminated from the rich-club set, reflecting a high level of specialization of these low degree nodes. In addition, the upper bound could reflect the absence of a fully interconnected core, which may be the result of the missing interhemispheric connections between nonhomotopic brain regions, or, alternatively, may reflect a certain level of differentiation between the highest hubs of the brain network. Future studies examining the upper and lower bounds of rich-club formation of the brain network are of particular interest. Second, a common problem in DTI studies is that, when the directional information at some point along the tract is not univocal, for example due to “crossing pathways,” the fiber tracking algorithm may become hindered in correctly tracing fiber streamlines. This may in turn result in an underrepresentation of the number of connections of the connectome. Ultra-high spatial resolution DTI and diffusion spectral imaging have been suggested to reduce some of these problems (Bassett et al., 2011); however, due to the longer acquisition times of these types of scans, the use of these techniques in large studies is currently limited. Second, in this study, the brain network was mainly examined at a relative low spatial resolution, representing the human connectome as a graph of 82 segmented brain regions. As mentioned, it has been noted that “graph resolution” can have an effect on the topological properties of the brain network (Wang et al., 2009; Fornito et al., 2010; Zalesky et al., 2010; Bassett et al., 2011, Wig et al., 2011). Indeed, using a regional approach, we were not able to replicate high positive assortativity of the connectome as reported in high-resolution networks (Eguíluzet al., 2005; Hagmann et al., 2008) (high-resolution connectome; see also Table 1). Nevertheless, similar to our regional-based results, also the exploratory high-resolution analysis (Fig. 8) showed a clear rich-club organization, suggesting the robust formation of a rich club in the connectome of the brain at a mesoscopic scale.
In this study, we report on rich-club organization of the human connectome. We demonstrate that highly connected and central brain hubs have a strong tendency to be mutually interconnected, together forming a central rich club, a focal point for whole-brain communication. Future studies assessing the functional impact of structural damage to this central rich club are of fundamental importance and may allow the formulation of novel hypotheses about how brain disorders result from damage to the human connectome.
Footnotes
This work was supported by a Rudolf Magnus Young Talent grant of the Rudolf Magnus Institute of Neuroscience (M.P.v.d.H.) and by the JS McDonnell Foundation (O.S.). We thank R. C. W. Mandl for providing the DTI preprocessing software.
- Correspondence should be addressed to Martijn P. van den Heuvel, Rudolf Magnus Institute of Neuroscience, University Medical Center Utrecht, Heidelberglaan 100, P.O. Box 85500, 3508 GA Utrecht, The Netherlands. m.p.vandenheuvel{at}umcutrecht.nl