## 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.

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\u2264i,j\u2264N$ (*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\u2264i\u2264N,\u20021\u2264j\u2264N$ of a graph $G$:

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),\u2002V=(vi)1\u2264i\u2264N$ defines all nodes and $A=(aij)1\u2264i\u2264N,\u20021\u2264j\u2264N$ 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),\u20021\u2264i\u2264N$.

#### 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).

#### 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,

In equation (4), $\mu j(i)$ is the *i*-th component of the *j*-th eigenvector of the adjacency matrix *A*, and $\lambda j$ is the corresponding *j*-th eigenvalue, that is, $A\mu j=\lambda j\mu 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\u22121/2AD\u22121/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):

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 $\lambda 1$ (called the principal eigenvalue):

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):

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:

*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):

In this equation, MS_{p} is the between-participants mean square and MS_{e} 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:

*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):

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).

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.

### 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.

### 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.

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.

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).

### 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).

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).

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.