Abstract

Functional hubs are brain regions that play a crucial role in facilitating communication among parallel, distributed brain networks. The developmental emergence and stability of hubs, however, is not well understood. The current study used measures of network topology drawn from graph theory to investigate the development of functional hubs in 99 participants, 10–20 years of age. We found that hub architecture was evident in late childhood and was stable from adolescence to early adulthood. Connectivity between hub and non-hub (“spoke”) regions, however, changed with development. From childhood to adolescence, the strength of connections between frontal hubs and cortical and subcortical spoke regions increased. From adolescence to adulthood, hub–spoke connections with frontal hubs were stable, whereas connectivity between cerebellar hubs and cortical spoke regions increased. Our findings suggest that a developmentally stable functional hub architecture provides the foundation of information flow in the brain, whereas connections between hubs and spokes continue to develop, possibly supporting mature cognitive function.

Introduction

The brain is a network formed by hierarchical interconnections among neurons, circuits, columns, functional regions, and neural systems. Connections at multiple levels of organization give rise to complex neural dynamics that support cognition and adaptive behaviors (Bassett and Gazzaniga 2011). Important refinements in brain connectivity continue to mature through adolescence, possibly supporting cognitive development (Fair et al. 2009; Kelly et al. 2009; Stevens et al. 2009; Asato et al. 2010; Dosenbach et al. 2010), and abnormal brain connectivity has been implicated in several psychiatric and developmental disorders (Geschwind and Levitt 2007; He et al. 2008; Seeley et al. 2009). Therefore, improving our understanding of how brain networks develop may yield important insights into the function and structure of the human brain, which could be of considerable theoretical and translational significance.

Functional brain networks characterize the pattern of functional interactions among brain regions using a single connectivity matrix, where each cell in the matrix describes the functional connectivity strength between one brain region and another (i.e. the statistical association between brain activity time-series). The network topology (i.e. patterns of interconnection) of the brain can be effectively studied using graph theory, which models brain regions as network nodes and connections among regions as edges (Bullmore and Sporns 2009) with varying connectivity strength. Graph-theory analyses of functional brain networks can therefore reveal how brain regions communicate and characterize global and local network properties that support inter-regional interactions.

Prior studies using graph theory have established that functional brain networks exhibit scale-free dynamics (Eguiluz et al. 2005; van den Heuvel et al. 2008), where connectivity distributions follow a power-law scaling (Barabasi and Albert 1999). Functional brain networks also exhibit small-world architecture (Achard et al. 2006; Bassett et al. 2006), in which networks are composed of highly clustered sub-networks interconnected by a few shortcuts (Watts and Strogatz 1998). These network characteristics have important functional consequences: Scale-free dynamics give rise to “hubs” that are critical in mediating information communication across distributed cortical regions (Sporns et al. 2007), whereas small-world organization supports efficient information processing both globally and locally (Sporns et al. 2004; Fair et al. 2009). Although there is evidence of small-world brain networks in infants, children, and adolescents (Fair et al. 2009; Supekar et al. 2009; Fan et al. 2011), how hub organization emerges developmentally is not yet known.

Similar to airline hubs, functional hubs are brain regions that are highly connected with a large number of regions, and they serve as way stations that direct high volumes of information traffic. As such, the hub architecture of the brain may serve as a foundational backbone supporting communication among functionally specialized networks. Although hub properties have been quantified and identified in adults using graph-theoretic centrality measures (Hagmann et al. 2008; Buckner et al. 2009; Cole et al. 2010), the development of functional hubs is not well understood. One study demonstrated that the spatial distribution of functional hubs in infants is significantly different from that of adults (Fransson et al. 2011). Yet, how functional hubs develop from childhood into adulthood, a period of time when functional brain networks undergo major development (Fair et al. 2009), has not been systematically investigated. The goal of the current study was to investigate developmental changes in the architecture of functional hubs, including the spatial distribution of functional hubs, the hub-to-hub interconnectivity pattern, and the functional connections between hubs and non-hub brain regions, which we will refer to as “spokes.”

To address these questions, we investigated the development of functional brain network's hub architecture using resting-state functional magnetic resonance imaging (RS-fMRI; see Fig. 1 for an overview of data analyses). RS-fMRI measures intrinsic, high-amplitude, low-frequency blood-oxygen-level dependence signal (BOLD) fluctuations of the brain (Fox and Raichle 2007), and correlated RS-fMRI signal reflects functional connectivity independent of any particular brain state (Van Dijk et al. 2010). Correlation matrices among RS-fMRI time-series were calculated to quantify the strength of functional connectivity among brain regions, which were further analyzed using graph-theory algorithms. To reduce potential biases resulting from selecting a single spatial scale to define brain regions or a single edge threshold to define functional connections among brain regions, we explored graph-theoretic measures across different spatial scales and empirically optimized edge thresholds.

Figure 1.

Schematic illustration of analyses. We first constructed functional networks (steps AD) with 2 different spatial scales (B), and optimal thresholds were estimated and applied (D). Once networks were constructed, graph measures (E) were used to identify hub regions for both the voxel-wise network (H) and the functional ROI-based network (F). The importance of hubs was tested by simulating the impact of removing hubs from the ROI-based network (G). A list of potential hubs were identified from the voxel-wise network (steps HI), which the inter-connectivity pattern between functional hubs was analyzed (J), and functional connections between hubs and non-hub regions were examined (K). Results from steps E, H, J, and K were submitted to group analyses to test for age-related differences (L).

Figure 1.

Schematic illustration of analyses. We first constructed functional networks (steps AD) with 2 different spatial scales (B), and optimal thresholds were estimated and applied (D). Once networks were constructed, graph measures (E) were used to identify hub regions for both the voxel-wise network (H) and the functional ROI-based network (F). The importance of hubs was tested by simulating the impact of removing hubs from the ROI-based network (G). A list of potential hubs were identified from the voxel-wise network (steps HI), which the inter-connectivity pattern between functional hubs was analyzed (J), and functional connections between hubs and non-hub regions were examined (K). Results from steps E, H, J, and K were submitted to group analyses to test for age-related differences (L).

We utilized several graph-theoretic centrality measures (Rubinov and Sporns 2010) that are sensitive to 2 key hub properties: 1) Higher numbers of connections with other brain regions (“hub connectivity”, HC) and 2) higher volumes of information passing through a hub region relative to non-hub regions (“hub traffic”, HT). Centrality measures were consolidated into 2 simple hub scores using principal component analysis (PCA) to yield summary estimates of HC and HT. By testing a range of edge definition thresholds, spatial scales, and graph metrics, our results provided convergent information about the development of functional hubs.

Materials and Methods

Participants

One hundred and thirty-nine participants aged 10–20 years participated in the study in accordance with University of Pittsburgh Institutional Review Board guidelines. Participants and their first-degree relatives had no history of psychiatric disorders. Data from 40 participants were excluded due to excessive head motion, poor coverage in the cerebellum, or peripheral physiological equipment failure. We report data from 99 participants: 28 Children aged 10–12 years (M = 11.38, standard deviation [SD] = 0.78), 41 adolescents aged 13–17 years (M = 15.66, SD = 1.31), and 30 adults aged 18–20 years (M = 18.82, SD= 0.39).

Data Acquisition

Data were acquired at the University of Pittsburgh Medical Center Magnetic Resonance Research Center using a Siemens 3T Tim Trio (Erlangen, Germany). All participants spent ∼15 min in a mock scanner to acclimate them to the MR environment before entering the research scanner. Structural images were acquired using a sagittal magnetization-prepared rapid gradient-echo sequence (repetition time [TR] = 1570 ms, echo time [TE] = 3.04 ms, flip angle = 8°, inversion time [TI] = 800 ms, voxel size = 0.78125 × 0.78125 × 1 mm). Functional images were acquired using an echo-planar sequence sensitive to BOLD contrast (T2*; TR = 1.5 s, TE = 29 ms, flip angle = 70°, voxel size = 3.125 × 3.125 mm in-plane resolution, 29 contiguous 4-mm axial slices). We collected a 5-min (200 TRs) resting-state scan for each participant. During the resting-state scan, participants were asked to close their eyes, relax, but not fall asleep. Respiration and heartbeat were continuously recorded using a respiration belt and a pulse oximeter attached to the left index finger, which were later used to attenuate physiological noise.

MRI Data Preprocessing

Imaging data were preprocessed using AFNI (Cox 1996) and FSL (Smith et al. 2004). Freesurfer (Fischl et al. 2002) was used to segment gray matter, white mater, ventricle, and non-brain tissue (NBT) voxels. Information based on subjects' anatomical parcellations was used to remove non-neural noise from the resting-state BOLD signal and to restrict network analyses to gray matter voxels. Preprocessing steps included: 1) the Removal of sudden spikes caused by MR artifacts or large head motions, 2) slice-timing correction, 3) motion correction, 4) co-registration, 5) scaling each voxel time-series to a mean value of 100, and 6) linear detrending. We further utilized the AFNI's ANATICOR tool (Jo et al. 2010) to reduce hardware noise, draining vessel effect, and motion artifacts in each gray matter voxel via regression of the following nuisance variables: 1) 6-parameter motion regressors, 2) local white matter regressors averaged from white matter voxels within a spherical mask (radius = 30 mm) centered at each gray matter voxel of interest, 3) ventricle signal regressors, and 4) NBT regressors. Note that we did not include the global signal into the regression model to remove physiological noise. Instead we reduced the effect of physiological noise via modeling the effect of respiration and heart rhythms from the recorded physiological parameters, using a AFNI's RetroTS program (Glover et al. 2000; Birn et al. 2008). Controlling for physiological noise is particularly important for developmental studies because respiration and heart rate could systematically differ across age groups, potentially reflecting age-related differences in physiological reactivity to the MR environment. Time-series were then bandpass filtered at 0.009 Hz < f < 0.08 Hz and voxels were spatially smoothed within the gray matter mask (5-mm full-width at half-maximum Gaussian blur). Preprocessed data were spatially normalized to the Montreal Neurological Institute (MNI) 152 template using FSL's non-linear registration procedure, and resampled to 4 mm cubic voxels.

