Individual variability in the anatomical distribution of nodes participating in rich club structural networks

With recent advances in computational analyses of structural neuroimaging, it is possible to comprehensively map neural connectivity, i.e., the brain connectome. The architectural organization of the connectome is believed to play an important role in several biological processes. Central to the conformation of the connectome are connectivity hubs, which are likely to be organized in accordance with the rich club phenomenon, as evidenced by graph theory analyses of neural architecture. It is yet unclear whether rich club connectivity hubs are consistently organized in the same anatomical framework across healthy adults. We constructed the brain connectome from 43 healthy adults, based on T1-weighted and diffusion tensor MRI data. Probabilistic fiber tractography was used to evaluate connectivity between each possible pair of cortical anatomical regions of interest. Connectivity hubs were identified in accordance with the rich club phenomenon applied to binarized matrices, and the variability in frequency of hub participation was assessed node-wise across all subjects. The anatomical location of nodes participating in rich club networks was fairly consistent across subjects. The most common locations for rich club nodes were identified in integrative areas, such as the cingulate and pericingulate regions, medial aspect of the occipital areas and precuneus; or else, they were found in important and specialized brain regions (such as the oribitofrontal cortex, caudate, fusiform gyrus, and hippocampus). Marked anatomical consistency exists across healthy brains in terms of nodal participation and location of rich club networks. The consistency of connections between integrative areas and specialized brain regions highlights a fundamental connectivity pattern shared among healthy brains. We propose that approaching brain connectivity with this framework of anatomical consistencies may have clinical implications for early detection of individual variability.


Introduction
Recent advances in neuroimaging now make it possible to chart the organization of neural connectivity across the entire human brain using magnetic resonance (MRI) diffusion tensor imaging (DTI) (Hagmann et al., 2008). Known as brain connectomes (Sporns et al., 2005;Sporns, 2011), the generation of whole brain maps of neural architecture is becoming increasingly popular as evidence accumulates in favor of the view that the disruption of structural connectivity is central to the neurobiology of many neurological and psychiatric illnesses (Buckner et al., 2009;Seeley et al., 2009;van den Heuvel and Sporns, 2011).
A common challenge encountered during connectome analysis is the identification of connectivity hubs, that is, nodes that are essential to the network, granting it a structured and nonrandom conformation. In general, hubs are expected to function as the basis of the brain's integrative capacity (van den Heuvel and Sporns, 2011;Sporns, 2013), and possess a high degree of connectivity, short neuronal path lengths, and high centrality (van den Heuvel and Sporns, 2011).
A promising strategy in identifying network hubs involves the assessment of nodal participation in so-called rich club subnetworks (Sporns, 2014). Assuming that central hubs of the connectome are more likely to be highly interconnected, rich-club nodes characteristically show higher connectivity with each other, beyond what would be expected by chance given their degrees (i.e., number of connections with other nodes).
Despite many reports demonstrating the presence of the rich club phenomenon in the brain (van den Heuvel and Sporns, 2011;van den Heuvel et al., 2012;van den Heuvel and Sporns, 2013;Collin et al., 2014) the actual anatomical variability of nodes pertaining to these kinds of networks across the healthy populations has yet to be established. Studies have so far demonstrated that the number and approximate locations of the hub nodes are similar (van den Heuvel and Sporns, 2011), while the exact position of each node differs slightly (Mueller et al., 2013). Thus, the normative definition of regional variability in rich-club nodes can aid the interpretation of future studies aiming to assess conformational differences in network architecture and hub sites, and their relationship with neurological and psychiatric diseases.
In this study, we aimed to measure the anatomical variability of nodal participation in structural rich club networks. We computed the whole brain connectome from a relatively large cohort of healthy individuals and we assessed the frequency with which anatomical sites were involved in the rich club. We hypothesized that, while the overall size of the rich club network (i.e., the number of nodes in the rich-club network) remained relatively stable across individuals, anatomical differences could be observed in its composition, for instance, in terms of nodal participation in rich club configuration.

