## Abstract

The network architecture of functional connectivity within the human brain connectome is poorly understood at the voxel level. Here, using resting state functional magnetic resonance imaging data from 1003 healthy adults, we investigate a broad array of network centrality measures to provide novel insights into connectivity within the whole-brain functional network (i.e., the functional connectome). We first assemble and visualize the voxel-wise (4 mm) functional connectome as a functional network. We then demonstrate that each centrality measure captures different aspects of connectivity, highlighting the importance of considering both global and local connectivity properties of the functional connectome. Beyond “detecting functional hubs,” we treat centrality as measures of functional connectivity within the brain connectome and demonstrate their reliability and phenotypic correlates (i.e., age and sex). Specifically, our analyses reveal age-related decreases in degree centrality, but not eigenvector centrality, within precuneus and posterior cingulate regions. This implies that while local or (direct) connectivity decreases with age, connections with hub-like regions within the brain remain stable with age at a global level. In sum, these findings demonstrate the nonredundancy of various centrality measures and raise questions regarding their underlying physiological mechanisms that may be relevant to the study of neurodegenerative and psychiatric disorders.

## Introduction

The emergence of resting-state functional magnetic resonance imaging (R-fMRI) has made it feasible to attain a high-resolution map of the human functional connectome (Biswal et al. 1995; Sporns et al. 2005; Fox and Raichle 2007; Biswal et al. 2010). However, with new advances come new challenges (DeFelipe 2010; Friston and Dolan 2010). In particular, the brain’s functional connectome is necessarily dynamic as it underpins a multitude of brain states involving emotion, cognition, action, perception, and sensation (Carhart-Harris and Friston 2010; Friston 2010; Deco et al. 2011). To understand the origin of these complex functions in terms of the brain’s functional connectivity, a first step consists of creating maps of potential connections in the brain and characterizing the complexity of their functional interactions.

Two fundamental approaches have dominated the functional connectivity literature to date: seed-based correlation analysis and independent component analysis (Cole, Smith, et al. 2010). Seed-based approaches estimate the strength and significance of pairwise relationships between regions of interest and all other voxels in the brain. In contrast, independent component analysis attempts to identify sets of brain regions that are separable on the basis of statistical patterns in their dynamic time series. Both approaches have enjoyed significant success in leading to the discovery and detailed mapping of large-scale brain intrinsic connectivity networks (ICNs) (Beckmann et al. 2005; Damoiseaux et al. 2006; Margulies et al. 2007, 2009). However, neither of these approaches fully characterize the brain’s functional connectome. While the seed-based method can be used to map the functional connectome in its entirety, it only provides a series of relationships between any given region and all other regions without taking into account the full pattern of connections. Likewise, independent component analysis provides information about how regions may be related within large-scale networks but does not reveal patterns of connectivity between these networks.

Recently, researchers have begun to apply graph theory–based network analysis to explore brain connectivity (Bullmore and Sporns 2009; Rubinov and Sporns 2010; Wang et al. 2010). This approach can be used to characterize functional connectivity within the whole-brain network (i.e., the functional connectome). In particular, centrality, a class of graph theory–based network measures assessing the centrality or functional importance, is receiving substantial attention (Koschützki et al. 2005). Unlike seed-based approaches or independent component analysis, centrality measures take into account a given region’s relationship with the entire functional connectome and not just its relation to individual regions or to separate larger components. As such, centrality measures allow us to capture the complexity of the functional connectome as a whole.

Initial examinations of centrality for the functional connectome have primarily focused on the identification of “functional hubs” (Achard et al. 2006; Hagmann et al. 2008; Buckner et al. 2009; He et al. 2009; Joyce et al. 2010; Lohmann et al. 2010; Tomasi and Volkow 2010; Fransson et al. 2011; Tomasi and Volkow 2011a). While some concordance of findings has been noted across studies (e.g., high centrality in posterior cingulate, ventral precuneus, medial prefrontal/parietal cortex, insula, temporal cortex, and parahippocampal gyrus), there are also some notable variations in the specific identification of regions as hubs. In part, such variations may reflect differences in the specific measures of centrality or the large-scale parcellation templates employed (Wang et al. 2009). While an increasing number of studies have begun to employ specific centrality measures (Joyce et al. 2010; Lohmann et al. 2010), mainly focused on detection of functional hubs, we lack a more concrete understanding of how multiple centrality measures are related to different aspects of connectivity within the functional connectome.