Network Analyses

Constructing Functional Networks

We constructed network matrices by calculating correlation matrices of resting-state BOLD signal fluctuations among network nodes. Correlation values were Fisher's z-transformed to improve normality. Given that the spatial resolution chosen to define network nodes (e.g. single voxel nodes, nodes based on anatomical landmarks, or nodes defined by functional regions of interests, ROIs) could strongly influence the results (Smith et al. 2011), we compared network metrics across different node definitions. In the present study, network “nodes” were defined at 2 spatial scales: 1) A voxel-wise approach in which network matrices were constructed for each participants' gray matter voxels (4 mm3 cubic voxels) and 2) a functional ROI approach consisting of 160 10-mm diameter functional brain regions defined based on a prior meta-analysis of fMRI activation studies (Dosenbach et al. 2010).

Voxel-wise sampling approach

Each participant's gray matter mask (derived from FreeSurfer) was spatially transformed into the MNI stereotactic space, and an intersection mask, consisting of 13 667 gray matter voxels present across participants, was created.

ROI sampling approach

We constructed networks using 160 functional ROIs identified from prior meta-analytic studies of cognitive control, error processing, default mode, memory, language, and sensorimotor functions (Dosenbach et al. 2010). This approach parcellated the brain into segregated functionally defined ROIs that covered most of the cortex.

Network Thresholding

Anatomical brain networks are sparse, having an estimated connection density ranging from 0.001 to 0.01 (Sporns et al. 2005). In contrast, functional brain networks derived from RS-fMRI data are densely connected before a proper threshold is applied to separate present versus absent functional connections among regions. Since there is no ideal biologically salient threshold that definitively identifies functional connections, the use of multiple thresholds are recommended to ensure that results reflect true network dynamics and not mathematical artifacts (Powers et al. 2010).

For graph metrics that required thresholding, we removed connections whose correlation coefficients fell below the threshold as well as all negative weights. To ensure that differences in graph metrics across age groups could not be attributed to differences in network density, we thresholded network matrices by edge density (i.e. the fraction of present connections relative to the total number of possible connections) across participants. For the ROI sampling approach, all graph analyses were performed at each density (D), ranging from D = 0.10 to 1.0. For voxel-wise analyses, because of the high computational burden, we used 3 different thresholding optimization approaches to derive a range of optimal thresholds: 1) Maximum clustering coefficient relative to a random graph (Elo et al. 2007), spectral clustering (Perkins and Langston 2009), and significance test of connection strength (Supplementary Methods). The estimated clustering coefficient density threshold was D = 0.008 and the spectral clustering threshold was D = 0.0004, consistent with previous estimates of the sparsity of human brain networks. When thresholded by connection weights that were significantly different from zero (at P < 0.001), the network density was 0.14. For graph metrics that can be calculated using edge weights (e.g. eigenvector centrality), unthresholded networks were used.

Graph-Theory Analyses