Subjects
We studied 43 right-handed, healthy subjects [mean age ± standard deviation (SD) = 37.1 ± 11.7 years, 28 females] with no previous history of neurological or psychiatric disorders. The Institutional Review Board of the Medical University of South Carolina approved this study. All subjects signed an informed consent prior to their inclusion in the present study.

Tractography and Connectome Measurement
White matter maps were obtained from segmentation of T1weighted images using FreeSurfer (Martinos Center for Biomedical Imaging, Harvard-MIT, Boston USA) (Fischl et al., 2002). Gray matter regions of interest (ROI) were obtained by the segmentation of the cortex into 41 ROIs in each hemisphere in accordance with the Lausanne anatomical atlas, distributed as part of the Connectome Mapping Toolkit [http://www.cmtk. org/] (listed in Table 1). Gray matter ROIs and white matter maps were transformed onto the image space of DWIs through linear registration using FSL FLIRT (Xue et al., 1999;Jenkinson et al., 2002). Whole brain tractography was reconstructed in DWI space. Tractography was performed using the software Diffusion Toolkit (DTK) (Wang et al., 2007). The acquisition geometry and gradients were obtained from DICOM images using dcm2nii, as part of the software MRIcron (Rorden et al., 2007). The parameters for tractography were as follows: (1) maximal angle threshold of 45 • (Mori and van Zijl, 2002); (2) inclusion mask from dynamic contrast range mask based on the diffusion-weighted image (DWI), as part of default DTK parameters; (3) inclusion mask from white matter map (in DWI space). Every white matter fiber was evaluated regarding its extreme points. If the extreme point of the fiber was located in the boundary between gray and white matter region (i.e., within approximately one or fewer voxels from the gray matter ROI), this extreme was considered to be linking this ROI. If both extremes were linking different gray matter ROIs, the fiber was counted as connecting these ROIs. For each possible pair of ROIs, the number of fibers connecting the pair was computed and recorded in a connectivity matrix. The connectivity was corrected due to the biases arising from the length of streamlines and volume sizes of the different ROIs (Hagmann et al., 2007). These steps were performed through in-house scripts written in MATLAB.
For each subject, the resulting connectome was represented as a weighted connectivity matrix A, symmetrical along its main diagonal (i.e., undirected), where each entry Aij corresponded to the weighted link between ROIs i and j.

Calculation of the Rich Club Network
For each subject, the nodes participating in the rich club network were defined as those with high degree, which were more densely connected with other highly-interconnected nodes than it would be expected just by chance based on their high degree (Sporns, 2013(Sporns, , 2014. We employed the same methodological approach previously described by Sporns (Sporns, 2013) in order to identify such rich-club nodes. We binarized each individual network based on a range of fixed density thresholds, whereby all links with weight higher than 60th to the 95th percentile of individual weight distribution were maintained. This step maintained only the reproducible core of the network, which was associated with high number of streamlines in each subject.
For each fixed density threshold, the rich club coefficient was iteratively calculated across all degrees (k) observed in the network, as the sum of connections of the sub-network composed of nodes with higher degrees (>k), divided by the number of all possible connections within this sub-network. For each network degree, the rich club coefficient is denoted as , and a corresponding vector is obtained for each subject representing the coefficient across all degrees, (k). A normalized coefficient was then obtained by calculating the element by element ratio between (k) and random (k), where random (k) denotes the rich club coefficient calculated (also across the range k) for 1000 random networks obtained by shuffling the links in the individual network, preserving the degree distribution compared with the individual network. For each subject, nodes participating in the rich club network were defined as those with degree k such that (k)/ random (k) > 1. We examined the variability across subjects of nodal participation in the rich club network by evaluating the relative frequency with which each node was part of the rich club network. We compared regional participation in the rich club networks across all fixed density thresholds through a series of chi-squares. Nonetheless, we focused our analyses on the data obtained from the networks including only links within the highest 5th percentile of link weight, which is a sparsity level more likely to represent the core framework of structural organization of the networks obtained with probabilistic tractography (Hinne et al., 2015).
Finally, we also evaluated whether there was a correlation between the frequency of participation in the rich club network and age.

Results
There were no statistically significant differences in regional rich club participation when comparing different fixed density thresholds. The relative participation of nodes in the rich club network across all thresholds can be observed in Supplementary Figure 1.
Focusing on networks composed of links with weight greater than the 95th percentile (in order to better represent network cores) we observed a relatively left vs. right symmetrical distribution of nodal degrees. As demonstrated in Figure 1, there was not a statistically significant difference in node-wise distribution of degrees across hemispheres, except for the nucleus accumbens (ROI 39), which was associated with a higher degree on the right hemisphere (t = 3.98, p = 0.0003).
The hemispheric distribution of node-wise rich club participation was also fairly symmetrical across hemispheres. Nonetheless, the lateral occipital region was more commonly part of the rich club network in the right (t = 2.95, p = 0.005). Notably, employing a more liberal statistical threshold, (i.e., 0.01 < p < 0.05), we observed a trend toward higher participation of the right pericalcarine region (t = 2.07, p = 0.04), the right nucleus accumbens (t = 2.67, p = 0.01), and the left transverse temporal region (t = −2.35, p = 0.02).
Group wise average curves demonstrating (k) (i.e., the rich club coefficient per degree k), random (k), and their ratios are shown in Figure 2. As described in the methods, the criterion for defining nodes participating in the rich club network is defined as (k)/ random (k) > 1. The nodal degrees where the ratio (k)/ random (k) was greater than 1 was variable among subjects, as demonstrated in Supplementary Figure 2, which demonstrates (k), random (k), and (k)/ random (k) for all subjects.
The average number of nodes pertaining to the rich club, per subject, was 45 ± 13 nodes. A chi-square comparison of the observed versus expected frequencies for each node FIGURE 1 | Nodal degree (mean ± SD) for each ROI in the left (blue) and right (red) hemispheres. A significant difference between the right and left hemispheres was found exclusively for ROI 39 (nucleus accumbens), marked with an asterisk. Error bars are SD.
FIGURE 2 | Rich club coefficient as a function of nodal degree (x-axis) for the average connectome of participants in this study ( , blue) and for 1000 random networks with similar degree distribution ( random , green). The left Y-axis demonstrates the rich club coefficient. The shaded areas represent the interval within 1 SD above and below the mean, which is demonstrated by the continuous colored line. The red line represents the proportion between and random as a function of degree (X-axis), and the ratio is demonstrated on the right Y-axis.  The node-wise percentage of participation in the rich club for each node is demonstrated anatomically in Figure 3 and in Supplementary Figure 3. The individual variability can also be appreciated in Supplementary Figure 4, which demonstrates the nodes participating in the rich club network for each individual assessed in this study.
We did not observe a relationship between age and the number of nodes pertaining to the rich club network. However, we observed a trend toward a higher participation of the left entorhinal cortex in the rich club network among older individuals (i.e. older than the median age in our sample, 35 years), Chi 2 = 9.382, p < 0.005.

Discussion
In this study, we demonstrated that nodes participating in rich club networks are consistently located in relatively equivalent anatomical locations across healthy adults. We also observed that the rich club nodes are located on all brain lobes in a fairly symmetrical distribution either in their medial or lateral aspects.
Interestingly, the anatomical areas more commonly associated with rich club nodes were those located in the proximity of medial integrative areas, such as the cingulate and pericingulate regions, medial aspect of the occipital areas and precuneus; or else, they were found in specialized brain regions (such as the oribitofrontal cortex, caudate, fusiform gyrus and hippocampus). These results suggest that these areas not only play an important functional role, but they also play an important part in the configuration of the structural brain networks in healthy individuals (Sporns, 2013). Specifically, they may serve as structural hubs and provide the structural framework for integration across different functional domains (Sporns, 2014). Our observations are in accordance with the theory that informational Frontiers in Neural Circuits | www.frontiersin.org  processing associated with a non-randomly and consistently organized structural pattern (van den Heuvel et al., 2012). We observed mild hemispheric asymmetry in regional rich club participation due to a higher prevalence of the right nucleus accumbens. This finding may be related to functional specialization of the nucleus accumbens in motor (Budilin et al., 2008), behavioral control (Savjani et al., 2014) and/or learning and memory (Downar et al., 2011). Furthermore, the higher participation of the entorhinal cortex in the rich club networks related to aging may be signifying increased medial temporal lobe network centrality in relationship with compensatory aging memory processing (Du et al., 2003).
In this study, we employed an atlas featuring 82 ROIs, since this parcellation method is commonly used in the neuroimaging literature (Hagmann et al., 2008) and our results can be related to other studies employing a similar approach. Nonetheless, it is yet unknown whether cortical parcellation methods exert a significant influence on the ensuing network architecture. In theory, a higher number of nodes could lead to a better evaluation of finer grained network complexity. With current limitations in spatial resolution of DTI, building the connectome with smaller ROIs can lead to an inability to resolve micro-connectivity.
Furthermore, our results are related to structural networks constructed using probabilistic tractography and should be interpreted in this context. Probabilistic tractography is a commonly used approach to evaluate DTI data, with the ability to partially resolve complex fiber anatomy and fiber crossings (Behrens et al., 2007). The streamline count obtained from probabilistic tractography is considered to be a biophysical representation of the underlying axonal bundles and myelin sheath (Jones et al., 2013). However, since the calculation of tractography streamlines is assessed from the distributions on voxel-wise principal diffusion directions, there is an inherent possibility that noise or low probability connections are included in the connectome. For this reason, the evaluation of networks based on density thresholds permits a more accurate assessment of the core framework of connectivity. Thus, we used the number of streamlines as a binary threshold method-i.e., regional links with a very high number of fibers were likely to be considered a core in the connectome framework, and therefore maintained in the network analyses. This "fixed density threshold" approach does not use number of fibers as weighted links in the network in the subsequent rich club analyses. We focused on a 95th percentile threshold to maintain consistency with the methodology from previous studies (van den Heuvel and Sporns, 2011;Cammoun et al., 2012;Hinne et al., 2015). Nonetheless, our exploratory analyses using different fixed density thresholds demonstrate a fair reproducibility in regional rich club participation across multiple sparsity levels.
There are many potential and practical implications from observations from this study, especially with regard to gain or loss of hub function by brain regions in specific pathological or developmental processes. For example, while studies to date have evaluated the number of rich club networks in particular disease states (e.g., see van den Heuvel and Sporns, 2011), it is conceivable that positioning of hubs may also be affected by network rearrangement as a consequence of neuropathological or neurodevelopmental processes. This may also aid in the further characterization of phenotypical manifestations of neurological diseases, as, for example, different network rearrangements may be associated with neurological or psychiatric symptoms. Thus, by better evaluating the plastic changes associated with neural architecture organization, it may be possible to correctly identify biomarkers related to neurological and neuropsychiatric problems. As an example, different neural networks may be involved in the generation and maintenance of seizures in patients with temporal lobe epilepsy (Bonilha et al., 2012). However, the clinical and electrophysiological manifestations of epilepsy may be indistinguishable regarding the neuronal networks that originate seizures. This is clinically relevant, as response to treatment of epilepsy may be directly associated with the underlying neurobiological mechanisms of seizure generation and propagation (Bonilha et al., 2013). By evaluating changes in core configuration of temporal lobe network hubs, it may be possible to better identify the conformation of different phenotypes of temporal lobe epilepsy, and to potentially contribute with similar approaches to other pathological models.
Another example of direct application of anatomical positioning of rich club nodes is the evaluation of loss and reestablishment of hubs in relationship with brain injury. For example, it is commonly observed that subjects with language processing problems (aphasia) after cortical and subcortical injury after strokes may recover spontaneously or with speech therapy (Holland et al., 1996;Robey, 1998). It is commonly assumed that recovery involves the re-establishment of new networks to support the regaining of function. Nonetheless, it is still undefined whether a greater loss of network hubs may actually prevent a more full recovery. Similarly, it is unclear whether recovery may be associated with the reclaiming of hub status by areas that are typically not associated with the rich club.
In this study, we aimed to provide a description of the anatomical positioning and variability of nodes participating in rich club networks. The findings from this study may provide a reference and framework for future studies that assess changes in this pattern in association with neurobiological mechanisms related to health and disease.