Accordingly, to gain a more comprehensive view of the nature of connectivity in the brain independent of particular parcellation templates, we examined multiple centrality measures of the voxel-wise functional connectome in terms of their consistency across participants (i.e., voxel-wise network centrality maps [VNCM], at group level), their test–retest (TRT) reliability across time, and their phenotypic relationships. Specifically, we carried out our analyses using a large sample (i.e., 1003 R-fMRI data sets from the 1000 Functional Connectomes Project: FCP, http://fcon_1000.projects.nitrc.org) to build up VNCM. We used the New York University (NYU) TRT data set (http://www.nitrc.org/projects/nyu_trt) to assess the TRT reliability of the VNCM. Finally, given age- and sex-related differences in the properties of regions commonly cited as functional hubs (Fair et al. 2008; Biswal et al. 2010; Power et al. 2010; Zuo, Kelly, Di Martino, et al. 2010; Tomasi and Volkow 2011b), we again used the FCP data set to examine age- and sex-related variations in VNCM.

## Materials and Methods

### Participants and Brain Imaging Procedures

#### NYU TRT Data

A total of 25 individuals participated in this TRT experiment (mean age: 30.7 ± 8.8, 9 males). Three resting-state scans were obtained for each participant using a Siemens Allegra 3.0 T scanner. Each scan consisted of 197 contiguous echo planar imaging (EPI) functional volumes. Scans 2 and 3 were conducted in a single scan session, 45 min apart, and were acquired 5–16 months (mean 11 ± 4) after scan 1. More details can be found in our previous studies employing this data set (Shehzad et al. 2009; Zuo, Di Martino, et al. 2010; Zuo, Kelly, Adelstein, et al. 2010; Zuo, Kelly, Di Martino, et al. 2010). None of the participants had a history of psychiatric or neurological illness as confirmed by clinical assessment. Informed consent was obtained prior to participation. Data were collected according to protocols approved by the institutional review boards of NYU and the NYU School of Medicine. These TRT data sets are fully available for downloading from the website at NITRC (http://www.nitrc.org/projects/nyu_trt).

#### 1000 FCP Data

A total of 1003 participants (mean age: 28.1 ± 12.7 years; 434 males) from 21 centers, included in the 1000 FCP, completed a R-fMRI scan. In contrast to the original FCP publication (Biswal et al. 2010), we only included 21 centers in this study to ensure uniform brain-wide coverage. Table 1 provides information on the 21 centers. Each center’s respective ethics committee approved the submission of de-identified data obtained with written informed consent from each participant. Details of each center’s R-fMRI data sets can be found on the FCP website at http://www.nitrc.org/projects/fcon_1000.

Table 1

The “1000 Functional Connectomes” project data for current study

 Center PI N 1. Baltimore, MD, USA James J. Pekar/Stewart H. Mostofsky 23 2. Bangor, UK Stan Colcombe 20 3. Berlin, Germany Daniel Margulies 26 4. Beijing, China Yu-Feng Zang 192 5. Cambridge, MA, USA Randy L. Buckner 198 6. Cleveland, OH, USA Mark J. Lowe 31 7. Dallas, TX, USA Bart Rypma 24 8. Leiden, Netherlands Serge A R B. Rombouts 31 9. Leipzig, Germany Arno Villringer 37 10. Milwaukee, WI, USA Shi-Jiang Li 46 11. Montreal, Canada Alan C. Evans 51 12. Munich, Germany Christian Sorg/Valentin Riedl 15 13. New York City, NY, USA Michael Milham/F. Xavier Castellanos 59 14. New York City, NY, USA Michael Milham/F. Xavier Castellanos 20 15. Newark, NJ, USA Bharat B. Biswal 19 16. Orangeburg, NY, USA Matthew J. Hoptman 20 17. Oulu, Finland Vesa J. Kiviniemi/Juha Veijola 103 18. Oxford, UK Steve M. Smith/Clare Mackay 22 19. Palo Alto, CA, USA Michael Greicius 17 20. Queensland, Australia Katie McMahon 18 21. Saint Louis, MO, USA Bradley L. Schlaggar/Steven E. Petersen 31
 Center PI N 1. Baltimore, MD, USA James J. Pekar/Stewart H. Mostofsky 23 2. Bangor, UK Stan Colcombe 20 3. Berlin, Germany Daniel Margulies 26 4. Beijing, China Yu-Feng Zang 192 5. Cambridge, MA, USA Randy L. Buckner 198 6. Cleveland, OH, USA Mark J. Lowe 31 7. Dallas, TX, USA Bart Rypma 24 8. Leiden, Netherlands Serge A R B. Rombouts 31 9. Leipzig, Germany Arno Villringer 37 10. Milwaukee, WI, USA Shi-Jiang Li 46 11. Montreal, Canada Alan C. Evans 51 12. Munich, Germany Christian Sorg/Valentin Riedl 15 13. New York City, NY, USA Michael Milham/F. Xavier Castellanos 59 14. New York City, NY, USA Michael Milham/F. Xavier Castellanos 20 15. Newark, NJ, USA Bharat B. Biswal 19 16. Orangeburg, NY, USA Matthew J. Hoptman 20 17. Oulu, Finland Vesa J. Kiviniemi/Juha Veijola 103 18. Oxford, UK Steve M. Smith/Clare Mackay 22 19. Palo Alto, CA, USA Michael Greicius 17 20. Queensland, Australia Katie McMahon 18 21. Saint Louis, MO, USA Bradley L. Schlaggar/Steven E. Petersen 31

### MRI Preprocessing

#### 1000 FCP Data

A standard data preprocessing strategy was carried out using both FSL (http://www.fmrib.ox.ac.uk/fsl/) and AFNI (http://afni.nimh.nih.gov/afni). Scripts containing the processing commands employed here have been released as part of the FCP. Specifically, it comprised: 1) discarding the first 5 EPI volumes from each resting-state scan to allow for signal equilibration, 2) 3D motion correction, 3) time series despiking, 4) 4D mean-based intensity normalization, 5) band-pass temporal filtering (0.01–0.1 Hz), 6) removing linear and quadratic trends, 7) estimating a 12 parameters affine linear transformation from individual functional space to Montreal Neurological Institute (MNI) 152 space, and 8) regressing out 9 nuisance signals (global mean, white matter, CSF signals, and 6 motion parameters). The output of these preprocessing steps was a 4D residual functional volume in native functional space, for each participant. The 4D native data were registered to the MNI152 space with 4 mm resolution based on the affine transformation.

#### NYU TRT Data

The preprocessing steps were similar to those described above except for skipping step 1 (discarding the first 5 EPI volumes) due to the scanner’s default setting of disregarding the first 4 volumes. Slice-timing correction for interleaved acquisitions using Sinc interpolation with a Hanning windowing kernel was additionally added prior to 3D motion correction (step 2). Step 7 was refined by a B-spline basis nonlinear transformation from an individual functional space to MNI152 standard brain space as implemented in FSL FNIRT (http://www.fmrib.ox.ac.uk/fsl/fnirt/index.html). Finally, the residuals were registered to the MNI152 standard space with 4 mm resolution. There were 2 considerations taken into account when choosing 4 mm as our voxel size: 1) the average voxel size of the EPI raw data (i.e., the maximum among x/y/z resolutions) across the 21 sites is approximately 4 mm (3.9497 mm), thereby deceasing the potential utility of using a smaller voxel size and 2) the graph has 42 835 nodes at 3 mm voxel size, which greatly increases the requirement of both physical memory (e.g., ∼50 GB for subgraph centrality [SC]) and computational time for centrality estimates.

For both 1000 FCP data and NYU TRT data, all preprocessing did not include spatial smoothing because of 3 considerations (van den Heuvel, Stam, et al. 2008): 1) the spatial registration of 4D residuals involves the interpolation of EPI voxels (i.e., somewhat like spatial smoothing), 2) spatial smoothing could introduce artifactual local correlations between voxels unrelated to their functional connectivity, and 3) spatial smoothing can blur the boundaries and artificially increase the correlations between 2 regions with respect to their time series.

### Graph Formation of ICNs

Prior to any centrality or other graph theoretical analyses, a graph has to be abstracted and formatted from the real data. This step is crucial before conducting centrality analyses (Smith et al. 2011). We constructed graphs for ICNs at voxel level. The time series for each voxel were extracted from the preprocessed R-fMRI data to calculate a correlation matrix $R=(rij), 1≤i,j≤N$ (N is the number of voxels), where $rij$ is the temporal Pearson’s correlation of time series between the i- and j-th voxels and measures the similarity of R-fMRI signals between the 2 voxels. To build a fully connected graph, a threshold of correlation $r0$ was estimated with statistical significance P = 0.0001 (uncorrected). It was used to threshold the correlation matrix into an adjacency matrix $A=(aij)1≤i≤N, 1≤j≤N$ of a graph $G$:

(1a)
$aij={0,rij

(1b)
$aij={0,rij

A binary graph was constructed as in equation (1a); its weighted version is equation (1b). To exclude artifactual correlations from non-gray matter voxels, we restricted our voxel-wise centrality analyses to a predefined gray matter mask with gray matter tissue probability greater than 20%. The gray matter tissue probability template was released as part of tissue priors in SPM8 (http://www.fil.ion.ucl.ac.uk/spm/software/spm8).

To explore other less conservative ways of constructing the graphs, we also used correlation thresholds with a significance criterion of P = 0.005 for voxel wise analysis. The findings obtained from these graphs were highly similar and thus are not presented.

### Centrality Indices

Given a graph $G=(V,E), V=(vi)1≤i≤N$ defines all nodes and $A=(aij)1≤i≤N, 1≤j≤N$ represents the graph’s adjacency matrix, in which all nonzero elements define graph $G$’s edges set $E$. The element $aij$ of the adjacency matrix represents the connection or edge from node $vi$ to node $vj$, which is 0 if no edge exists and nonzero for an edge with a weight $aij$. In the current work, we only focus on undirected graphs, that is, the adjacency matrix$A$ is symmetric. A graph is binary if only 0 and 1 appear within its adjacency matrix; otherwise, it is a weighted graph. To quantify nodal contribution (i.e., nodal centrality) for a graph, we examined 4 different centrality measures: degree, eigenvector, page-rank, and subgraph. They can be categorized into local (degree), mesoscale (subgraph), and global centralities (eigenvector and page-rank). Each of them is a function of nodes $f(i)=f(vi), 1≤i≤N$.

#### Degree Centrality

For a binary graph, degree centrality (DC) is the number of edges connecting to a node. For a weighted graph, it is defined as the sum of weights from edges connecting to a node (also sometimes referred to as the node strength). According to the adjacency matrix of a graph, DC can be computed as in equation (3). It represents the most local and directly quantifiable centrality measure. This measure has been widely used to examine node characteristics of ICNs (Buckner et al. 2009; Bullmore and Sporns 2009; He et al. 2009; Cole, Pathak, et al. 2010; Wang et al. 2010; Fransson et al. 2011).

(3)
$DC(i)=∑j=1Naij.$

#### Subgraph Centrality

SC characterizes the participation of a node in all subgraphs comprised in a network (Estrada and Rodriguez-Velazquez 2005). The “mesoscale” nature of this index is reflected in that smaller subgraphs receive more weight than larger subgraphs. SC is typically calculated based on the spectral representation of a graph, that is, the eigenvalues and eigenvectors of the adjacency matrix of the graph,

(4)
$SC(i)=∑j=1N[μj(i)]2sinh(λj).$

In equation (4), $μj(i)$ is the i-th component of the j-th eigenvector of the adjacency matrix A, and $λj$ is the corresponding j-th eigenvalue, that is, $Aμj=λjμj$. SC can be considered as a mesoscale property of order in networks, which characterizes the odd-cyclic subgraph or closed walk centrality of the i-th node in brain networks. Motif and clustering coefficient are 2 network metrics based on subgraph abundances that are widely used to characterize a complex network (Sporns 2006; Bullmore and Sporns 2009) (e.g., its small-worldness). Although the mesoscale properties of ICNs are likely to be of interest (Sporns et al. 2005), SC has not been widely employed in investigating ICNs. SC takes into account all possible subgraphs. For a weighted graph, the adjacency matrix A is converted to $D−1/2AD−1/2$ using the degree diagonal matrix D before computing SC (Crofts and Higham 2009). To calculate SC at the voxel level, we use only the first 20 largest eigenvalues (i.e., j = 1 … 20 in eq. 4) due to the computational limitation on solving the eigensystem of the large sparse adjacency matrix. The MATLAB function “eigs” was employed to compute eigenvalues/eigenvectors for the large-scale sparse matrix, which implements deflation techniques for an implicitly restarted Arnoldi iteration (Lehoucq and Sorensen 1996). We also modify (eq. 4) to its large network version because the weight of closed walks can produce an infinite SC (i.e., beyond the machine accuracy):

$SC(i)=∑j=120[μj(i)]2(N−1N−1−λj).$

This general SC, named resolvent centrality, was recently proposed for revealing network properties in terms of matrix functions (Estrada and Higham 2010).

#### Eigenvector Centrality

Introduced by Bonacich (Bonacich 1972), eigenvector centrality (EC) is simply the first eigenvector of the adjacency matrix, which corresponds to the largest eigenvalue $λ1$ (called the principal eigenvalue):

(5)
$EC(i)=μ1(i)=1λ1Aμ1=1λ1∑j=1Naijμ1(j).$

Clearly, EC not only depends on the degree of all adjacent nodes but also on their EC. Because of this recursive property, EC is able to capture an aspect of centrality that extends to global features of the graph. Of note, similar to SC, its implementation was also based on eigs in MATLAB. Neuroscientists have started to employ EC for measuring functional connectivity in both resting and task states of brain activity (Lohmann et al. 2010).

#### Page-rank centrality

The well-known Google page-rank centrality (PC) algorithm is a variant of EC. It introduces a small probability (1 − d = 0.15) of random damping (i.e., damping factor) to handle walking traps on a graph (Boldi et al. 2009):

(6)
$PC(i)=r(i)=1−d+d∑j=1Naijr(j)DC(j).$

PC has been shown to be similar to DC for undirected graphs (Fortunato et al. 2008). However, no study has yet directly examined PC of ICNs (i.e., R-fMRI networks) in a neurobiological context. Both EC and PC are solvable for binary and weighted graphs. To accelerate the computation of PC at the voxel level, we employed an inner–outer iteration strategy (Gleich et al. 2010).

### Individual-level VNCM Analyses

For each of 21 sites, a group mask was first generated to include all voxels presented across all subjects in the site. Specifically, an individual mask included all voxels with positive standard deviation (SD) of the residual time series. All individual masks were multiplied to produce the site’s mask. All 21 sites’ group masks were then combined into a final group mask including voxels present across all 1003 subjects from the 21 sites. To optimize comparability of network size across the 21 sites, all individual network analyses were constrained within the final group mask at voxel level. Specifically, only gray matter voxels (i.e., >20% gray matter tissue probability) within the final mask were employed to calculate individual centralities.

Voxel-based graphs were generated for each subject. Based on the graph, 4 centrality measures were calculated at the individual level. These individual centrality voxel-wise maps were standardized to z-scores to make them comparable across subjects, sites, and centralities (Buckner et al. 2009; Zuo, Di Martino, et al. 2010). The z-score standardization is:

$zi=xi−μσ, 1≤i≤N.$
where μ and σ are mean and SD of the centrality measure x across all N nodes (i.e., voxels in the group mask).

### Group-level VNCM Analyses

All individual centrality z-score maps were first spatially smoothed with a Gaussian smoothing kernel (full-width at half-maximum = 6 mm) and used as inputs of subsequent group-level analyses.

#### 1000 FCP Data: One-Way Analysis of Variance

As in our previous study (Biswal et al. 2010), to examine the factors of interest simultaneously, we developed a unified general linear model frame for center-level statistical analyses. The unified statistical model is a one-way analysis of variance (ANOVA) treating centers as the between-factor. F-contrasts were used to measure the effect of centers. Overall group mean contrasts across all centers were also modeled. Specifically, we used a 1-factor 21-level ANOVA (factor: center; 1003 participants) with age and sex as covariates to examine the effects of age, sex, and center on the centrality measures. In considering the potential sensitivity of age-related differences in centrality to head motion and registration quality, the root mean squares of overall head displacement/rotation, their temporal derivative, and the registration error were included to adjust for head motion and registration. The computational details of head motion and registration error can be found in our previous study (Zuo, Kelly, Di Martino, et al. 2010). Multiple comparisons were corrected at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected).

#### NYU TRT Data: TRT Reliability Analyses

To assess the TRT reliability of each centrality measure, we used intraclass correlation coefficients (ICCs) as in our previous work (Zuo, Kelly, Adelstein, et al. 2010):

$rTRT=MSp−MSeMSp+MSe.$

In this equation, MSp is the between-participants mean square and MSe is the error mean square. We examined both intrasession (i.e., short term) and intersession (i.e., long term) ICC values based upon centralities from scans 2/3 and scans 1/2, respectively. To estimate significantly reliable centralities, ICC values were converted to Z statistics:

$Z=N−3×tanh−1(rTRT),$
where N is the number of participants. A minimum Z = 2.3 (P < 0.01) was applied to threshold Z-statistical values. Voxel-wise Z maps were corrected for multiple comparisons at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected).

#### Human Brain Functional Connectome

To visualize the full brain functional connectome and further examine relationship between different centralities, for all 21 sites, we constructed a mean functional brain network. Specifically, all correlation matrices from participants were first converted into Fisher z-values and then averaged across participants. The averaged Fisher z-matrices were back-transformed to average correlation matrices by using inverse Fisher z-transformation. Based on the averaged correlation matrix, the same graph or network construction approach as described above (see “Graph Formation of ICNs”) was used to generate the mean functional network. The mean of the correlation thresholds across all 1003 subjects was used to threshold the mean correlation matrix. Specifically, corresponding to the same P value = 0.0001, the mean of the correlation thresholds was 0.2982 for voxel-derived graphs.

We assessed 1) the visualization of the full brain functional connectome based on the ability of MATLAB (http://www.mathworks.com) to handle large sparse matrices and Gephi (http://gephi.org) to explore a large network; 2) distributions of centrality measures, and 3) correlations between centrality measures.

Finally, for each voxel, to characterize the concordance of centrality measures across 1003 subjects, we proposed another summary statistic—Kendall’s W (Legendre 2005):

(7)
$W=12Sm2(n3−n),E=−Wln(W).$

Entropy E was used to measure the richness of information encoded across centrality measures. Here, m is the number of centralities, n is the number of subjects, and S is the squared SD of ranks.

## Results

In the present work, centrality analyses were performed for both binary and weighted graphs representing the functional connectome. We primarily report findings based on binarized networks except for cases where weighted centrality analyses revealed significant different findings.

### The Human Brain Functional Connectome: Communities and Network Layout

Given the challenge of visualizing the voxel-wise connectivity matrix graph (i.e., for the entire brain functional connectome, 59 960 050 connections were detected on the mean connectivity matrix, representing a sparsity of 1.2%), we first depict the sparse adjacency matrix for the entire functional connectome (Supplementary Fig. S1A) and its reordered version (Supplementary Fig. S1B). To build a network layout of the entire functional connectome, we employed a fast force-directed layout algorithm—OpenOrd (Martin et al. 2011)—to emphasize different clusters or hierarchy in the connectome, resulting in 20 detected communities (Fig. 1A; see Table 2 for details of brain regions). We applied a repeated iterative algorithm for fast unfolding communities in large networks (Blondel et al. 2008). Specifically, each iteration includes 2 phases: one where modularity is optimized by allowing only local changes of communities; one where the communities found are aggregated in order to build a new network of communities. Such phases are repeated iteratively until no increase of modularity is possible. The number of final communities was thus found in such automatically optimal way. The voxel-wise whole-brain network layout is illustrated in Figure 1B (see its full resolution from http://lfcd.psych.ac.cn/pdfs/kfc_layout_22knodes.pdf). All visualization procedures were completed with a graphics workstation (Lenovo ThinkStation D20 armed with 16 CPU cores, 16 GB physical memory, and Nivida Quadro FX4800 video card) in the Neuropsychology and Applied Cognitive Neuroscience laboratory at Institute of Psychology, Chinese Academy of Sciences. To provide a public resource for exploring further properties of large brain networks, we uploaded the adjacency matrix of the entire functional connectome (N = 1003) to the LFCD website (http://lfcd.psych.ac.cn/vncm/sp_fullbrain_4mm_1003subs_p0001.mat).

Table 2

The functional parcellation including 20 functional communities

 Functional community 1. Hippocampus (B)—Brain-stem (B)—Cerebellum Anterior Loble (B) 2. Amygdala (B)—Parahippocampus (B)—Posterior Temporal Gyrus (R) 3. Orbital Frontal Cortex (B) 4. Striatum (B)—Medial Prefrontal Cortex (B) 5. Anterior Temporal Gyrus (R) 6. Temporal Gyrus (L) 7 Visual Cortex (B) 8. Thalamus (B) 9. Posterior Cingulate Cortex (B)—Precuneus (B)—Angular Cortex (B)—Superior Frontal Gyrus (B)—Parahippocampal Gyrus (B)—Medial Prefrontal Cortex (B)—Middle Temporal Gyrus (B) 10. Dorsal Lateral Prefrontal Cortex (B-L) 11. Middle Frontal Gyrus (R)—IFG (R) 12. Middle Frontal Gyrus (L)—IFG (L) 13. Supplementary Motor Cortex (B)—Middle Cingulate Gyrus (B) 14. Precentral Gyrus (B)—Postcentral Gyrus (B) 15. Anterior Cingulate Cortex (B)—Supramarginal Gyrus (B)—Inferior Parietal Lobule (B)—Superior Temporal Gyrus (B)—Insular Lobe (B) 16. Superior Frontal Gyrus (L)—Middle Frontal Cortex (L) 17. Middle Frontal Gyrus (R)—Inferior Parietal Lobule (B) 18. Lateral Occipital Cortex (B)—Pecuneus (B) 19. Superior Temporal Gyrus (L)—Middle Temporal Gyrus (R) 20. Inferior Temporal Gyrus (B)
 Functional community 1. Hippocampus (B)—Brain-stem (B)—Cerebellum Anterior Loble (B) 2. Amygdala (B)—Parahippocampus (B)—Posterior Temporal Gyrus (R) 3. Orbital Frontal Cortex (B) 4. Striatum (B)—Medial Prefrontal Cortex (B) 5. Anterior Temporal Gyrus (R) 6. Temporal Gyrus (L) 7 Visual Cortex (B) 8. Thalamus (B) 9. Posterior Cingulate Cortex (B)—Precuneus (B)—Angular Cortex (B)—Superior Frontal Gyrus (B)—Parahippocampal Gyrus (B)—Medial Prefrontal Cortex (B)—Middle Temporal Gyrus (B) 10. Dorsal Lateral Prefrontal Cortex (B-L) 11. Middle Frontal Gyrus (R)—IFG (R) 12. Middle Frontal Gyrus (L)—IFG (L) 13. Supplementary Motor Cortex (B)—Middle Cingulate Gyrus (B) 14. Precentral Gyrus (B)—Postcentral Gyrus (B) 15. Anterior Cingulate Cortex (B)—Supramarginal Gyrus (B)—Inferior Parietal Lobule (B)—Superior Temporal Gyrus (B)—Insular Lobe (B) 16. Superior Frontal Gyrus (L)—Middle Frontal Cortex (L) 17. Middle Frontal Gyrus (R)—Inferior Parietal Lobule (B) 18. Lateral Occipital Cortex (B)—Pecuneus (B) 19. Superior Temporal Gyrus (L)—Middle Temporal Gyrus (R) 20. Inferior Temporal Gyrus (B)

Note: B, Bilateral; L, Left; and R, Right.

Figure 1.

The whole-brain functional connectome: communities and layout. The full brain contains about 22 387 4-mm cubic voxels. We detected 20 functional communities (A), which are colored distinctly within multiple axial slices and rendered on a MNI152 standard brain surface (see Table 2 for details on communities and regions). To highlight the overall layout, the voxel-wise connectome was further visualized as a network layout with the same colors (B).

Figure 1.

The whole-brain functional connectome: communities and layout. The full brain contains about 22 387 4-mm cubic voxels. We detected 20 functional communities (A), which are colored distinctly within multiple axial slices and rendered on a MNI152 standard brain surface (see Table 2 for details on communities and regions). To highlight the overall layout, the voxel-wise connectome was further visualized as a network layout with the same colors (B).

### Voxel-wise Centrality Strengths

Voxel-wise centrality analyses allowed us to depict full brain patterns for 4 centrality measures at the 4-mm level. As shown in Figure 2, voxels in the insula exhibited high centrality and high similarity across centrality measures. We plotted voxels against one another with respect to centrality scores (Fig. 3). Except for DC and PC, which again were highly correlated, we observed more complex relationships, with subpopulations of voxels exhibiting notably different relationships. For example, while one population of voxels exhibited relatively high correlation between EC and both DC and PC, another population of voxels exhibited little to no relationship between EC and either of these measures. Such distinctions among voxels likely reflect differences with respect to whether or not a given voxel has a high degree of global connectivity, a feature to which EC is particularly sensitive. The difference map between EC and DC in Supplementary Figure S3 illustrates the distinction between these 2 populations (e.g., medial visual cortex versus posterior cingulate cortex). Finally, the relationship between EC and SC was more parabolic (i.e., “U” shaped), indicating that SC is intermediate between DC and EC.

Figure 2.

VNCM. For each of 4 centralities including DC, SC, EC, and PC, a one-way ANOVA model including center, age, sex, head motion, and registration error as covariates of interest was used to calculate the voxel-wise centrality map across 1003 participants. Multiple comparisons were corrected at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 2.

VNCM. For each of 4 centralities including DC, SC, EC, and PC, a one-way ANOVA model including center, age, sex, head motion, and registration error as covariates of interest was used to calculate the voxel-wise centrality map across 1003 participants. Multiple comparisons were corrected at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 3.

Scatter plots of voxel-wise network centrality measures. For each pair of 4 centrality measures including DC, SC, EC, and PC, a subscatter plot of all voxel centrality scores was used to show the relationship between pairs of centrality measures. The fitted lines were also plotted in red color according to the best least square fit of the centrality scores. The histograms of each centrality score were depicted as in the diagonal scatter plots.

Figure 3.

Scatter plots of voxel-wise network centrality measures. For each pair of 4 centrality measures including DC, SC, EC, and PC, a subscatter plot of all voxel centrality scores was used to show the relationship between pairs of centrality measures. The fitted lines were also plotted in red color according to the best least square fit of the centrality scores. The histograms of each centrality score were depicted as in the diagonal scatter plots.

### Regional Concordance of Centrality

Two analytic approaches were employed to assess concordance of centrality measures assigned to a given region. First, we calculated Kendall’s W and its entropy version by treating the different centrality measures as raters and subjects as samples. Core regions of the default network and dorsal attention network exhibited lower concordance (Fig. 4A) and higher entropy (Fig. 4B) compared with lower order regions. Of note, all centrality measures were highly concordant with respect to differentiating the gray matter/white matter boundary.

Figure 4.

Regional variability of concordance among centrality measures across 1003 participants. Kendall’s W (A) and entropy (B) were calculated to measure concordance among different centrality measures. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 4.

Regional variability of concordance among centrality measures across 1003 participants. Kendall’s W (A) and entropy (B) were calculated to measure concordance among different centrality measures. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

To understand the unique aspects of functional connectivity captured by different centrality measures, we compared the results obtained with each of the centrality measures with those obtained for DC, using voxel-wise (Fig. 5) paired t-tests carried out across participants. We used DC as the reference measure due to its relative simplicity and wide usage. Relative to DC, EC demonstrated significantly higher centrality for subcortical regions, suggesting that while these regions may not be as widely connected throughout the brain, they are connected with key hubs, thereby increasing their own EC. In contrast, DC was greater in cortical regions, suggesting more diffuse patterns of connectivity. It is worth noting that PC, which by definition is intermediate to the more globally focused EC and more locally focused DC measures, indeed yielded results intermediate to EC and DC. Like EC, PC was notably more sensitive to the centrality of subcortical regions. PC also exhibited significantly lower centrality in cortical regions, though not to the extent of EC.

Figure 5.

Distinction between voxel-wise network centrality measures. Using DC as the reference centrality measure, a two-sample paired t-test on centrality scores was performed for each voxel between DC and each of other 3 centrality scores: SC, EC, and PC across 1003 participants. Multiple comparisons were corrected at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 5.

Distinction between voxel-wise network centrality measures. Using DC as the reference centrality measure, a two-sample paired t-test on centrality scores was performed for each voxel between DC and each of other 3 centrality scores: SC, EC, and PC across 1003 participants. Multiple comparisons were corrected at the cluster level using Gaussian random field theory (min Z > 2.3; cluster significance: P < 0.05, corrected). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

When comparing SC with DC, findings were notably distinct from those obtained for both EC and PC. In particular, SC emphasized the centrality of parietal and frontal regions commonly implicated in attentional control (Behrmann et al. 2004), particularly the posterior inferior frontal gyrus (IFG) known to underlie attentional biasing and selection (MacDonald et al. 2000; Miller and Cohen 2001; Milham et al. 2003; Banich 2009).

### TRT Reliability

All voxel-wise centrality measures yielded spatially extended short- (Supplementary Fig. S2) and long-term (Fig. 6) TRT reliability, with short-term reliability outperforming long-term reliability. The 4 measures of centrality showed highly similar pattern of TRT reliability for both binarized (Fig. 6: intersession; Supplementary Fig. S2: intrasession) and weighted networks (Supplementary Fig. S8: intersession; Supplementary Fig. S9: intrasession). These voxel-level findings are consistent with those from a recent study on TRT reliability of large-scale functional brain networks (Wang et al. 2011). Once again, the specific centrality measures were not all equivalent with respect to which regions exhibited higher or lower reliability. DC and PC showed near maximal correlation with each other with respect to short- (r = 0.99) and long-term (r = 0.99) ICC values across voxels (i.e., spatially correlating the ICC maps for DC and PC). This implies that DC and PC have nearly identical reliability. The other correlations with DC were much lower (for EC, short term: r = 0.69; long term: r = 0.70; for SC, short term: r = 0.82; long term: r = 0.84). SC exhibited the lowest degree of correlation of ICC values with the 3 other centrality measures (short term: average r = 0.61; long term: average r = 0.62).

Figure 6.

TRT reliability of voxel-wise network centrality measures. For each of 4 centrality measures including DC, SC, EC, and PC, intersession or long-term ICCs were computed for each voxel. For each centrality measure, all significantly reliable voxels were visualized as multiple axial slices (cluster significance: P < 0.05, corrected for multiple comparisons). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 6.

TRT reliability of voxel-wise network centrality measures. For each of 4 centrality measures including DC, SC, EC, and PC, intersession or long-term ICCs were computed for each voxel. For each centrality measure, all significantly reliable voxels were visualized as multiple axial slices (cluster significance: P < 0.05, corrected for multiple comparisons). The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

### Age and Sex Effects

Before examining specific findings, we first calculated the concordance among the 4 measures with respect to the statistical volumes produced for age and sex effects. For both age and sex, overall statistical maps exhibited a high degree of concordance among DC, EC, and PC measures (age: average Kendall’s W = 0.93; sex: average Kendall’s W = 0.92); more moderate concordance was observed between SC and the other measures, particularly for age effects (DC: Kendall’s W = 0.85; EC: Kendall’s W = 0.72; and PC: Kendall’s W = 0.87).

When examining voxel-level findings related to age (Fig. 7), age-related increases in centrality within medial temporal lobe and caudate regions were highly concordant across all centrality measures. Additionally, age-related decreases in centrality were detected in superior parietal and posterior cingulate regions, though only by DC, SC, and PC. EC was relatively insensitive to age, suggesting preservation of global organization over the age range we examined. Beyond highly similar findings, for both DC and SC, weighted networks (Supplementary Fig. S5) seemed more robust in the detection of age effect than binarized networks (Fig. 7: binarized network).

Figure 7.

Areas exhibiting significant age-related variation in voxel-wise network centrality measures. Group-level maps were derived from one-way ANOVA across 1003 participants from 21 centers (factor: center; covariates: age, sex, head motion, and registration error). All group-level maps depicted were corrected for multiple comparisons at the cluster level using Gaussian random-field theory (min Z > 2.3; cluster significance P < 0.05, corrected). For each of 4 centrality measures including DC, SC, EC, and PC, voxels exhibiting significant effects of age as detected by the one-way ANOVA are depicted. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 7.

Areas exhibiting significant age-related variation in voxel-wise network centrality measures. Group-level maps were derived from one-way ANOVA across 1003 participants from 21 centers (factor: center; covariates: age, sex, head motion, and registration error). All group-level maps depicted were corrected for multiple comparisons at the cluster level using Gaussian random-field theory (min Z > 2.3; cluster significance P < 0.05, corrected). For each of 4 centrality measures including DC, SC, EC, and PC, voxels exhibiting significant effects of age as detected by the one-way ANOVA are depicted. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Only the hippocampus showed consistent sex differences across all 4 measures with greater centrality being observed in females than males (Fig. 8). Consistent effects were detected in multiple regions for 3 of 4 centrality measures. The specific combination of measures varied somewhat. For example, lateral precentral and postcentral lobule regions exhibited greater DC, EC, and PC, but not SC centrality in males. In contrast, medial occipital areas exhibited greater DC, SC, and PC, but not EC centrality in females. Intriguingly, the posterior cingulate cortex showed detectable but inconsistent sex findings across 3 of the 4 measures—that is, females exhibited greater DC and EC, but lower SC. In contrast, weighted SC did not detect the significantly higher centrality of PCC in females than males (Supplementary Fig. S7).

Figure 8.

Areas exhibiting significant sex-related variation in voxel-wise network centrality measures. Group-level maps were derived from one-way ANOVA across 1003 participants from 21 centers (factor: center; covariates: age, sex, head motion, and registration error). All group-level maps depicted were corrected for multiple comparisons at the cluster level using Gaussian random-field theory (min Z > 2.3; cluster significance P < 0.05, corrected). For each of 4 centrality measures including DC, SC, EC, and PC, voxels exhibiting significant effects of sex as detected by one-way ANOVA are depicted. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Figure 8.

Areas exhibiting significant sex-related variation in voxel-wise network centrality measures. Group-level maps were derived from one-way ANOVA across 1003 participants from 21 centers (factor: center; covariates: age, sex, head motion, and registration error). All group-level maps depicted were corrected for multiple comparisons at the cluster level using Gaussian random-field theory (min Z > 2.3; cluster significance P < 0.05, corrected). For each of 4 centrality measures including DC, SC, EC, and PC, voxels exhibiting significant effects of sex as detected by one-way ANOVA are depicted. The final statistical maps are visualized as 4 hemispheric surfaces (cortical regions: 2 for left hemisphere [L] and 2 for right hemisphere [R]) and 8 coronal slices (subcortical regions: y coordinates from −23.5 to 11.5 mm with 5 mm space).

Our ANOVA model includes head motion and registration error as covariates. As expected, the spatial extents of age effects were markedly reduced, implying the potential impacts of the 2 factors on centrality changes with age. However, significant age-centrality relationships were still observed (Fig. 7: head motion and registration error as covariates; Supplementary Fig. S4: no head motion and registration error as covariates). Of note, weighted versions of centrality seem more robust against these 2 confounding factors (Fig. 7: binary network; Supplementary Fig. S5: weighted network). In contrast, the 2 factors had no obvious impact on sex-centrality relationship (Fig. 8: head motion and registration error as covariates; Supplementary Fig. S6: no head motion and registration error as covariates).

## Discussion

Investigating voxel-wise centrality maps provided novel insights into the patterns and complexity of functional connectivity throughout the huge human functional connectome. Despite marked differences in definitions among centrality measures, prior studies of centrality have treated them as if they were interchangeable, and generally, similar patterns of results have been reported (Fransson and Marrelec 2008; Hagmann et al. 2008; Buckner et al. 2009; He et al. 2009; Fransson et al. 2011). In the present work, we first assembled and visualized the 4-mm resolution functional connectome as a functional network with 22 387 nodes and then analyzed the contribution of each node to connectome-wide functional connectivity. Our findings highlight the importance of considering multiple (i.e., local, meso, and global) scales of connectivity properties in the functional connectome and distinguishing them in terms of their interindividual variability, TRT reliability, and phenotypic sensitivity. We argue that using voxel-wise centrality as measures of whole-brain functional connectivity will facilitate discovery of physiological mechanisms underlying the brain’s intrinsic functional organization.

### Functional Communities within the Functional Connectome

The construction of whole-brain functional network is the first step of computing centrality measures. Particularly, at the voxel level, the network and its organization are not commonly visualized due to its size. It is important to illustrate the complex network and its organization for readers to get an appreciation of the graph on which the centrality and/or functional connectivity measures will be calculated and to guide its further explorations. We first built and visualized the network layout of the functional connectome at 4 mm spatial resolution, displaying the high complexity of the network. Community detection revealed 20 functional modules, parcellating the brain in terms of functional interactions among voxels on the basis of spontaneous low-frequency fluctuations.

Beyond the networks consistent with those found in previous studies on functional brain parcellations (Peltier et al. 2003; Bellec et al. 2006; van den Heuvel, Mandl, et al. 2008; He et al. 2009; Mezer et al. 2009) (e.g., default network and control network), the communities detected appear to reflect the fundamental organization of the functional connectome. The algorithm produced 20 communities or modules spanning across 40 separate brain regions. Clearly, such a functional parcellation is different from a structural parcellation (Tzourio-Mazoyer et al. 2002), reflecting both topological and hierarchical distinctions between structure and function of the brain. Compared with 3 recent studies on functional parcellation of the human brain (Craddock et al. 2011; Doucet et al. 2011; Yeo et al. 2011), there is notable convergence regarding the basic organization of the brain. For example, all these studies revealed several common functional modules of sensorimotor, auditory, and visual cortex as well as default network and attention control network, which are bilaterally distributed (i.e., homotopic communities). Note that more accurate comparison with these previous studies is difficult, as they employed various clustering methods to build modular organization while we considered a community detection algorithm from network science. In the sense of hierarchy in algorithm, the work of Doucet et al. (2011) was more comparable with ours. We found a close match for the number of modules (20) to theirs (23). However, our method revealed some unique lateralized communities (i.e., nonhomotopic communities) that are commonly associated with functions, such as language, working memory, and biasing attention (Toga and Thompson 2003; Yan et al. 2009; Zuo, Kelly, Di Martino, et al. 2010; Gee et al. 2011). We believe that this difference is likely related to a methodological choice. In sum, the current functional parcellation based on a community detection algorithm provides a unique way of mapping and deepening understanding of functional hierarchy within the human brain (Smith et al. 2009).

### The Importance of Multiple-scale Connectivity Properties: Network Centrality Mapping

While DC, a measure of local (direct) network connectivity, highlights the importance of higher order cortical association regions, it is less sensitive to paralimbic and subcortical regions. However, when global properties are taken into account by centrality measures, such as EC and PC, hippocampal, striatal, and thalamic regions become prominent, likely due to their extensive global associations (Achard et al. 2006; Sporns et al. 2007).

At mesoscale connectivity, SC particularly emphasizes the centrality of parietal and frontal regions commonly implicated in attentional control (Behrmann et al. 2004). Interestingly, the centrality of the posterior IFG, a key region implicated in top-down control of attention (Milham et al. 2003; Behrmann et al. 2004; Banich 2009), is only appreciated by SC. Further research will have to indicate whether such divergent findings reflect diminished informational importance of this region or rather something unique about the nature of its information flow.

### The Importance of Multiple-scale Connectivity Properties: Relating Phenotypic Information to Network Centrality

A major feature of centrality measures is their inherently exploratory nature—no a priori predictions are required. This makes them ideally suited for exploration of the neural correlates of both dimensional and categorical phenotypic data (Di Martino et al. 2009; Cox et al. 2010). Here, we were able to detect robustly significant if complex associations for both age and sex.

Even when the focus is limited to the cerebral cortex, the information emerging from global and local measures can differ substantially. For instance, when examining age effects, DC revealed age-related decreases in centrality within precuneus and posterior cingulate regions, but EC did not. Although such disparities seem problematic at first glance, they are informative. In this case, we interpret these findings as supporting the hypothesis that while the local or (direct) connectivity of these regions decreases with age, they maintain their connections with hub-like regions within the brain and thus their centrality at a global level (Park and Reuter-Lorenz 2009). Interestingly, the PC measure was potentially most effective for capturing centrality at both the global and local level, likely due to the inclusion of “random jumps” that allow the otherwise global algorithm to better appreciate local distinctions.

Assessing sex-related differences in centrality revealed a complex range of results at the mesoscale. Females exhibited greater centrality than males for all 4 measures in the hippocampus. However, for other cortical regions, the direction of sex-related differences depended on the specific centrality measure. If our results are replicated independently, they would motivate an improved understanding of how to integrate higher order measures and models across multiple scales.

### Looking Beyond Functional Hubs

Following the structural literature, the functional connectivity community has begun its own hunt for hubs (reviewed in Wang et al. 2010). Yet, the definition of a functional hub is poorly operationalized, at least when functional connectivity is based upon simple cross-correlation. Our results highlight the potential utility of using centrality as a way of characterizing information flow through the full brain functional connectome rather than attempting to identify functional hubs. We did not find functional hubs emerging consistently across centrality measures—at least not at the individual subject level. Heteromodal association areas in the cortex showed the lowest concordances (i.e., highest entropy) across centrality measures. This attests to the complexity and richness of these regions—presumably not to their lack of importance in information flow.

### Limitations and Future Directions

We note several limitations to be considered when interpreting our findings. A potential criticism concerns discrepancies in analytical strategies and graph theory approaches. For example, various strategies can be used for thresholding a functional connectivity matrix to format a graph adjacency matrix (Wang et al. 2010): 1) correlation criterion using a fixed correlation value for all subjects, 2) sparsity criterion using a fixed sparsity value for all subjects, or 3) area under the curve criterion to combine results based on multiple thresholds. In the current study, we use a fourth alternative that we consider a “significance criterion”—that is, maintaining a constant probability of false positives across samples. This criterion is motivated by our use of 1003 R-fMRI data sets from 21 centers, each comprising a different number of scan volumes. Although we conduct analyses at 2 different thresholds of significance, the degree to which our results will be robust to a broader range of significance criteria requires additional investigation. Of note, more consistent results might be obtained if graph analyses could be carried out on fully weighted (i.e., nonthresholded) functional networks. Initial efforts in this regard have been directed at characterizing functional modules (Rubinov and Sporns 2011) and may be extendable to assessing nodal centrality and nodal influence given a module partition.