The separate analyses of network properties were conducted for the voxel-wise and functional ROI networks using graph-theoretic algorithms. Note that all networks were symmetric, undirected networks. Graph analyses were carried out using the C++ port of the Brain Connectivity Toolbox (http://code.google.com/p/bct-cpp/).

Definition and notations

A network (or graph) G can be described as G =(N,E), where N indicates the set of n of nodes in the network and E is the collection of connections (edges) between nodes. For binary networks, aij denotes the connection status between nodes i and j (aij= 1 if connection exists, otherwise aij = 0). For weighted networks, wij is the connection weight between nodes i and j.

Clustering coefficient

Clustering coefficient (C) is a measure of how nodes tend to cluster together to form connected local structures within a network, defined as 

formula
where Ci is the local clustering coefficient of node i, defined as 
formula
where Ei is the number of edges connected to node i, and ki is the number of neighbors of node i.

Path length

Path length measures, on average, how many steps one node must traverse to reach another node within the network. The path length (L) of a network is the average minimal distance between 2 nodes, defined as 

formula
where d(i,j) is the minimal distance (number of edges) between nodes i and j.

Efficiency

Efficiency (Eg) is the average inverse path length of the network. A network with larger path length is less efficient, as information will require more paths to reach its destination. Efficiency is defined as 

formula
where Eloci is the local efficiency of node i, defined as 
formula

Small-worldness

Small-world networks are characterized by high local clustering with short path length. Compared with random networks, small-world networks have higher clustering coefficients, but short path lengths comparable with those of random networks (Watts and Strogatz 1998). A single small-worldness index (Sw) can be calculated using 

formula
where small-world networks will have a Sw value greater than one.

Centrality measures

Graph-theoretic centrality measures estimate the functional significance of a node to a network, where different centrality measures are sensitive to different hub properties. We examined degree centrality, node strength centrality, and eigenvector centrality to estimate the degree of connectivity for each node (HC) and examined weighted and binary betweenness centrality to estimate the volume of information passing through each node (HT).

Degree centrality

For each node, degree centrality can be calculated by summing the number of edges connected to each node 

formula

Node strength

Node strength is defined as the sum of connection weights connected to each node 

formula

Eigenvector centrality

Eigenvector centrality weights were determined by the number of edges among nodes (similar to degree) and the connectedness of neighboring nodes. In this way, the eigenvector centrality estimate for a given node is proportional to the sum of the centrality estimates of neighboring nodes. This metric uses the full information of unthresholded weighted networks (after removing negative weights) and is defined by the eigenvector of the largest eigenvalue of a symmetric matrix 

formula
where λ is the largest eigenvalue, and e is the corresponding eigenvector.

Betweenness centrality

Betweenness centrality (bi) is defined as the fraction of the shortest paths between any pair of nodes that travel through the node 

formula
where ghj is the number of shortest path between node h and j, and ghj(i) is the number of shortest paths between h and j that pass through i. When calculating weighted betweenness centrality, weighted distance (the inverse of the correlation coefficient) was used to calculate the shortest path length.

Hub connectivity and hub traffic scores—aggregated centrality measures

PCA was used to consolidate the centrality measures across different thresholds into 2 “hub scores” reflecting the largest proportion of shared variance among measures. Given large differences in scaling, centrality measures were first normalized to a zero mean and unit variance. Degree centrality from all thresholds, node strength, and eigenvector centrality were consolidated into the HC score, while binary betweenness centrality (all thresholds) and weighted betweenness centrality were consolidated into the HT score. These data were then entered as variables into a single data matrix representing all participants and nodes (voxels or ROIs), with participants and nodes concatenated across rows (i.e. as observations). Principal components were generated by singular value decomposition, and the first principal component score of each node for each participant was retained as hub scores. The use of a single data matrix for dimensionality reduction ensured that the scale of hub scores was comparable across participants and nodes. For voxel-wise analyses, the HC accounted for 72% of variance among degree centrality, node strength, and eigenvector centrality, while the HT accounted for 83% of the variance for betweenness centrality measures. For ROI analyses, the HC accounted for 74% of the variance, and the HT accounted for 91%. Thus, the PCA-derived centrality measures summarize the metrics well and provide summary estimates that are not biased by the selection of edge threshold.

Identifying candidate hub regions

From the whole-brain voxel-wise analysis, a conjunction map of the top 5% of voxels from the 2 hub scores was created, and an automated (local extrema) search algorithm was used to locate peaks of the HT score within this map. Peaks separated by <10 mm were consolidated by averaging coordinates. Peak coordinates were then referenced with atlases included with FSL to determine each candidate hub's anatomical label and Brodmann area assignment (Supplementary Methods).

Hub-to-hub network

After identifying the coordinates of candidate hubs, we generated 94 hub ROIs using spheres with an 8 mm radius around the peak coordinates. The BOLD time-series were averaged within each hub ROI and a 94 × 94 cross-correlation matrix was calculated. Graph metrics were then calculated across density thresholds (D= 0.1–1). The strengths of connections between functional hubs were also examined across age.

Functional Networks Associated with Functional Hubs

We performed voxel-wise seed-based functional connectivity analyses that correlated mean time-series from each cortical hub with voxels across the whole brain. This procedure generated one voxel-wise connectivity map for each cortical hub, representing the strength of connections between the cortical hub and non-hub “spoke” regions. These images were then submitted to group analyses.

Developmental Differences

Age was treated as a continuous variable in regression models, where both linear and quadratic age effects were examined. In addition age was also treated as a categorical variable to compare graph metrics and connectivity strength between childhood (aged 8–12), adolescence (aged 13–17), and adults (aged 18 or older). No standards exist for defining child and adolescent age ranges, thus groups were defined based on our past behavioral studies indicating differential cognitive performance (Luna et al. 2004). The directions of group differences were analyzed using an independent sample t-test to compare children versus adolescents and adolescents versus adults.

Multiple Comparisons Correction

To explore the effects of age on voxel-wise hub scores, we used a liberal uncorrected voxel threshold of P < 0.05 and a minimum cluster size of 10 contiguous voxels. For the effects of age on ROI hub scores, graph metrics, and the strength of hub-to-hub connections, we corrected for family-wise error using the Holm-Bonferroni method (Holm 1979). For the voxel-wise hub–spoke analyses, we controlled for multiple comparisons by correcting for the number of tests within statistical parametric maps (i.e. voxelwise) and between maps, where 2 maps (adults vs. adolescents and adolescents vs. children) were generated for each hub. The minimum cluster size of contiguous significant voxels was determined via Monte Carlo simulation using AFNI's AlphaSim. Spatial smoothness was determined empirically from the residual time-series of each hub–spoke functional connectivity analysis, and averaged across subjects (full-width at half-maximum x = 4.08 mm, y = 4.06 mm, z = 6.23 mm). Note that during preprocessing we only performed spatial smoothing within gray matter voxels. Cluster thresholds were Bonferroni corrected to further control for the large number of comparisons (94 hubs × 2 age group comparisons = 198) that were performed (α = 0.05/198). A voxel threshold of P < 0.01 and a minimum cluster size of 18 contiguous voxels (1152 mm3) yielded a corrected cluster threshold of P< 0.00025, thereby maintaining an overall α level of 0.05 across voxel-wise maps for all hubs. The center coordinates of corrected clusters were reported.

Bayesian Analysis of the Association Between Age and Hub Scores

To obtain an estimate of the strength of the association between age and hub HC and HT scores, we conducted a Bayesian correlation analysis for each of the 160 ROIs. The advantage of the Bayesian approach over frequentist statistics, in this case, is that one can derive a posterior probability that the age–hub score association is weak, here |r| < 0.2, rather than interpreting a non-significant hypothesis test, which may be attributable to a weak effect, insufficient power, or between-subjects variability (Morey and Rouder 2011).

For each ROI, a Bayesian linear regression model was estimated using a Gibbs sampler to simulate draws from the posterior distribution of the regression of hubs scores on age, where 5000 iterations were discarded for Markov Chain Monte Carlo (MCMC) burn-in and 500 000 values were drawn from the sampler with a thinning interval of 5, which yielded a posterior distribution consisting of 100 000 correlation values. The prior distribution for regression coefficients was assumed to be Gaussian, whereas the prior for conditional error variance was an inverse Gamma distribution. Models were estimated using the MCMCpack library for R, and trace plots indicated good convergence of the MCMC chain toward stationarity. The posterior probability of a weak correlation was computed as the number of draws from the posterior distribution that fell in the range −0.2 < r < 0.2. This correlation range was chosen because it corresponds to a small-to-medium effect according to established guidelines for associational effect sizes (Cohen 1988).

Head Movement Analysis

To mitigate the effects of head motion on estimates of functional connectivity, we first excluded subjects with excessive head movement (average root mean square, RMS >1 mm or 1°). Secondly, as described below, before calculating correlation estimates, we censored volumes within each retained subject's fMRI time-series that were associated with sudden head motion. The measures of head movement were obtained from AFNI's rigid-body rotation and translation realignment algorithm (3dvolreg; Cox 1996). Translations and rotations in the x, y, and z dimensions were averaged across volumes, and total RMS of linear and angular precision measures was calculated. Participants whose average RMS movement exceeded 1 mm or 1° were discarded (n = 17). One-way analysis of variance (ANOVA) did not identify a significant difference of head motion among age groups, F2, 96 = 1.94, P = 0.15, and average movement for all age groups was considerably below the cutoff of 1 mm or 1° (average RMS = 0.23 for adults, 0.29 adolescents, 0.33 for children).

For subjects included for analyses, we calculated 2 data quality metrics proposed by Power et al. (2012): Frame-wise displacement (FD) and the RMS variance of the temporal derivative of fMRI time-series (DVARS). FD summarizes instantaneous head motion, whereas DVARS measures the average signal change across all voxels from one volume to the next (large head motion is typically associated with large-amplitude signal fluctuations). These metrics were used to identify which fMRI volumes should be censored from data analyses due to excessive head motion artifact at those timepoints. Adopting the threshold suggested by Power et al., we censored volumes from each subject's time-series where FD >0.5 mm and DVARS >0.5 (0.5% signal change). Overall, 7% of volumes were censored for children, 2% of volumes for adolescents, and 0.2% for adults. Two correlation matrices were then computed for each subject: One where connectivity estimates did not include censored frames and the other where all time points in the fMRI data were used to compute connectivity. We then compared centrality measures of all ROIs before and after censoring.

Results

Centrality Measures Identified a Set of Regions As Functional Hubs

Voxel-Wise Analyses

We first mapped brain regions that displayed hub properties (HC and HT) using a whole-brain voxel-wise approach (Fig. 2). The voxel-wise analyses sampled resting-state functional connections in the entire brain at high spatial resolution without predefined anatomical or functional parcellations, allowing us to identify potential hub regions in an unbiased and data-driven manner. Although the spatial distribution of regions with high HC versus HT largely overlapped, there were also important differences. Across all age groups, regions displaying high HC were primarily located in midline structures including the precuneus, the cingulate gyrus, the medial prefrontal cortex, and unimodal sensory cortices, including visual, motor, and auditory cortices. In contrast, regions showing high HT were located primarily in associative frontal-parietal cortices, and in subcortical structures including the thalamus and basal ganglia. Regions in the anterior temporal pole, the superior temporal cortex, frontal and posterior midline structures, and the cerebellum exhibited both high HT and HC.

Figure 2.

Spatial distribution of functional hubs for each age group. Voxel-wise mapping of brain regions that displayed high HC and high HT. Voxels were thresholded at the 95% percentile or higher for each hub score.

Figure 2.

Spatial distribution of functional hubs for each age group. Voxel-wise mapping of brain regions that displayed high HC and high HT. Voxels were thresholded at the 95% percentile or higher for each hub score.

The spatial separability of HC versus HT properties coincided with lower correlations between HC and HT scores relative to correlations among graph-theory metrics comprising each of these summary indices (Fig. 3, Supplementary Table S1). We further assessed the concordance of centrality measures using Kendall's W (Supplementary Fig. S1; Legendre 2005; Zuo et al. 2011). Across all voxels, concordance among all centrality measures (median W = 0.62) was significantly lower than concordances among measures comprising the HC score (median W= 0.89; Wilcoxon signed-rank test P < 0.001) and the HT score (median W = 0.81; Wilcoxon signed-rank test P < 0.001). Further, the spatial patterns of hub scores compared with their constituent centrality measures were highly similar (Supplementary Figs S2 and S3). Altogether, these results corroborate the validity of using PCA to summarize node centrality properties. Centrality measures across different density thresholds were also highly correlated (Fig. 3, Supplementary Table S1), further suggesting that the sparse networks derived from our threshold optimization procedures still faithfully captured functional brain network's architecture and that our results were not sensitive to edge threshold definition.

Figure 3.

Correlation between hub scores and individual centrality measures. Lower left blue box encircles correlations between centrality measures constitute the HC score. Upper right blue box encircles correlations between centrality measures constitute the HT score. D, density.

Figure 3.

Correlation between hub scores and individual centrality measures. Lower left blue box encircles correlations between centrality measures constitute the HC score. Upper right blue box encircles correlations between centrality measures constitute the HT score. D, density.

After identifying voxels across participants that ranked in the 95th percentile or greater of either HC or HT score distributions, we used a peak extraction procedure (see Materials and Methods) to identify 94 regions as putative functional hubs (Supplementary Table S2). The most prominent hubs included the precuneus, the right superior parietal cortex, bilateral anterior temporal pole, bilateral medial dorsal thalamus, bilateral superior temporal gyrus, medial prefrontal cortex, bilateral lingual gyrus/visual cortex, bilateral caudate, left middle frontal cortex, and bilateral cerebellum (Table 1).

Table 1

Putative hubs that showed strong HT or HC scores

MNI coordinates
 
Region BA Hub traffic score Hub connectivity score 
x y z Hemisphere 
Top ranking hubs identified from voxel-wise analysis, ranked by the HT score 
−40 50 Right Precuneus 3.71 2.11 
28 −76 46 Right Superior parietal 3.23 2.07 
16 −80 42 Right Superior parietal 3.06 2.85 
52 16 −10 Right Anterior temporal pole 38 2.96 1.60 
−14 Right Thalamus (medial dorsal)  2.85 0.98 
−8 −14 Left Thalamus (medial dorsal)  2.82 0.85 
−4 −60 Left Lingual 17 2.80 2.57 
−8 Left Caudate  2.80 1.23 
−52 16 −14 Left Anterior temporal pole 38 2.76 2.19 
52 16 26 Left Middle frontal 2.69 1.35 
−56 Right Lingual 18 2.66 3.38 
−76 34 Right Precuneus 2.60 2.77 
12 10 Right Caudate  2.37 1.01 
48 −64 −26 Right Cerebellum (Crus I)  2.26 2.02 
−32 −80 −30 Left Cerebellum (Crus I)  2.23 2.33 
64 −8 Right Superior temporal 41 2.10 2.68 
−64 −16 Left Superior temporal 41 2.10 2.06 
−16 54  Medial frontal 2.10 1.92 
High connectivity hubs identified from ROI analyses 
−40 50 Right Precuneus 14.27 15.62 
−54 −22 Left Superior temporal 41 13.61 11.52 
−1 52  Supplementary motor 13.14 12.13 
17 −68 20 Right Cuneus 17 11.16 3.38 
15 −77 32 Right Precuneus 11.11 2.35 
19 −66 −1 Right Lingual 19 10.66 2.76 
15 45 Right Medial frontal 10.60 11.82 
14 −75 −21 Right Cerebellum (VI)  10.24 2.73 
High traffic hubs identified from ROI analyses 
−40 50 Right Precuneus 14.27 15.62 
−1 52  Supplementary motor 13.14 12.13 
−3 −38 45 Left Precuneus 10.07 12.09 
15 45  Medial frontal 10.60 11.82 
−54 −22 Left Superior temporal 41 13.61 11.52 
−1 28 40 Left Anterior cingulate 32 7.42 11.32 
−32 −58 46 Left Inferior parietal 4.48 9.12 
−29 57 10 Left Middle frontal 10 6.79 9.00 
MNI coordinates
 
Region BA Hub traffic score Hub connectivity score 
x y z Hemisphere 
Top ranking hubs identified from voxel-wise analysis, ranked by the HT score 
−40 50 Right Precuneus 3.71 2.11 
28 −76 46 Right Superior parietal 3.23 2.07 
16 −80 42 Right Superior parietal 3.06 2.85 
52 16 −10 Right Anterior temporal pole 38 2.96 1.60 
−14 Right Thalamus (medial dorsal)  2.85 0.98 
−8 −14 Left Thalamus (medial dorsal)  2.82 0.85 
−4 −60 Left Lingual 17 2.80 2.57 
−8 Left Caudate  2.80 1.23 
−52 16 −14 Left Anterior temporal pole 38 2.76 2.19 
52 16 26 Left Middle frontal 2.69 1.35 
−56 Right Lingual 18 2.66 3.38 
−76 34 Right Precuneus 2.60 2.77 
12 10 Right Caudate  2.37 1.01 
48 −64 −26 Right Cerebellum (Crus I)  2.26 2.02 
−32 −80 −30 Left Cerebellum (Crus I)  2.23 2.33 
64 −8 Right Superior temporal 41 2.10 2.68 
−64 −16 Left Superior temporal 41 2.10 2.06 
−16 54  Medial frontal 2.10 1.92 
High connectivity hubs identified from ROI analyses 
−40 50 Right Precuneus 14.27 15.62 
−54 −22 Left Superior temporal 41 13.61 11.52 
−1 52  Supplementary motor 13.14 12.13 
17 −68 20 Right Cuneus 17 11.16 3.38 
15 −77 32 Right Precuneus 11.11 2.35 
19 −66 −1 Right Lingual 19 10.66 2.76 
15 45 Right Medial frontal 10.60 11.82 
14 −75 −21 Right Cerebellum (VI)  10.24 2.73 
High traffic hubs identified from ROI analyses 
−40 50 Right Precuneus 14.27 15.62 
−1 52  Supplementary motor 13.14 12.13 
−3 −38 45 Left Precuneus 10.07 12.09 
15 45  Medial frontal 10.60 11.82 
−54 −22 Left Superior temporal 41 13.61 11.52 
−1 28 40 Left Anterior cingulate 32 7.42 11.32 
−32 −58 46 Left Inferior parietal 4.48 9.12 
−29 57 10 Left Middle frontal 10 6.79 9.00 

BA, Brodmann areas.

ROI Analyses

Given that the spatial scale one uses to define network nodes could influence graph analyses (Power et al. 2010), we repeated the above analyses using functional brain regions of a lower spatial resolution, using a partition scheme that consisted of 160 functionally defined ROIs (10-mm diameter spheres, ROIs list in Supplementary Table S3) associated with various cognitive functions derived from a previous meta-analysis (Dosenbach et al. 2010). Across all subjects, we identified ROIs (n = 8) that showed the strongest HC or HT (95th percentile or greater; Table 1). Consistent with the voxel-wise analyses, ROIs in the precuneus, the superior temporal gyrus, and the medial prefrontal cortex showed both high HC and HT. ROIs in unimodal sensory cortices including the visual cortices and bilateral superior temporal gyrus had high HC, whereas ROIs in frontal-parietal associative cortices including anterior cingulate, the left middle prefrontal cortex, and the left inferior parietal cortex showed high HT.

Hub Architecture is Stable Across Development

Voxel-Wise Analyses

As depicted in Figure 2, the spatial distribution of regions with high HT and HC was highly consistent across age groups. To formally test developmental differences in HT and HC, we treated age as a grouping variable (children, adolescents, and adults) in ANOVA and as a continuous variable in regression models. Normality of HT and HC scores was verified (Supplementary Fig. S4) prior to parametric tests, and both hub scores were approximately normally distributed (HC skewness = 0.26, excess kurtosis = 0.02; HT skewness = 0.42, excess kurtosis = 0.25; Glass et al. 1972). Using a liberal threshold (uncorrected voxel threshold P < 0.05, minimal cluster size = 10 contiguous voxels), both ANOVA and regression analyses revealed no significant age effects in the location of functional hubs, the mean strength of HC and HT (Fig. 4A), and in all individual centrality measures.

Figure 4.

HC scores and HT scores showed no developmental differences across age. (A) Regression of HT and HC scores by age from 4 representative functional hubs identified from the voxel-wise analyses. For each subject, hub scores were averaged within an 8-mm hub ROI centered on the identified coordinate (see Materials and Methods). No significant (P < 0.05 corrected) age effect was found for these 4 functional hubs and all other hub regions that are not shown in this figure. ANOVA analyses also did not find any significant age differences in hub scores. (B) Regression of HT and HC scores by age from 3 example ROIs identified from the ROI analyses (Table 1). No significant (P < 0.05 corrected) age effect was found for these ROIs and all other ROIs not shown in this figure. ANOVA analyses also did not find any significant age differences.

Figure 4.

HC scores and HT scores showed no developmental differences across age. (A) Regression of HT and HC scores by age from 4 representative functional hubs identified from the voxel-wise analyses. For each subject, hub scores were averaged within an 8-mm hub ROI centered on the identified coordinate (see Materials and Methods). No significant (P < 0.05 corrected) age effect was found for these 4 functional hubs and all other hub regions that are not shown in this figure. ANOVA analyses also did not find any significant age differences in hub scores. (B) Regression of HT and HC scores by age from 3 example ROIs identified from the ROI analyses (Table 1). No significant (P < 0.05 corrected) age effect was found for these ROIs and all other ROIs not shown in this figure. ANOVA analyses also did not find any significant age differences.

ROI Analyses

Consistent with the voxel-wise analyses, individual centrality measures and hub scores of all ROIs did not change with age (P > 0.05, Holm-Bonferroni corrected, Fig. 4B).

Effect Size Estimates of Developmental Differences

To ensure that the lack of age effects was not attributable to low statistical power, we estimated voxel-wise effect sizes of HC and HT scores as a function of age group (children, adolescents, and adults). The average standardized mean difference for adults versus adolescents contrast for HC was Cohen's d = 0.21, SD = 0.17; HT d = 0.2, SD= 0.17. The effect size of the adolescents versus children contrast for HC was d = 0.22, SD = 0.18; HT d = 0.21, SD = 0.17. Effect sizes in this range correspond to an overlap of approximately 85% in the distribution of scores between groups, which is a “small” effect size (Cohen 1988). Thus, our effect size analyses indicated that any age-related differences in hub properties are very subtle, if not trivial.

Bayesian Analysis of Developmental Differences

The above results suggested that age-related changes in hub properties were minimal, but a non-significant hypothesis test cannot be used directly to argue for a weak or non-existent effect. Thus, to augment our effect size analyses above, we also correlated age with hub scores for each ROI using Bayesian linear regression models. Across 160 ROIs, the median posterior probabilities that the correlations between age and HC and age and HT were weak, −0.2 < r < 0.2, were 0.9 and 0.86, respectively. For the vast majority of ROIs, the 95% credible intervals of these correlations included zero, supporting the conclusion that age–hub score associations were weak or absent in most cases (Supplementary Figs S5 and S6), and across a population of correlation values, one would expect a few significant correlations by chance.

Functional Hubs Support Functional Integration

To test whether hubs play a larger role in information integration than non-hub regions, we simulated the effects of “lesioning” functional hubs on the brain network's topology (Sporns et al. 2007). Specifically, we first removed all connections to the 12 hubs identified in the aforementioned ROI analyses (Table 1) and then recalculated the global efficiency and clustering coefficient of the network. Efficiency is a function of the average shortest distance among nodes, and shorter path lengths yield greater efficiency of information transfer. The clustering coefficient summarizes the ability of the network to integrate information within tightly connected clusters. To produce confidence intervals on the network estimates following hub lesioning, we removed all connections to 12 randomly selected non-hub ROIs and re-estimated efficiency and clustering coefficient. This procedure was repeated 1000 times in order to generate an empirical cumulative distribution function amenable to statistical testing. We found that at all ages, lesioning functional hubs significantly reduced efficiency and clustering coefficients, regardless of network density, compared with lesioning non-hub regions (Fig. 5). These results corroborate the role of functional hubs in facilitating information integration.

Figure 5.

Functional hubs are critical for information integration. Efficiency and clustering coefficient measure was calculated for each individual subject, and averaged within age group. To estimate confidence interval, a null network was constructed by randomly leisioning connections to 12 randomly selected ROIs 1000 times. At each density threshold, upper and lower confidence interval values were determined by finding the 5th and 95th percentile values. The impact of removing functional hubs was examined across different density thresholds (X axis). Across network density thresholds, the efficiency and clustering coefficient of children, adolescents, and adults' brain network were reduced when connections to functional hubs were removed.

Figure 5.

Functional hubs are critical for information integration. Efficiency and clustering coefficient measure was calculated for each individual subject, and averaged within age group. To estimate confidence interval, a null network was constructed by randomly leisioning connections to 12 randomly selected ROIs 1000 times. At each density threshold, upper and lower confidence interval values were determined by finding the 5th and 95th percentile values. The impact of removing functional hubs was examined across different density thresholds (X axis). Across network density thresholds, the efficiency and clustering coefficient of children, adolescents, and adults' brain network were reduced when connections to functional hubs were removed.

Functional Hubs Form a Stable, Efficient Small-World Network Across Development

From voxel-wise analyses, we identified 94 regions as putative functional hubs (Supplementary Table S2). We further analyzed how hubs communicate with other hubs by examining the network topology of the hub-to-hub network. We found that this core hub-to-hub network displayed efficient small-world properties from ages 10 to 20 (Fig. 6), as indicated by high clustering coefficient and low path length, as well as a small-worldness index greater than one. We then analyzed age-related differences in path length, clustering coefficient, efficiency, and small-worldness in the hub-to-hub network and tested whether the strength of connections among functional hubs changed over development. We found that the strength of functional connections among hubs did not change with age (P > 0.05 Holm-Bonferroni corrected), nor did the path length, clustering coefficient, small-worldness, or efficiency of the hub-to-hub network (Fig. 6; P > 0.05 uncorrected). These results suggest that functional hubs form a small-world network that is stable from ages 10 to 20.

Figure 6.

Network topology of the core hub-to-hub network. Graph metrics of the identified core hub-to-hub network (see Material and Methods) were examined across development and density thresholds (X axis). Clustering coefficient, path length, efficiency, and small-worldness of the hub-to-hub network were highly stable across development, showing highly similar pattern across age groups. This core hub-to-hub network showed high clustering and low path length. Small-worldness measure was greater than one across all thresholds and age groups, indicating networks having small-world architecture.

Figure 6.

Network topology of the core hub-to-hub network. Graph metrics of the identified core hub-to-hub network (see Material and Methods) were examined across development and density thresholds (X axis). Clustering coefficient, path length, efficiency, and small-worldness of the hub-to-hub network were highly stable across development, showing highly similar pattern across age groups. This core hub-to-hub network showed high clustering and low path length. Small-worldness measure was greater than one across all thresholds and age groups, indicating networks having small-world architecture.

Strength of Connections Between Functional Hubs and Distributed Non-Hub Brain Regions (Spokes) Changes with Age

To examine the development of functional connections between functional hubs and non-hub spoke regions, we constructed whole-brain functional connectivity maps for each cortical hub. These functional connectivity maps were then tested with age as a group variable using planned t-tests to characterize developmental transitions across specific stages. Twenty-six of the 94 functional hubs showed age-related increases or decreases in the functional connectivity strength of “hub–spoke” connections (Table 2). Note that the predominant direction of information flow cannot be discerned because connectivity was calculated using the Pearson correlation coefficient. The majority of the developmental changes occurred during the transition from childhood to adolescence, and fewer age-related differences were evident from adolescence to adulthood (Fig. 7).

Table 2

List of connections between functional hubs and spokes that showed developmental change

Hub
 
Target non-hub region
 
x y z  BA Label x y z  BA Label 
Developmental pattern C > T 
−32 −80 −30 Left  Cerebellum (Crus I) 16 −88 38 Right 19 Superior occipital/cuneus 
−4 −68 54 Left Precuneus 44 −68 40 Right 39 Inferior parietal 
−56 Right 10 Lingual −21 −49 −5 Left 19 Lingual 
Developmental pattern T > C 
48 30 Right Superior medial frontal −30 −76 −39 Left  Cerebellum (Crus II) 
      31 −75 −43 Right  Cerebellum (Crus II) 
−4 60 30 Left 10 Medial frontal −44 −74 28 Left 22 Middle temporal 
48 14 Right 10 Medial frontal 16 −82 −39 Right  Cerebellum 
46 Right SMA 60 30 Right Precentral 
      48 −60 −14 Right 37 Inferior temporal 
16 38 Right Middle cingulate −32 22 24 Left 44 Inferior frontal 
      35 15 −14 Right 47 Insula 
−24 40 42 Left Superior frontal 42 −77 26 Right 39 Middle temporal 
      67 Right 10 Medial frontal 
−28 60 Left 10 Middle frontal 14 36  32 Middle cingulate 
52 16 26 Right 44/45 Inferior frontal 46 −40 51 Right 40 Inferior parietal 
      −42 −41 45 Left 40 Inferior parietal 
      58 −58 −9 Right 37 Inferior temporal 
−20 −22 Left 47 Orbital frontal 38 14 Right 13 Insula 
      −32 −71 −34 Left  Cerebellum (Crus II) 
48 −8 54 Right Precentral 57 −5 32 Right Precentral 
      −32 −13 −8 Left  Putamen 
56 −24 50 Right Postcentral 51 10 Right 44 Inferior frontal 
48 −36 54 Right 40 Inferior parietal 46 38 17 Right 45 Inferior frontal 
      53 13 16 Right 44 Inferior frontal 
−52 −32 50 Left 40 Inferior parietal 51 11 Right 44 Inferior frontal 
64 −28 30 Right 40 Inferior parietal −34 −60 −34 Left  Cerebellum (Crus II) 
      −34 21 Left 13 Insula 
52 16 −10 Right 38 Temporal pole 56 −17 36 Right Postcentral 
      −21 30  23 Middle cingulate 
      −58 −35 39 Left 40 Inferior parietal 
      40 10 −4 Right 13 Insula 
      −39 −10 Left 13 Insula 
−52 16 −14 Left 38 Temporal pole −34 38 27 Left 10 Middle frontal 
      −39 43 15 Left 10 Middle frontal 
−64 −16 Left 22 Superior temporal 56 −9 44 Right Precentral 
−28 −8 −2 Left  Putamen −32 32 Left Precentral 
48 −52 −34 Right  Cerebellum (Crus I) −2 46 13 Left 32 Anterior cingulate 
Developmental pattern T > A 
−28 −8 −2 Left  Putamen −35 −40 64 Left Postcentral 
      −58 −59 10 Left 22 Middle temporal 
      51 −64 −13 Right 19 Inferior occipital 
28 −72 −34 Right  Cerebellum (Crus I) −32 23 45 Left Middle frontal 
      37 20 46 Right Middle frontal 
      −31 −23 61 Left Central sulculs/precentral gyrus 
Developmental pattern A > T 
−44 −68 −30 Left  Cerebellum (Crus I) 32 −84 30 Right 19 Middle occipital 
      49 −42 Right 22 Middle temporal 
48 −64 −26 Right  Cerebellum (Crus I) −25 −63 −4 Left 17 Lingual 
Hub
 
Target non-hub region
 
x y z  BA Label x y z  BA Label 
Developmental pattern C > T 
−32 −80 −30 Left  Cerebellum (Crus I) 16 −88 38 Right 19 Superior occipital/cuneus 
−4 −68 54 Left Precuneus 44 −68 40 Right 39 Inferior parietal 
−56 Right 10 Lingual −21 −49 −5 Left 19 Lingual 
Developmental pattern T > C 
48 30 Right Superior medial frontal −30 −76 −39 Left  Cerebellum (Crus II) 
      31 −75 −43 Right  Cerebellum (Crus II) 
−4 60 30 Left 10 Medial frontal −44 −74 28 Left 22 Middle temporal 
48 14 Right 10 Medial frontal 16 −82 −39 Right  Cerebellum 
46 Right SMA 60 30 Right Precentral 
      48 −60 −14 Right 37 Inferior temporal 
16 38 Right Middle cingulate −32 22 24 Left 44 Inferior frontal 
      35 15 −14 Right 47 Insula 
−24 40 42 Left Superior frontal 42 −77 26 Right 39 Middle temporal 
      67 Right 10 Medial frontal 
−28 60 Left 10 Middle frontal 14 36  32 Middle cingulate 
52 16 26 Right 44/45 Inferior frontal 46 −40 51 Right 40 Inferior parietal 
      −42 −41 45 Left 40 Inferior parietal 
      58 −58 −9 Right 37 Inferior temporal 
−20 −22 Left 47 Orbital frontal 38 14 Right 13 Insula 
      −32 −71 −34 Left  Cerebellum (Crus II) 
48 −8 54 Right Precentral 57 −5 32 Right Precentral 
      −32 −13 −8 Left  Putamen 
56 −24 50 Right Postcentral 51 10 Right 44 Inferior frontal 
48 −36 54 Right 40 Inferior parietal 46 38 17 Right 45 Inferior frontal 
      53 13 16 Right 44 Inferior frontal 
−52 −32 50 Left 40 Inferior parietal 51 11 Right 44 Inferior frontal 
64 −28 30 Right 40 Inferior parietal −34 −60 −34 Left  Cerebellum (Crus II) 
      −34 21 Left 13 Insula 
52 16 −10 Right 38 Temporal pole 56 −17 36 Right Postcentral 
      −21 30  23 Middle cingulate 
      −58 −35 39 Left 40 Inferior parietal 
      40 10 −4 Right 13 Insula 
      −39 −10 Left 13 Insula 
−52 16 −14 Left 38 Temporal pole −34 38 27 Left 10 Middle frontal 
      −39 43 15 Left 10 Middle frontal 
−64 −16 Left 22 Superior temporal 56 −9 44 Right Precentral 
−28 −8 −2 Left  Putamen −32 32 Left Precentral 
48 −52 −34 Right  Cerebellum (Crus I) −2 46 13 Left 32 Anterior cingulate 
Developmental pattern T > A 
−28 −8 −2 Left  Putamen −35 −40 64 Left Postcentral 
      −58 −59 10 Left 22 Middle temporal 
      51 −64 −13 Right 19 Inferior occipital 
28 −72 −34 Right  Cerebellum (Crus I) −32 23 45 Left Middle frontal 
      37 20 46 Right Middle frontal 
      −31 −23 61 Left Central sulculs/precentral gyrus 
Developmental pattern A > T 
−44 −68 −30 Left  Cerebellum (Crus I) 32 −84 30 Right 19 Middle occipital 
      49 −42 Right 22 Middle temporal 
48 −64 −26 Right  Cerebellum (Crus I) −25 −63 −4 Left 17 Lingual 

Coordinates are in MNI space. BA, Brodmann areas; SMA, supplementary motor areas.

Figure 7.

Developmental changes in connections connecting functional hubs and distributed regions. Hubs and non-hub regions, along with the interconnections that were significantly modulated by age are placed on 2 brains depicting the sagittal (left) and axial plane (right), where either the X- or Z-axes were collapsed. Placing of hubs and non-hub regions were determined by the MNI coordinate. Coordinates of each region listed in Table 2, and legends of hubs versus non-hub regions on the right. Blue lines marked strength of connections that decreased with age (adolescents < children or adults < adolescents), red lines marked connections increased with age (adolescents > children, adults > adolescents). (A) Age-related changes in connections connecting hubs and non-hub regions from childhood to adolescence. (B) Age-related changes in connections connecting hubs and non-hub regions from adolescence to adulthood. Note that all cerebellar regions were located in Crus I and Crus II. A, anterior; P, posterior; I, Inferior; S, superior.

Figure 7.

Developmental changes in connections connecting functional hubs and distributed regions. Hubs and non-hub regions, along with the interconnections that were significantly modulated by age are placed on 2 brains depicting the sagittal (left) and axial plane (right), where either the X- or Z-axes were collapsed. Placing of hubs and non-hub regions were determined by the MNI coordinate. Coordinates of each region listed in Table 2, and legends of hubs versus non-hub regions on the right. Blue lines marked strength of connections that decreased with age (adolescents < children or adults < adolescents), red lines marked connections increased with age (adolescents > children, adults > adolescents). (A) Age-related changes in connections connecting hubs and non-hub regions from childhood to adolescence. (B) Age-related changes in connections connecting hubs and non-hub regions from adolescence to adulthood. Note that all cerebellar regions were located in Crus I and Crus II. A, anterior; P, posterior; I, Inferior; S, superior.

From childhood to adolescence, we found developmental decreases in functional connectivity in the posterior part of the brain, connecting hubs in the cerebellum, the precuneus, and visual cortex with spoke regions in the occipital and parietal cortices. Developmental increases from childhood to adolescence were more widespread, primarily including cortical–cortical and intra-cortical connections. The frontal cortex was heavily represented in age-related connectivity increases in hub–spoke connections, particularly the middle frontal gyrus, the inferior frontal gyrus, the insula, the medial frontal cortex, and the supplementary motor area (Fig. 7A). For example, the strengths of frontal–parietal, frontal–temporal, frontal–cerebellar, and temporal–frontal connections increased with age (Table 2). From adolescence to adulthood, only connections between cerebellar hubs, occipital, and temporal spoke regions increased (Table 2, Fig. 7B), whereas connections between cerebellar hubs and frontal spoke regions, and between the putamen hub and occipital, temporal and parietal spoke regions decreased. Hub and spoke cerebellar regions were located in the lateral cerebellum, Crus I and II.

Control Analyses: Effects of Head Motion

Two recent studies found that head movement during fMRI acquisition introduces systematic artifacts into RS-fMRI connectivity estimates that cannot be adequately controlled for by conventional signal processing methods (Power et al. 2012; Van Dijk et al. 2012). To ensure that our results were not unduly influenced by motion artifact, we calculated FD and DVARS to identify data volumes that were probably contaminated by motion artifact, censored these volumes, and repeated all ROI analyses. For HC scores, the range of correlations for all ROIs (across subjects) calculated from connectivity estimates with and without censoring was 0.9997–1.0. The range of correlations between all ROIs' HC scores before and after censoring was 0.9959–1.0. Together, these results suggest that our results are likely not attributable to the effects of head motion.

Discussion

A Combination of Different Centrality Measures Identified Sets of Brain Regions as Functional Hubs That Support Information Integration

We sought to provide a more comprehensive characterization of functional hub development using a plurality of centrality measures, different spatial scales, and several plausible definitions of network edges. Most past studies have focused on one centrality measure to identify functional hubs, yet different centrality measures may yield discrepant results (Zuo et al. 2011) due to their varying sensitivities to HC and HT. For instance, it is possible that some hubs are densely connected within a sub-network or a local cluster but do not mediate communication among sub-networks (i.e. high HC but low HT). Conversely, hubs could be placed along paths connecting distant portions of the network, thereby supporting a high volume of information traffic, but may have a lower degree of connectivity (i.e. high HT but lower HC). Indeed, we found that hubs in the midline, occipital, temporal, and motor regions showed a higher degree of connectivity (high HC), whereas hubs in the association cortices and subcortical regions supported higher information traffic (high HT). The high HC observed in temporal, motor, occipital, and midline structures may reflect the large size and/or large number of regions/voxels forming the sensory, motor, and default mode sub-networks (Dosenbach et al. 2010; Sepulcre et al. 2010). In contrast, hubs located in subcortical, associative and multimodal frontal–parietal regions had high HT, implicating their importance in mediating information transmission. This finding is consistent with the known anatomical architecture of the human brain. The thalamus and the basal ganglia are key structures that relay sensory inputs, cortical–cortical communications, and motor outputs (Alexander et al. 1986; Guillery and Sherman 2002). The associative frontal and parietal regions receive and integrate information from distributed cortical regions (Goldman-Rakic 1988). The cerebellum has been shown to have polysynaptic loops with other associative regions including the prefrontal cortex (Middleton and Strick 2001; Kelly and Strick 2003). The anterior temporal pole also had high HT, which could reflect its role as a semantic hub that integrates representations distributed across sensory and motor cortices (Patterson et al. 2007). Hubs in the superior temporal cortex, the temporal pole, and midline structures showed both high HC and HT, consistent with a prior structural connectivity study that classified these regions as the “structural cores” of the human brain (Hagmann et al. 2008).

Note that HC and HT scores were continuous measures, and we acknowledge that the 95th percentile threshold we used to identify hubs is potentially arbitrary. Nevertheless, our simulation results suggested that the functional hubs we identified are critical in facilitating information integration at both global and local scales relative to regions we designated as non-hubs. Removing (lesioning) connections to functional hubs drastically reduced the efficiency of functional brain networks, indicating that long-range transmission of information would occur more slowly and at a higher metabolic cost. The clustering coefficient was also reduced when hubs were lesioned, suggesting that the removal of functional hubs impairs the brain's ability to form tightly connected local clusters or sub-networks, a hallmark of small-world brain networks (Achard et al. 2006). The functional importance of hubs is further underscored by previous studies demonstrating that deficits in hubs could contribute to the symptoms of neurologic and psychiatric disorders (Bassett et al. 2008; Buckner et al. 2009; Cole et al. 2011).

The Development of Functional Brain Networks and Hub Architecture

No age-related differences were found in the location of functional hubs. Neither did we observe differences in the strengths of HC, HT, and node centrality measures in functional hub regions. Our simulated hub lesioning findings suggest that across all age groups, functional hubs are critical in integrating information within tightly connected local clusters and in maintaining the efficiency of information transmission. Furthermore, functional hubs formed a small-world hub-to-hub network in which clustering coefficients, small-worldness, and efficiency of the hub-to-hub network were high regardless of age. It is important to note that non-significant age-related changes in hub scores and graph metrics could be due to the lack of statistical power or high variability across subjects. Although null hypothesis significance tests cannot conclusively prove that these measures do not change with age, we also conducted traditional effect size analyses and Bayesian correlational analyses to estimate the strength of the association between age and hub scores. These analyses corroborated that age-related changes in hub scores are likely small or absent in the age range observed.

Throughout development, brain networks demonstrate small-world organization (Fair et al. 2009; Supekar et al. 2009; Hagmann et al. 2010; Fan et al. 2011), which suggests that certain large-scale network properties may be a foundation that supports developmental increases in the efficiency of information processing. Our results extend this observation by demonstrating the existence of a hub architecture that is stable from ages 10 to 20. We further propose that the existence of a core hub network may be inherent to the functional network architecture of the human brain. This notion is in agreement with data demonstrating that functional hubs are already present in infancy (Fransson et al. 2011), albeit with different spatial locations, and with evidence that centrality measures from structural connectivity do not change with age (Hagmann et al. 2010).

Although the hub architecture was stable from childhood to early adulthood, the strength of hub–spoke connections changed with age. From childhood to adolescence, the strength of connections between frontal hubs and distributed frontal, parietal, temporal, cerebellar spoke regions increased, while connections in the posterior part of the brain decreased. We speculate that this pattern may reflect the gradual maturation of the frontal lobe's ability to coordinate distributed cortical functions for goal-directed behaviors (Miller and Cohen 2001; Luna et al. 2010). For example, we found developmental increases in the strength of connections between the inferior frontal hub and the inferior parietal cortex. These connections are part of the frontal–parietal networks associated with cognitive control functions (Cole and Schneider 2007; Corbetta et al. 2008; Hwang et al. 2010), which develop rapidly from childhood to adolescence (Luna et al. 2004).

Developmental increases in hub–spoke connectivity from adolescence and adulthood were fewer and tended to be more posterior, involving connections between subcortical hubs (the putamen and cerebellum) and frontal, occipital, temporal spoke regions. The development of cerebellar–cortical connections may provide additional resources to assist the fine-tuning of goal-directed behaviors. We identified cerebellar hubs in the lateral cerebellum Crus I and Crus II, which prior RS-fMRI studies have found to be functionally coupled with frontal–parietal association areas (Dosenbach et al. 2007; Krienen and Buckner 2009), and have a homotopic representation of the frontal–parietal cognitive control networks (Buckner et al. 2011). Past research has found that the lateral cerebellum is functionally connected with regions within the frontal–parietal and cingulo-opercular cognitive control networks, and suggests that these cerebellar–cortical connections could be involved in transmitting error-related signals to cognitive control networks for online performance monitoring and optimization (Dosenbach et al. 2007). We speculate that the development of cerebellar–cortical connections could enhance the ability to adapt internal representations supporting goal-directed behavior through refinements in error processing, timing, and learning (Strick et al. 2009), processes which have a protracted development through adolescence into adulthood (Velanova et al. 2008; Luna et al. 2010).

Some hub–spoke connections weakened over time, which may reflect in part the effects of gray matter thinning and synaptic pruning that occur from childhood to early adulthood (Huttenlocher 1979; Rakic et al. 1986; Gogtay et al. 2004). Appropriate pruning of spurious connections may increase the efficacy in information processing, whereas excessive and/or insufficient pruning could contribute to the emergence of psychiatric disorders during adolescence and early adulthood (Rolls and Deco 2011). However, future studies are needed to pinpoint the functional role of these connections.

Although developmental changes in functional connectivity may partially reflect the development of anatomical connectivity, it is important to note that functional connectivity and structural connectivity do not have a direct one-to-one mapping. The development of functional connections could reflect the availability of efficient structural connections (Stevens et al. 2009), or may be shaped by experience (Lewis et al. 2009), such that strong functional connections reflect a history of co-activation even without direct monosynaptic links (Van Dijk et al. 2010). Similarly, the generalizability of graph metrics such as small-worldness and efficiency may be limited when trying to relate network properties observed in functional networks to physical structure, where increased functional connectivity could either reflect a direct anatomical connection or an indirect polysynaptic pathway. Future studies are needed to further explore the bidirectional influence between anatomical structure and functional brain dynamics throughout development.

Conclusion

In conclusion, our results provide new insights into the development of functional brain networks. We found that functional hubs form an efficient network that provides a communication backbone for information transmission in the brain that is present by childhood and remains stable into adulthood. Our results suggest that functional hubs may be an inherent part of human brain architecture, such that the stability of small-world organization and functional hubs likely provides an efficient information processing architecture for maturing brain networks. Hub–spoke connectivity changed considerably over development, suggesting that the fine-tuning of connections between functional hubs and non-hub regions may help form distributed brain networks that support developmental improvements in cognition.

Supplementary Material

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

Funding

This research was supported by the National Institutes of Mental Health (NIMH R01 MH067924 and R01 MH080243 to B.L., and F32 MH090629 to M.N.H.). Computational resources were provided in part by the Center for Simulation and Modeling at the University of Pittsburgh, and by the National Institutes of Health through resources provided by the National Resource for Biomedical Supercomputing (P41 RR06009), which is part of the Pittsburgh Supercomputing Center.

Notes

Conflict of Interest: None declared.

References

Achard
S
Salvador
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
)
Alexander
GE
DeLong
MR
Strick
PL
Parallel organization of functionally segregated circuits linking basal ganglia and cortex
Ann Rev Neurosci
 , 
1986
, vol. 
9
 (pg. 
357
-
381
)
Asato
MR
Terwilliger
R
Woo
J
Luna
B
White matter development in adolescence: a DTI study
Cerebr Cortex
 , 
2010
, vol. 
20
 (pg. 
2122
-
2131
)
Barabasi
AL
Albert
R
Emergence of scaling in random networks
Science
 , 
1999
, vol. 
286
 (pg. 
509
-
512
)
Bassett
DS
Bullmore
E
Verchinski
BA
Mattay
VS
Weinberger
DR
Meyer-Lindenberg
A
Hierarchical organization of human cortical networks in health and schizophrenia
J Neurosci
 , 
2008
, vol. 
28
 (pg. 
9239
-
9248
)
Bassett
DS
Gazzaniga
MS
Understanding complexity in the human brain
Trends Cogn Sci
 , 
2011
, vol. 
15
 (pg. 
200
-
209
)
Bassett
DS
Meyer-Lindenberg
A
Achard
S
Duke
T
Bullmore
E
Adaptive reconfiguration of fractal small-world human brain functional networks
Proc Natl Acad Sci USA
 , 
2006
, vol. 
103
 (pg. 
19518
-
19523
)
Birn
RM
Smith
MA
Jones
TB
Bandettini
PA
The respiration response function: the temporal dynamics of fMRI signal fluctuations related to changes in respiration
Neuroimage
 , 