Centrality analyses also raise challenging issues regarding structure–function relationships. Although high centrality regions (i.e., network hubs) have been consistently identified by both structural and functional data (Fransson and Marrelec 2008; Hagmann et al. 2008; Buckner et al. 2009; Gong et al. 2009; Fransson et al. 2011), differences have also been noted (He et al. 2009). Such differences are to be expected since structural and functional connectivity do not show simple 1:1 correspondence (Vincent et al. 2007; Margulies et al. 2009; Roy et al. 2009). In addition, graph metrics applied to structural and functional networks are amenable to different neurobiological interpretations. Future studies should clarify centrality relationships by combining structural and functional connectomes based on high-resolution neuroimaging data.

There are 2 potential confounds in analyzing multisite neuroimaging data and examining age effects. First, artificial effects could be introduced if the age range is not equally distributed across different sites, which is particularly relevant in the current situation as most sites had young participants while only a few sites had a small number of old participants. Such a site effect cannot be totally ruled out by including center as a factor in the model. We thus performed the same analysis based upon the data from a single site (Montreal), which has the largest age span (19–85 years). Limited by the number of samples (51 subjects), a liberal uncorrected Z-statistic threshold (Z > 1.0) was used to allow visualization of subthreshold trends present for the age effect on centrality. As indicated in Supplementary Figs S10 and S11, an overall trend of age-centrality relationship was observed, mitigating this concern. Second, there is growing realization that motion has an impact on ICNs (Van Dijk et al. 2011). We had demonstrated similar effects of head motion and registration quality on a life span developmental trajectory study (Zuo, Kelly, Di Martino, et al. 2010). Similarly, we updated our centrality analyses by including head motion and registration error as covariates. As expected, the spatial extents of age effects were notably reduced, implying a potential impact of these 2 factors on centrality changes with age. However, the significant age-centrality relationship was kept, implying such an effect is unlikely to be fully explained by the 2 confounding factors.