2008
, vol. 
40
 (pg. 
644
-
654
)
Buckner
RL
Krienen
FM
Castellanos
A
Diaz
JC
Yeo
BT
The organization of the human cerebellum estimated by intrinsic functional connectivity
J Neurophysiol
 , 
2011
, vol. 
106
 (pg. 
2322
-
2345
)
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
)
Cohen
J
Statistical power analysis for the behavioral sciences
 , 
1988
2nd ed. Hillsdale (NJ): Erlbaum
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
)
Cole
MW
Schneider
W
The cognitive control network: integrated cortical regions with dissociable functions
Neuroimage
 , 
2007
, vol. 
37
 (pg. 
343
-
360
)
Corbetta
M
Patel
G
Shulman
GL
The reorienting system of the human brain: from environment to theory of mind
Neuron
 , 
2008
, vol. 
58
 (pg. 
306
-
324
)
Cox
RW
AFNI: software for analysis and visualization of functional magnetic resonance neuroimages
Comput Biomed Res
 , 
1996
, vol. 
29
 (pg. 
162
-
173
)
Dosenbach
NU
Fair
DA
Miezin
FM
Cohen
AL
Wenger
KK
Dosenbach
RA
Fox
MD
Snyder
AZ
Vincent
JL
Raichle
ME
, et al.  . 
Distinct brain networks for adaptive and stable task control in humans
Proc Natl Acad Sci USA
 , 
2007
, vol. 
104
 (pg. 
11073
-
11078
)
Dosenbach
NU
Nardos
B
Cohen
AL
Fair
DA
Power
JD
Church
JA
Nelson
SM
Wig
GS
Vogel
AC
Lessov-Schlaggar
CN
, et al.  . 
Prediction of individual brain maturity using fMRI
Science
 , 
2010
, vol. 
329
 (pg. 
1358
-
1361
)
Eguiluz
VM
Chialvo
DR
Cecchi
GA
Baliki
M
Apkarian
AV
Scale-free brain functional networks
Phys Rev Lett
 , 
2005
, vol. 
94
 pg. 
018102
 
Elo
LL
Jarvenpaa
H
Oresic
M
Lahesmaa
R
Aittokallio
T
Systematic construction of gene coexpression networks with applications to human T helper cell differentiation process
Bioinformatics
 , 
2007
, vol. 
23
 (pg. 
2096
-
2103
)
Fair
DA
Cohen
AL
Power
JD
Dosenbach
NU
Church
JA
Miezin
FM
Schlaggar
BL
Petersen
SE
Functional brain networks develop from a "local to distributed" organization
PLoS Comput Biol.
 , 
2009
, vol. 
5
 pg. 
e1000381
 
Fan
Y
Shi
F
Smith
JK
Lin
W
Gilmore
JH
Shen
D
Brain anatomical networks in early human brain development
Neuroimage
 , 