Future work on centrality analysis with directed graphs (i.e., graphs in which temporal order or causal relationships can be inferred or demonstrated) derived from R-fMRI data would yield crucial insights into the information flow through the functional connectome. We performed our analyses using undirected graphs. Although establishing a “true” network model for the brain connectome is difficult (Smith et al. 2011), recent work suggests that it may be possible to map out a directed graphical representation of the causal relationships between brain regions using resting fMRI data (Friston et al. 2011). Our approach of detecting a limited number of centrality communities may provide candidates for directed graphical analyses using data-driven dynamical causal models (Friston et al. 2011). Such analyses can be first carried out for region of interest–derived directed networks. Of note, voxel-level centrality analyses on directed graphs are particularly challenging. Region of interest–based analyses, on the other hand, can open the door to graph theoretic analyses that leverage advances in centrality-based characteristics of information flow in directed networks (Shepelyansky and Zhirov 2010).

## Conclusion

We demonstrate that each of 4 centrality measures captures a different aspect of functional connectivity, yielding significant differences in both their relative sensitivities (e.g., global vs. local properties) and utility (e.g., TRT reliability). Our findings highlight the importance of considering multiple-scale connectivity properties when examining the functional connectome. In going beyond detecting functional hubs, we treat centrality as a way to chart information flow in the functional connectome and demonstrate its power to capture phenotypic correlates (i.e., age and sex). In sum, our work shows that combining multiple centrality measures allows a more comprehensive characterization of the complexities of the functional connectome and could serve to motivate studies on the physiological mechanisms underlying functional organization in neurodegenerative and psychiatric disorders (Greicius 2008; Broyd et al. 2009; Seeley et al. 2009; Seeley 2010; Cole et al. 2011).

## Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/

## Funding

Dr X.-N.Z. acknowledges support from Startup Foundation for Distinguished Research Professor of Institute for Psychology (Y0CX492S03), the Natural Science Foundation of China (81171409), and Dr O.S. acknowledges support from the J.S. McDonnell Foundation.

The authors would like to thank Dr Ernesto Estrada for valuable discussion on subgraph centrality of very large graphs, Drs Adriana Di Martino and Stan Colcombe for editorial comments, and Drs Clare Kelly and Christine Cox for feedback on using centrality analyses scripts. Conflict of Interest : None declared.

## References

Achard
S
R
Whitcher
B
Suckling
J
Bullmore
E
A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs
J Neurosci
,
2006
, vol.
26
(pg.
63
-
72
)
Banich
MT
Executive function: the search for an integrated account
Curr Dir Psychol Sci
,
2009
, vol.
18
(pg.
89
-
94
)
Beckmann
CF
DeLuca
M
Devlin
JT
Smith
SM
Investigations into resting-state connectivity using independent component analysis
Philos Trans R Soc Lond B Biol Sci
,
2005
, vol.
360
(pg.
1001
-
1013
)
Behrmann
M
Geng
JJ
Shomstein
S
Parietal cortex and attention
Curr Opin Neurobiol
,
2004
, vol.
14
(pg.
212
-
217
)
Bellec
P
Perlbarg
V
Jbabdi
S
Pelegrini-Issac
M
Anton
JL
Doyon
J
Benali
H
Identification of large-scale networks in the brain using fMRI
Neuroimage
,
2006
, vol.
29
(pg.
1231
-
1243
)
Biswal
B
Yetkin
FZ
Haughton
VM
Hyde
JS
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
Magn Reson Med
,
1995
, vol.
34
(pg.
537
-
541
)
Biswal
BB
Mennes
M
Zuo
XN
Gohel
S
Kelly
C
Smith
SM
Beckmann
CF
JS
Buckner
RL
Colcombe
S
, et al.  .
Toward discovery science of human brain function
Proc Natl Acad Sci U S A
,
2010
, vol.
107
(pg.
4734
-
4739
)
Blondel
VD
Guillaume
JL
Lambiotte
R
Lefebvre
E
Fast unfolding of communities in large networks
J Stat Mech
,
2008
, vol.
2008
pg.
P10008