2011
, vol. 
54
 (pg. 
1862
-
1871
)
Fischl
B
Salat
D
Busa
E
Albert
M
Dieterich
M
Haselgrove
C
van der Kouwe
A
Killiany
R
Kennedy
D
Klaveness
S
Whole brain segmentation automated labeling of neuroanatomical structures in the human brain
Neuron
 , 
2002
, vol. 
33
 (pg. 
341
-
355
)
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
Aden
U
Blennow
M
Lagercrantz
H
The functional architecture of the infant brain as revealed by resting-state fMRI
Cerebr Cortex
 , 
2011
, vol. 
21
 (pg. 
145
-
154
)
Geschwind
DH
Levitt
P
Autism spectrum disorders: developmental disconnection syndromes
Curr Opin Neurobiol
 , 
2007
, vol. 
17
 (pg. 
103
-
111
)
Glass
GV
Peckham
PD
Sanders
JR
Consequences of failure to meet assumptions underlying the fixed effects analyses of variance and covariance
Rev Educ Res
 , 
1972
, vol. 
42
 (pg. 
237
-
288
)
Glover
GH
Li
TQ
Ress
D
Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR
Magn Reson Med
 , 
2000
, vol. 
44
 (pg. 
162
-
167
)
Gogtay
N
Giedd
JN
Lusk
L
Hayashi
KM
Greenstein
D
Vaituzis
AC
Nugent
TF
3rd
Herman
DH
Clasen
LS
Toga
AW
, et al.  . 
Dynamic mapping of human cortical development during childhood through early adulthood
Proc Natl Acad Sci USA
 , 