Boldi
P
Santini
M
Vigna
S
PageRank: functional dependencies
ACM Trans Inf Syst
,
2009
, vol.
27
(pg.
1
-
23
)
Bonacich
P
Factoring and weighting approaches to clique identification
J Math Sociol
,
1972
, vol.
2
(pg.
113
-
120
)
Broyd
SJ
Demanuele
C
Debener
S
Helps
SK
James
CJ
Sonuga-Barke
EJ
Default-mode brain dysfunction in mental disorders: a systematic review
Neurosci Biobehav Rev
,
2009
, vol.
33
(pg.
279
-
296
)
Buckner
RL
Sepulcre
J
Talukdar
T
Krienen
FM
Liu
H
Hedden
T
Andrews-Hanna
JR
Sperling
RA
Johnson
KA
Cortical hubs revealed by intrinsic functional connectivity: mapping, assessment of stability, and relation to Alzheimer’s disease
J Neurosci
,
2009
, vol.
29
(pg.
1860
-
1873
)
Bullmore
E
Sporns
O
Complex brain networks: graph theoretical analysis of structural and functional systems
Nat Rev Neurosci
,
2009
, vol.
10
(pg.
186
-
198
)
Carhart-Harris
RL
Friston
KJ
The default-mode, ego-functions and free-energy: a neurobiological account of Freudian ideas
Brain
,
2010
, vol.
133
(pg.
1265
-
1283
)
Cole
DM
Smith
SM
Beckmann
CF
Advances and pitfalls in the analysis and interpretation of resting-state FMRI data
Front Syst Neurosci
,
2010
, vol.
4
pg.
8

Cole
MW
Anticevic
A
Repovs
G
Barch
D
Variable global dysconnectivity and individual differences in schizophrenia
Biol Psychiatry.
,
2011
, vol.
70
(pg.
43
-
50
)
Cole
MW
Pathak
S
Schneider
W
Identifying the brain’s most globally connected regions
Neuroimage
,
2010
, vol.
49
(pg.
3132
-
3148
)
Cox
CL
Gotimer
K
Roy
AK
Castellanos
FX
Milham
MP
Kelly
C
PLoS One
,
2010
, vol.
5
pg.
e12296

RC
James
GA
Holtzheimer
PE
Hu
XP
Mayberg
HS
A whole brain fMRI atlas generated via spatially constrained spectral clustering
Hum Brain Mapp
,
2011

doi:10.1002/hbm.21333
Crofts
JJ
Higham
DJ
A weighted communicability measure applied to complex brain networks
J R Soc Interface
,
2009
, vol.
6
(pg.
411
-
414
)
Damoiseaux
JS
Rombouts
SA
Barkhof
F
Scheltens
P
Stam
CJ
Smith
SM
Beckmann
CF
Consistent resting-state networks across healthy subjects
Proc Natl Acad Sci U S A
,
2006
, vol.
103
(pg.
13848
-
13853
)
Deco
G
Jirsa
VK
McIntosh
AR
Emerging concepts for the dynamical organization of resting-state activity in the brain
Nat Rev Neurosci
,
2011
, vol.
12
(pg.
43
-
56
)
DeFelipe
J
From the connectome to the synaptome: an epic love story
Science
,
2010
, vol.
330
(pg.
1198
-
1201
)
Di Martino
A
Z
Kelly
C
Roy
AK
Gee
DG
Uddin
LQ
Gotimer
K
Klein
DF
Castellanos
FX
Milham
MP
Relationship between cingulo-insular functional connectivity and autistic traits in neurotypical adults
Am J Psychiatry
,
2009
, vol.
166
(pg.
891
-
899
)
Doucet
G
Naveau
M
Petit
L
Delcroix
N
Zago
L
Crivello
F
Jobard
G
Tzourio-Mazoyer
N
Mazoyer
B
Mellet
E
, et al.  .
Brain activity at rest: a multiscale hierarchical functional organization
J Neurophysiol
,
2011
, vol.
105
(pg.
2753
-
2763
)
E
Higham
DJ
Network properties revealed through matrix functions
SIAM Rev
,
2010
, vol.
52
(pg.
696
-
714
)
E
Rodriguez-Velazquez
JA
Subgraph centrality in complex networks
Phys Rev E Stat Nonlin Soft Matter Phys
,
2005
, vol.
71
pg.
056103

Fair
DA
Cohen
AL
Dosenbach
NU
Church
JA
Miezin
FM
Barch
DM
Raichle
ME
Petersen
SE
Schlaggar
BL
The maturing architecture of the brain’s default network
Proc Natl Acad Sci U S A
,
2008
, vol.
105
(pg.
4028
-
4032
)
Fortunato
S
Boguñá
M
Flammini
A
Menczer
F
Aiello
W
Broder
A
Janssen
J
Milios
E
Approximating pagerank from in-degree
Fourth international workshop on algorithms and models for the Web-Graph
,
2008
Springer-Verlag
(pg.
59
-
71
)
Fox
MD
Raichle
ME
Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging
Nat Rev Neurosci
,
2007
, vol.
8
(pg.
700
-
711
)
Fransson
P
U
Blennow
M
Lagercrantz
H
The functional architecture of the infant brain as revealed by resting-state fMRI
Cereb Cortex
,
2011
, vol.
21
(pg.
145
-
154
)
Fransson
P
Marrelec
G
The precuneus/posterior cingulate cortex plays a pivotal role in the default mode network: evidence from a partial correlation network analysis
Neuroimage
,
2008
, vol.
42
(pg.
1178
-
1184
)
Friston
K
The free-energy principle: a unified brain theory?
Nat Rev Neurosci
,
2010
, vol.
11
(pg.
127
-
138
)
Friston
KJ
Dolan
RJ
Computational and dynamic models in neuroimaging
Neuroimage
,
2010
, vol.
52
(pg.
752
-
765
)
Friston
KJ
Li
B
Daunizeau
J
Stephan
KE
Network discovery with DCM
Neuroimage
,
2011
, vol.
56
(pg.
1202
-
1221
)
Gee
DG
Biswal
BB
Kelly
C
Stark
DE
Margulies
DS
Z
Uddin
LQ
Klein
DF
Banich
MT
Castellanos
FX
, et al.  .
Low frequency fluctuations reveal integrated and segregated processing among the cerebral hemispheres
Neuroimage
,
2011
, vol.
54
(pg.
517
-
527
)
Gleich
DF
Gray
AP
Greif
C
Lau
T
An inner-outer iteration for computing pagerank
SIAM J Sci Comput
,
2010
, vol.
32
(pg.
349
-
371
)
Gong
G
Rosa-Neto
P
Carbonell
F
Chen
ZJ
He
Y
Evans
AC
Age- and gender-related differences in the cortical anatomical network
J Neurosci
,
2009
, vol.
29
(pg.
15684
-
15693
)
Greicius
M
Resting-state functional connectivity in neuropsychiatric disorders
Curr Opin Neurol
,
2008
, vol.
21
(pg.
424
-
430
)
Hagmann
P
Cammoun
L
Gigandet
X
Meuli
R
Honey
CJ
Wedeen
VJ
Sporns
O
Mapping the structural core of human cerebral cortex
PLoS Biol
,
2008
, vol.
6
pg.
e159

He
Y
Wang
J
Wang
L
Chen
ZJ
Yan
C
Yang
H
Tang
H
Zhu
C
Gong
Q
Zang
Y
, et al.  .
Uncovering intrinsic modular organization of spontaneous brain activity in humans
PLoS One
,
2009
, vol.
4
pg.
e5226

Joyce
KE
Laurienti
PJ
Burdette
JH
Hayasaka
S
A new measure of centrality for brain networks
PLoS One
,
2010
, vol.
5
pg.
e12200

Koschützki
D
Lehmann
KA
Peeters
L
Richter
S
Tenfelde-Podehl
D
Zlotowski
O
Brandes
U
Erlebach
T
Centrality indices
Network analysis: methodological foundations
,
2005
New York
Springer-Verlag
(pg.
16
-
61
)
Legendre
P
Species associations: Kendall coefficient of concordance revisited
J Agric Biol Environ Stat
,
2005
, vol.
10
(pg.
226
-
245
)
Lehoucq
RB
Sorensen
DC
Deflation techniques for an implicitly restarted Arnoldi iteration
SIAM J Matrix Anal Appl
,
1996
, vol.
17
(pg.
789
-
821
)
Lohmann
G
Margulies
DS
Horstmann
A
Pleger
B
Lepsien
J
Goldhahn
D
Schloegl
H
Stumvoll
M
Villringer
A
Turner
R
Eigenvector centrality mapping for analyzing connectivity patterns in fMRI data of the human brain
PLoS One
,
2010
, vol.
5
pg.
e10232