2004
, vol. 
101
 (pg. 
8174
-
8179
)
Goldman-Rakic
PS
Topography of cognition: parallel distributed networks in primate association cortex
Ann Rev Neurosci
 , 
1988
, vol. 
11
 (pg. 
137
-
156
)
Guillery
RW
Sherman
SM
Thalamic relay functions and their role in corticocortical communication: generalizations from the visual system
Neuron
 , 
2002
, vol. 
33
 (pg. 
163
-
175
)
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
 
Hagmann
P
Sporns
O
Madan
N
Cammoun
L
Pienaar
R
Wedeen
VJ
Meuli
R
Thiran
JP
Grant
PE
White matter maturation reshapes structural connectivity in the late developing human brain
Proc Natl Acad Sci USA
 , 
2010
, vol. 
107
 (pg. 
19067
-
19072
)
He
Y
Chen
Z
Evans
A
Structural insights into aberrant topological patterns of large-scale cortical networks in Alzheimer's disease
J Neurosci
 , 
2008
, vol. 
28
 (pg. 
4756
-
4766
)
Holm
S
A simple sequentially rejective multiple test procedure
Scand J Stat
 , 
1979
, vol. 
6
 (pg. 
65
-
70
)
Huttenlocher
PR
Synaptic density in human frontal cortex – developmental changes and effects of aging
Brain Res
 , 
1979
, vol. 
163
 (pg. 
195
-
205
)
Hwang
K
Velanova
K
Luna
B
Strengthening of top-down frontal cognitive control networks underlying the development of inhibitory control: a functional magnetic resonance imaging effective connectivity study
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
15535
-
15545
)
Jo
HJ
Saad
ZS
Simmons
WK
Milbury
LA
Cox
RW
Mapping sources of correlation in resting state fMRI, with artifact detection and removal
Neuroimage
 , 
2010
, vol. 
52
 (pg. 
571
-
582
)
Kelly
AMC
Di Martino
A
Uddin
LQ
Shehzad
Z
Gee
DG
Reiss
PT
Margulies
DS
Castellanos
FX
Milham
MP
Development of anterior cingulate functional connectivity from late childhood to early adulthood
Cerebr Cortex
 , 
2009
, vol. 
19
 (pg. 
640
-
657
)
Kelly
RM
Strick
PL
Cerebellar loops with motor cortex and prefrontal cortex of a nonhuman primate
J Neurosci
 , 
2003
, vol. 
23
 (pg. 
8432
-
8444
)
Krienen
FM
Buckner
RL
Segregated fronto-cerebellar circuits revealed by intrinsic functional connectivity
Cereb Cortex
 , 
2009
, vol. 
19
 (pg. 
2485
-
2497
)
Legendre
P
Species associations: the Kendall coefficient of concordance revisited
J Agric Biol Environ Statistics
 , 
2005
, vol. 
10
 (pg. 
226
-
245
)
Lewis
CM
Baldassarre
A
Committeri
G
Romani
GL
Corbetta
M
Learning sculpts the spontaneous activity of the resting human brain
Proc Natl Acad Sci USA
 , 
2009
, vol. 
106
 (pg. 
17558
-
17563
)
Luna
B
Garver
KE
Urban
TA
Lazar
NA
Sweeney
JA
Maturation of cognitive processes from late childhood to adulthood
Child Develop
 , 
2004
, vol. 
75
 (pg. 
1357
-
1372
)
Luna
B
Padmanabhan
A
O'hearn
K
What has fMRI told us about the development of cognitive control through adolescence?
Brain Cogn
 , 
2010
, vol. 
72
 (pg. 
101
-
113
)
Middleton
FA
Strick
PL
Cerebellar projections to the prefrontal cortex of the primate
J Neurosci
 , 
2001
, vol. 
21
 (pg. 
700
-
712
)
Miller
EK
Cohen
JD
An integrative theory of prefrontal cortex function
Annu Rev Neurosci
 , 
2001
, vol. 
24
 (pg. 
167
-
202
)
Morey
RD
Rouder
JN
Bayes factor approaches for testing interval null hypotheses
Psyc Methods
 , 
2011
, vol. 
16
 (pg. 
406
-
419
)
Patterson
K
Nestor
PJ
Rogers
TT
Where do you know what you know? The representation of semantic knowledge in the human brain
Nat Rev Neurosci
 , 
2007
, vol. 
8
 (pg. 
976
-
987
)
Perkins
AD
Langston
MA
Threshold selection in gene co-expression networks using spectral graph theory techniques
BMC Bioinformatics
 , 
2009
, vol. 
10
 
Suppl 11
pg. 
S4
 
Power
JD
Barnes
KA
Snyder
AZ
Schlaggar
BL
Petersen
SE
Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion
Neuroimage
 , 
2012
, vol. 
59
 (pg. 
2142
-
2154
)
Power
JD
Fair
DA
Schlaggar
BL
Petersen
SE
The development of human functional brain networks
Neuron
 , 
2010
, vol. 
67
 (pg. 
735
-
748
)
Rakic
P
Bourgeois
JP
Eckenhoff
MF
Zecevic
N
Goldman-Rakic
PS
Concurrent overproduction of synapses in diverse regions of the primate cerebral cortex
Science
 , 
1986
, vol. 
232
 (pg. 
232
-
235
)
Rolls
ET
Deco
G
A computational neuroscience approach to schizophrenia and its onset
Neurosci Biobehav Rev
 , 
2011
, vol. 
35
 (pg. 
1644
-
1653
)
Rubinov
M
Sporns
O
Complex network measures of brain connectivity: uses and interpretations
Neuroimage
 , 
2010
, vol. 
52
 (pg. 
1059
-
1069
)
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
)
Sepulcre
J
Liu
H
Talukdar
T
Martincorena
I
Yeo
BTT
Buckner
RL
The organization of local and distant functional connectivity in the human brain
PLoS Comput Biol
 , 
2010
, vol. 
6
 pg. 
e1000808
 
Smith
SM
Jenkinson
M
Woolrich
MW
Beckmann
CF
Behrens
TE
Johansen-Berg
H
Bannister
PR
De Luca
M
Drobnjak
I
Flitney
DE
, et al.  . 
Advances in functional and structural MR image analysis and implementation as FSL
Neuroimage
 , 
2004
, vol. 
23
 
Suppl 1
(pg. 
S208
-
S219
)
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
Chialvo
DR
Kaiser
M
Hilgetag
CC
Organization, development and function of complex brain networks
Trends Cogn Sci (Regul Ed)
 , 
2004
, vol. 
8
 (pg. 
418
-
425
)
Sporns
O
Honey
CJ
Kötter
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
 
Stevens
MC
Skudlarski
P
Pearlson
GD
Calhoun
VD
Age-related cognitive gains are mediated by the effects of white matter development on brain network integration
Neuroimage
 , 
2009
, vol. 
48
 (pg. 
738
-
746
)
Strick
PL
Dum
RP
Fiez
JA
Cerebellum and nonmotor function
Annu Rev Neurosci
 , 
2009
, vol. 
32
 (pg. 
413
-
434
)
Supekar
K
Musen
M
Menon
V
Development of large-scale functional brain networks in children
PLoS Biol
 , 
2009
, vol. 
7
 pg. 
e1000157
 
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
 , 
2012
, vol. 
59
 (pg. 
431
-
438
)
Van Dijk
KRA
Hedden
T
Venkataraman
A
Evans
KC
Lazar
SW
Buckner
RL
Intrinsic functional connectivity as a tool for human connectomics: theory, properties, and optimization
J Neurophysiol
 , 
2010
, vol. 
103
 (pg. 
297
-
321
)
Velanova
K
Wheeler
ME
Luna
B
Maturational changes in anterior cingulate and frontoparietal recruitment support the development of error processing and inhibitory control
Cerebr Cortex
 , 
2008
, vol. 
18
 (pg. 
2505
-
2522
)
Watts
DJ
Strogatz
SH
Collective dynamics of ‘small-world’ networks
Nature
 , 
1998
, vol. 
393
 (pg. 
440
-
442
)
Zuo
X-N
Ehmke
R
Mennes
M
Imperati
D
Castellanos
FX
Sporns
O
, et al.  . 
Network centrality in the human functional connectome
Cerebr Cortex
 , 
2012
, vol. 
22
 (pg. 
1862
-
1875
)