MacDonald
AW
Cohen
JD
Stenger
VA
Carter
CS
Dissociating the role of the dorsolateral prefrontal and anterior cingulate cortex in cognitive control
Science
,
2000
, vol.
288
(pg.
1835
-
1838
)
Margulies
DS
Kelly
AM
Uddin
LQ
Biswal
BB
Castellanos
FX
Milham
MP
Mapping the functional connectivity of anterior cingulate cortex
Neuroimage
,
2007
, vol.
37
(pg.
579
-
588
)
Margulies
DS
Vincent
JL
Kelly
C
Lohmann
G
Uddin
LQ
Biswal
BB
Villringer
A
Castellanos
FX
Milham
MP
Petrides
M
Precuneus shares intrinsic functional architecture in humans and monkeys
Proc Natl Acad Sci U S A
,
2009
, vol.
106
(pg.
20069
-
20074
)
Martin
S
Brown
WM
Klavans
R
Boyack
K
Wong
PC
Park
J
Hao
MC
Chen
C
Börner
K
Kao
DL
Roberts
JC
OpenOrd: an open-source toolbox for large graph layout
SPIE conference on visualization and data analysis
,
2011
San Francisco (CA)
SPIE
Mezer
A
Yovel
Y
Pasternak
O
Gorfine
T
Assaf
Y
Cluster analysis of resting-state fMRI time series
Neuroimage
,
2009
, vol.
45
(pg.
1117
-
1125
)
Milham
MP
Banich
MT
V
Competition for priority in processing increases prefrontal cortex’s involvement in top-down control: an event-related fMRI study of the stroop task
Brain Res Cogn Brain Res
,
2003
, vol.
17
(pg.
212
-
222
)
Miller
EK
Cohen
JD
An integrative theory of prefrontal cortex function
Annu Rev Neurosci
,
2001
, vol.
24
(pg.
167
-
202
)
Park
DC
Reuter-Lorenz
P
The adaptive brain: aging and neurocognitive scaffolding
Annu Rev Psychol
,
2009
, vol.
60
(pg.
173
-
196
)
Peltier
SJ
Polk
TA
Noll
DC
Detecting low-frequency functional connectivity in fMRI using a self-organizing map (SOM) algorithm
Hum Brain Mapp
,
2003
, vol.
20
(pg.
220
-
226
)
Power
JD
Fair
DA
Schlaggar
BL
Petersen
SE
The development of human functional brain networks
Neuron
,
2010
, vol.
67
(pg.
735
-
748
)
Roy
AK
Z
Margulies
DS
Kelly
AM
Uddin
LQ
Gotimer
K
Biswal
BB
Castellanos
FX
Milham
MP
Functional connectivity of the human amygdala using resting state fMRI
Neuroimage
,
2009
, vol.
45
(pg.
614
-
626
)
Rubinov
M
Sporns
O
Complex network measures of brain connectivity: uses and interpretations
Neuroimage
,
2010
, vol.
52
(pg.
1059
-
1069
)
Rubinov
M
Sporns
O
Weight-conserving characterization of complex functional brain networks
Neuroimage
,
2011
, vol.
56
(pg.
2068
-
2079
)
Seeley
WW
Anterior insula degeneration in frontotemporal dementia
Brain Struct Funct
,
2010
, vol.
214
(pg.
465
-
475
)
Seeley
WW
Crawford
RK
Zhou
J
Miller
BL
Greicius
MD
Neurodegenerative diseases target large-scale human brain networks
Neuron
,
2009
, vol.
62
(pg.
42
-
52
)
Z
Kelly
AM
Reiss
PT
Gee
DG
Gotimer
K
Uddin
LQ
Lee
SH
Margulies
DS
Roy
AK
Biswal
BB
, et al.  .
The resting brain: unconstrained yet reliable
Cereb Cortex
,
2009
, vol.
19
(pg.
2209
-
2229
)
Shepelyansky
DL
Zhirov
OV
Phys Lett A
,
2010
, vol.
374
(pg.
3206
-
3209
)
Smith
SM
Fox
PT
Miller
KL
Glahn
DC
Fox
PM
Mackay
CE
Filippini
N
Watkins
KE
Toro
R
Laird
AR
, et al.  .
Correspondence of the brain’s functional architecture during activation and rest
Proc Natl Acad Sci U S A
,
2009
, vol.
106
(pg.
13040
-
13045
)
Smith
SM
Miller
KL
Salimi-Khorshidi
G
Webster
M
Beckmann
CF
Nichols
TE
Ramsey
JD
Woolrich
MW
Network modelling methods for FMRI
Neuroimage
,
2011
, vol.
54
(pg.
875
-
891
)
Sporns
O
Small-world connectivity, motif composition, and complexity of fractal neuronal connections
Biosystems
,
2006
, vol.
85
(pg.
55
-
64
)
Sporns
O
Honey
CJ
Kotter
R
Identification and classification of hubs in brain networks
PLoS One
,
2007
, vol.
2
pg.
e1049

Sporns
O
Tononi
G
Kotter
R
The human connectome: a structural description of the human brain
PLoS Comput Biol
,
2005
, vol.
1
pg.
e42

Toga
AW
Thompson
PM
Mapping brain asymmetry
Nat Rev Neurosci
,
2003
, vol.
4
(pg.
37
-
48
)
Tomasi
D
Volkow
ND
Functional connectivity density mapping
Proc Natl Acad Sci U S A
,
2010
, vol.
107
(pg.
9885
-
9890
)
Tomasi
D
Volkow
ND
Association between functional connectivity hubs and brain networks
Cereb Cortex
,
2011
, vol.
21
(pg.
2003
-
2013
)
Tomasi
D
Volkow
ND
Gender differences in brain functional connectivity density
Hum Brain Mapp
,
2011

doi:10.1002/hbm.21252
Tzourio-Mazoyer
N
Landeau
B
Papathanassiou
D
Crivello
F
Etard
O
Delcroix
N
Mazoyer
B
Joliot
M
Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain
Neuroimage
,
2002
, vol.
15
(pg.
273
-
289
)
van den Heuvel
M
Mandl
R
Hulshoff Pol
H
Normalized cut group clustering of resting-state fMRI data
PLoS One
,
2008
, vol.
3
pg.
e2001

van den Heuvel
MP
Stam
CJ
Boersma
M
Hulshoff Pol
HE
Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain
Neuroimage
,
2008
, vol.
43
(pg.
528
-
539
)
Van Dijk
KR
Sabuncu
MR
Buckner
RL
The influence of head motion on intrinsic functional connectivity MRI
Neuroimage
,
2011

doi:10.1016/j.neuroimage.2011.07.044
Vincent
JL
Patel
GH
Fox
MD
Snyder
AZ
Baker
JT
Van Essen
DC
Zempel
JM
Snyder
LH
Corbetta
M
Raichle
ME
Intrinsic functional architecture in the anaesthetized monkey brain
Nature
,
2007
, vol.
447
(pg.
83
-
86
)
Wang
J
Wang
L
Zang
Y
Yang
H
Tang
H
Gong
Q
Chen
Z
Zhu
C
He
Y
Parcellation-dependent small-world brain functional networks: a resting-state fMRI study
Hum Brain Mapp
,
2009
, vol.
30
(pg.
1511
-
1523
)
Wang
J
Zuo
X
He
Y
Graph-based network analysis of resting-state functional MRI
Front Syst Neurosci
,
2010
, vol.
4
pg.
16

Wang
J
Zuo
XN
Gohel
S
Milham
MP
Biswal
BB
He
Y
Graph theoretical analysis of functional brain networks: test-retest evaluation on short- and long-term resting-state functional MRI data
PLoS One
,
2011
, vol.
6
pg.
e21976

Yan
H
Zuo
XN
Wang
D
Wang
J
Zhu
C
Milham
MP
Zhang
D
Zang
Y
Hemispheric asymmetry in cognitive division of anterior cingulate cortex: a resting-state functional connectivity study
Neuroimage
,
2009
, vol.
47
(pg.
1579
-
1589
)
Yeo
BT
Krienen
FM
Sepulcre
J
Sabuncu
MR
Lashkari
D
M
Roffman
JL
Smoller
JW
Zollei
L
Polimeni
JR
, et al.  .
The organization of the human cerebral cortex estimated by functional connectivity
J Neurophysiol
,
2011
, vol.
106
(pg.
1125
-
1165
)
Zuo
XN
Di Martino
A
Kelly
C
ZE
Gee
DG
Klein
DF
Castellanos
FX
Biswal
BB
Milham
MP
The oscillating brain: complex and reliable
Neuroimage
,
2010
, vol.
49
(pg.
1432
-
1445
)
Zuo
XN
Kelly
C
JS
Klein
DF
Castellanos
FX
Milham
MP
Reliable intrinsic connectivity networks: test-retest evaluation using ICA and dual regression approach
Neuroimage
,
2010
, vol.
49
(pg.
2163
-
2177
)
Zuo
XN
Kelly
C
Di Martino
A
Mennes
M
Margulies
DS
Bangaru
S
R
Evans
AC
Zang
YF
Castellanos
FX
, et al.  .
Growing together and growing apart: regional and sex differences in the lifespan developmental trajectories of functional homotopy
J Neurosci
,
2010
, vol.
30
(pg.
15034
-
15043
)