Abstract
Age-related changes in human functional neuroanatomy are poorly understood. This is partly due to the limits of interpretation of standard fMRI. These limits relate to age-related variation in noise levels in data from different subjects, and the common use of standard adult brain parcellations for developmental studies. Here we used an emerging MRI approach called multiecho (ME)-fMRI to characterize functional brain changes with age. ME-fMRI acquires blood oxygenation level-dependent (BOLD) signals while also quantifying susceptibility-weighted transverse relaxation time (T2*) signal decay. This approach newly enables reliable detection of BOLD signal components at the subject level as opposed to solely at the group-average level. In turn, it supports more robust characterization of the variability in functional brain organization across individuals. We hypothesized that BOLD components in the resting state are not stable with age, and would decrease in number from adolescence to adulthood. This runs counter to the current assumptions in neurodevelopmental analyses of brain connectivity that the number of BOLD signal components is a random effect. From resting-state ME-fMRI of 51 healthy subjects of both sexes, between 8.3 and 46.2 years of age, we found a highly significant (r = −0.55, p ≪ 0.001) exponential decrease in the number of BOLD components with age. The number of BOLD components were halved from adolescence to the fifth decade of life, stabilizing in middle adulthood. The regions driving this change were dorsolateral prefrontal cortices, parietal cortex, and cerebellum. The functional network of these regions centered on the cerebellum. We conclude that an age-related decrease in BOLD component number concurs with the hypothesis of neurodevelopmental integration of functional brain activity. We show evidence that the cerebellum may play a key role in this process.
SIGNIFICANCE STATEMENT Human brain development is ongoing from childhood to at least 30 years of age. Functional MRI (fMRI) is key for characterizing changes in brain function that accompany development. However, developmental fMRI studies have relied on reference maps of adult brain organization in the analysis of data from younger subjects. This approach may limit the characterization of functional activity patterns that are particular to children and adolescents. Here we used an emerging fMRI approach called multi-echo fMRI that is not susceptible to such biases when analyzing the variation in functional brain organization over development. We hypothesized an integration of the components of brain activity over development, and found that the number of components decreases exponentially, halving from 8 to 35 years of age. The brain regions most affected underlie executive function and coordination. In summary, we show major changes in the organization and integration of functional networks over development into adulthood, with both methodological and neurobiological implications for future lifespan and disease studies on brain connectivity.
Introduction
Characterizing brain development from adolescence to adulthood is critical for understanding neuropsychiatric disease and healthy brain function. However, the trajectories of changes in functional organization during brain development are not yet well characterized. Developmental studies of white matter structural change based on diffusion weighted MRI report nonlinear trajectories involving faster changes at earlier ages, followed by stabilization at later ages (Wallace et al., 2006). Microstructural changes in white matter are also known to be tract and region specific, suggesting that trajectories are both global and regionally specific (Paus et al., 2001). In this study, we use advanced functional MRI (fMRI) techniques to address age-related trajectories of change in functional brain organization at both the whole-brain and regional levels.
Recent findings indicate that neurodevelopmental changes in the organization of functional networks are detectable as age-related increases in network coherence (Gu et al., 2015). It follows that networks of functional correlation may become more integrated from adolescence into adulthood, changing the broader organization of functional networks. However, age-related differences in anatomy complicate the comparison of functional activity data across development (Power et al., 2012). To more sensitively detect trajectories of age-related functional brain change, we used multiecho (ME) fMRI (Kundu et al., 2015) to image the resting state. ME-fMRI isolates the susceptibility-weighted transverse relaxation (T2*) component of blood oxygenation level-dependent (BOLD) fMRI signals without many of the arbitrary denoising models used for standard fMRI. Furthermore, using ME-fMRI, BOLD signal components can be characterized without standard brain parcellations derived from adult data, which are commonly used, but may bias results to represent normative adult anatomy (Power et al., 2011; Craddock et al., 2012).
In this study, we especially focused on characterizing the number of functional BOLD components in the resting state as a marker of brain development. The number of BOLD components in resting-state fMRI data may be considered to represent the dimensionality, degrees of freedom, or dimensionality of spontaneous brain activity (Friston et al., 1995). We considered that this parameter could be useful in representing the level of fragmentation versus the integration of functional networks across the brain. The variation of BOLD component number with age is evaluated in this article. Specifically, we had the following three aims: (1) to show that age affects the number of components in the fMRI signal, as detected experimentally using ME-independent component analysis (ICA); (2) to define data-driven brain regions that are particularly susceptible to the impact of age in terms of the number of BOLD components they express; and (3) to characterize the relationships among brain regions that show similar trajectories of change in regional BOLD component number. Altogether, we find that functional brain networks undergo a trajectory of functional integration with age, in a regionally specific way, which has implications for the development of normal cognition and behavior as well as neurodevelopmental disorders.
Materials and Methods
Overview
We assessed three levels of functional brain organization. First, we determined the number of BOLD signal components at the individual subject level. Second, for each subject, we computed a map of the number of components corepresented in each voxel. Third, we computed the relationship between component corepresentation and subject age, yielding regions of interest (ROIs) of age-related component change. We then used these ROIs as seed regions to estimate their functional connectivity and elucidate age-sensitive functional networks.
Experimental design
Participants.
Fifty-one healthy subjects (mean age, 21.9 years; age range, 8.3–46.2 years; 20 females) completed the study. This study was approved by the National Institutes of Health Institutional Review Board.
Image acquisition.
Data were acquired on a GE MR750 3 T Scanner using a 32-channel receive-only head coil (GE Healthcare). Each imaging session first involved acquiring a whole-brain anatomical MPRAGE scan with 1 mm isotropic resolution. The resting-state fMRI scan was 10 min long and involved acquisition of multiecho time-course EPI using the following parameters: 240 cm field of view; 64 × 64 resolution yielding 3.75 isotropic voxels; in-plane SENSE (sensitivity encoding) acceleration factor, 2; flip angle, 77°; repetition time (TR), 2.0 s; and echo times (TEs), 12.8, 28, and 43 ms. The ME-fMRI sequence was implemented using vendor EPI excitation and a modified EPI readout, and used on-line reconstruction (Poser et al., 2006). Each TR corresponded with the acquisition of three volumes having TE values of 15, 35, and 58 ms.
Anatomical and functional imaging processing.
Anatomical images were first processed using the FreeSurfer pipeline for skull stripping, segmentation, and cortical surface mesh construction. Separately, functional images were processed using the ME-ICA pipeline as implemented in the AFNI meica.py toolbox. This toolbox implements preprocessing, decomposition, and denoising steps for multiecho EPI data, which are detailed in the study by Kundu et al. (2015). Briefly, multiecho EPI time series datasets of each TE were aligned for slice-timing offsets, all volumes were aligned with rigid-body motion correction to the volume of the first TR, and functional images were skull stripped. These processed functional images were used to compute T2* and S0 maps by fitting signal means of different TEs for each voxel to a monoexponential decay model. Using meica.py, parameters of affine coregistration between (skull-stripped) anatomical and functional images were estimated, using the thresholded T2* map as a weight volume (Kundu et al., 2015). Anatomical images were nonlinearly warped to MNI standard space using the AFNI 3dQwarp technique (Cox, 2012), then a single warp combining deobliquing, motion correction, anatomical–functional coregistration, and nonlinear warp to MNI space was applied to each echo dataset separately. This procedure produced the preprocessed datasets for input into decomposition and denoising analyses.
Statistical analysis
BOLD components.
After preprocessing, the next stage of ME-ICA was to extract BOLD components. The processing steps of the ME-ICA pipeline are summarized here and have been detailed in prior publications (Kundu et al., 2015, 2017). The first step of ME-ICA was to compute a weighted average of the separate echo time series into a single optimized time series dataset. The weighting function was evaluated at each voxel, factoring in TE and voxel-specific estimates of T2* (Posse et al., 1999). The second step was to estimate the total number of components and to remove Gaussian distributed components from the data. This step is standard for all applications of ICA to fMRI data, and is usually performed using a variant of principal components analysis (PCA). The ME-ICA pipeline performs this step using ME-PCA (Kundu et al., 2013). The third step was to conduct the ICA on the resulting data and elucidate functional BOLD and artifact components (Kundu et al., 2012). Finally, the ICA components were organized into BOLD and non-BOLD categories using metrics of TE dependence and TE independence. Strong BOLD weighting of a component was expressed in high values of the pseudo-F statistic κ, and strong non-BOLD (artifact) weighting was expressed in high values of the pseudo-F statistic ρ. Functional BOLD components were interpreted as those with high κ values and low ρ values (Kundu et al., 2012, 2017). All remaining components are considered to be non-BOLD components. This component classification is also used to denoise fMRI time series. The non-BOLD components are projected out of the time series data, based on multiple least-squares fit of the entire mixing matrix. Maps of the signal-to-noise ratio (SNR) of these time series data were computed for each subject by dividing the voxel means by the standard deviation (SD).
Global age-related variation of BOLD components.
The total number of BOLD components per subject was compared against subject age. The resting-state fMRI dataset from each subject produced one such measure. Based on prior reports of structural brain changes with age (Giedd et al., 1999; Paus et al., 2001), the relationship between component number and age was computed using linear, quadratic, exponential, and power law functions.
The validity of values of BOLD dimensionality derived from ME-ICA was further evaluated. First, BOLD dimensionality values were correlated against values of framewise displacement derived from rigid-body head motion parameters. This was done to determine whether BOLD dimensionality was significantly influenced by head motion. Second, BOLD and non-BOLD dimensionality values (determined simultaneously in ME-ICA) were correlated against each other. This was done to assess whether greater non-BOLD (i.e., artifact) signal variance predetermined the number of BOLD components. Last, ME-ICA-denoised BOLD time series were analyzed with probabilistic PCA (PPCA) from the FSL program MELODIC. This give an independent estimate of dimensionality for comparison with ME-ICA estimates of BOLD dimensionality.
We sought further to establish BOLD dimensionality as a metric of integration of functional brain activity, via comparison to corresponding graph theoretic measures based on functional connectomes. First, for each subject, functional connectomes were constructed. Each subject's Deskian-Killany FreeSurfer brain parcellation was coregistered to their preprocessed ME-ICA-derived BOLD functional time series. Parcel average time series were then computed, and used to compute a matrix of Pearson correlations. These matrices were then thresholded. Several thresholds were sampled: r = 0.3–0.9, in Δr = 0.1 increments. Negative correlations were removed, so only positive correlations were considered. Thresholded time series correlation matrices were then Fisher R–Z transformed to create functional connectomes. Graph theoretic analysis of connectomes was performed using the Brain Connectivity Toolbox, implemented in MATLAB (Rubinov and Sporns, 2010). For each subject connectome (weighted, undirected), Louvain community detection was applied. Per-node participation coefficient values were then computed based on those communities. The mean participation coefficient was used as a per-subject measure of functional connectomic integration. These per-subject values of mean participation coefficient were then compared against BOLD dimensionality using Pearson (parametric) and Spearman (nonparametric) correlations. This procedure was repeated for each sampled threshold on time series correlation to check for robustness of the comparison.
Regional age-related variation of BOLD components.
The BOLD components identified using ME-ICA were used to delineate functional regions. A region was defined as a set of contiguous voxels with a common set of overlapping BOLD components. The effect of a component on a region was determined by the significance of BOLD TE dependence across the contiguous voxels in the component map. Component maps were rendered using the F-R2* statistic (Kundu et al., 2012). This statistic indicates the level of TE dependence [i.e., the susceptibility-weighted transverse relaxation rate (R2*), the inverse of T2*]. F-R2* is computed voxelwise and follows a standard F(1,2) distribution. BOLD components (i.e., functional networks) were each rendered in F-R2* units, thresholded (p < 0.01, uncorrected, as per prior work on statistical significance for TE dependence maps, Kundu et al., 2012), and binarized. Binary maps of all BOLD components were summed, producing a map of the number of components with significant weight at each voxel, in effect serving as a map of component overlap. Next, the relationship of how the number of overlapping BOLD components changed with age was mapped. A functional brain map showing the number of overlapping BOLD components per voxel was rendered for each subject, called a T2* component overlap map (T2*-COLAP). This map was normalized to standard cortical surface space using FreeSurfer cortical surface meshes and AFNI SUMA tools (Argall et al., 2006; Fischl, 2012), with subcortical regions from nonlinear warps merged to create a whole-brain standard space map.
Across subjects, for each voxel in the T2*-COLAP maps in standard space, we computed a nonparametric Spearman correlation of the number of overlapping BOLD components versus subject age. This produced a new map reflecting Spearman correlation values, which was thresholded to Spearman's ρ for p < 0.01 significance, as shown in Figure 4B. The resulting parametric map was then cluster corrected, based on Monte Carlo α probability simulations, to α < 0.05. First, the smoothness of preprocessed multiecho functional data was estimated in terms of spatial autocorrelation using AFNI 3dFWHMx (with the −ACF option). Then, the cluster probability table for nearest-neighbor clustering (considering both positive and negative values) was calculated using 3dClustSim (Cox et al., 2017). Finally, spatial clustering was applied using AFNI 3dmerge. For the target cluster probabilities, a cluster extent of 38 voxels was required.
A separate spatial-clustering strategy for regional BOLD dimensionality correlation maps, called density-based spatial clustering of applications with noise (DBSCAN; Ester et al., 1996), was also evaluated. This technique, a standard multivariate clustering technique, can deterministically cluster arbitrary feature spaces. Here we used DBSCAN on a feature space including both spatial coordinates (i.e., spatial clustering) as well as a data variate of the Spearman correlation values, as described above. DBSCAN formally distinguishes clusters and noise, defining clusters on density and noise in terms of being outside of “density reachability.” Practically, this application of DBSCAN is well suited to spatially clustering a densely populated statistical parametric map while rejecting voxels with neither high value nor contiguity with larger clusters. In such cases, otherwise distinct clusters tend to have some touching voxels, yielding apparent contiguity and thus very large clusters. Instead of applying an arbitrary refinement scheme like erosion and dilation after spatial clustering, we chose DBSCAN as an optimization-based clustering technique to solve the spatial contiguity problem. The Python scikit-learn implementation of DBSCAN was used. A parameter search for the single DBSCAN parameter η was conducted, based on a minimum cluster size of 40 to find the solution that yielded the maximum number of clusters. Each surviving cluster was then treated as a separate mask. For each subject map, component overlap values were averaged within each mask, and values were correlated against age for each mask region. The relationships between average regional BOLD component overlap and subject age was assessed for linear, exponential, and power law relationships.
Cross-validation of age as predicted by regional BOLD component overlap.
The extent to which patterns of regional BOLD component overlap were predictive of age was assessed using support vector machine (Scholkopf et al., 1997), cross-validation, and permutation testing, implemented in Python software (Pedregosa et al., 2011). The BOLD component overlap values from all voxels of T2*-COLAP maps across subjects were used as features for age prediction. Prediction was made using the support vector regression (SVR) framework, training on values of age. We confirmed that raw BOLD component counts gave inferior prediction compared with log-linearized counts, so the latter is shown in the results. The accuracy of the age prediction from BOLD component overlap was determined using leave-one-out (LOO) cross-validation. The mean and standard deviation of the prediction error were computed. The significance of classification was determined by permutation testing. Using 1000 randomly permutated assignments of ages to datasets, the SVR was trained and tested in LOO cross-validation, with median absolute error as the loss function. In effect, the significance of association between “true” age and BOLD dimensionality was established by comparison to chance association between an “age-like” variate and dimensionality values, and was expressed in a significance (p) value.
Seed-based functional connectivity of regions with age-related change in BOLD components.
Functional connectivity analysis was conducted based on time series extracted from those brain regions that showed a significant change in BOLD component overlap with age. For each subject, voxel time series within each respective cluster mask were averaged. The Pearson correlation between these time series was computed, and was scaled to correct for the effective degrees of freedom for correlation (i.e., the number of BOLD independent components comprising the denoised time series) according to the following (Kundu et al., 2013):
For each region, the Fisher Z map of functional connectivity for each subject was input to a one-sample t test group analysis, producing a group map of functional connectivity for that region. Each such group map was thresholded with cluster correction for multiple comparisons based on α probability simulation (as described above) to α < 0.01 (Cox et al., 2017). Then, the overlap of each thresholded group-level regional connectivity map was computed as a count of voxels. These corrected maps were used to determine which regions had seed-connectivity maps that overlapped with seed-connectivity maps of the other regions, which was represented in a graph.
Covariation of functional connectivity with age.
We tested the relationship between regional changes in component overlap versus functional connectivity and its change with age. The clusters of the map of Spearman correlation between T2*-COLAP maps and subject age were used to identify centroids to use as seed voxels for functional connectivity analysis. The Pearson correlations of each such time course against the time courses of all other voxels within the brain mask were computed, yielding whole-brain connectivity maps. Time series correlation values were normalized using the Fisher R–Z transform. The standard error (SE) term that accounted for the variability across datasets in terms of the total number of BOLD components (Eq. 1) was included, because false-positive error may be introduced when this factor is not accounted for (Kundu et al., 2013). Group analyses of these connectivity maps were then performed using mass-univariate Pearson correlation analysis of connectivity versus age, voxelwise. For each seed, the number of voxels with significant (p < 0.01) group-level correlation of connectivity with age was counted, separately for positive and negative values. A two-sample paired t test was then used to compare the number of positively versus negatively age-correlated voxels, to determine whether there was a net increase or decrease of functional connectivity with age across regions that showed decreasing regional BOLD dimensionality.
Results
ME-ICA of resting-state fMRI across adolescent and adult data
The separation of ME-ICA components into BOLD and non-BOLD categories was conducted successfully for each subject. κ and ρ metrics for datasets of a representative subsample of subjects between 8.25 and 46.17 years are shown in Figure 1, A and B, respectively. κ metrics are plotted in Scree plots in descending order, juxtaposed with corresponding ρ values (i.e., in the order of κ values). Together, a change is apparent from high-κ/low-ρ components to low-κ/high-ρ components. These two categories represent the two distinct BOLD and non-BOLD sets of components, from which respective total component numbers were determined. Figure 1 shows for, an adolescent and an adult subject, maps of the “projections” of all components into a single map, in terms of the most statistically significant clusters (Fig. 1C,D). On visual inspection of these maps, the adolescent subject dataset shows a wider coverage of gray matter by nodes of detected functional networks, including in subcortex and cerebellum.
A, B, κ and ρ values (i.e., pseudo-F statistics) for component-level TE dependence and TE independence of ME-ICAs for youngest (A) and oldest (B) subjects expressed as Scree plots, sorted by κ values. Functional network components are characterized by high κ and low ρ values. Note the transition from the κ ≫ ρ (BOLD) regime to the κ ≤ ρ (non-BOLD) regime. C, D, Overlap of BOLD components (thresholded α < 0.01) from subjects of the youngest and oldest ages, respectively (8.25 and 46.17 years of age). The colored regions are derived from thresholded BOLD components. These component maps are overlaid in order of decreasing component size (in terms of number of significant voxels). C, Overlap of 53 BOLD components. D, Overlap of 13 BOLD components.
Global BOLD component number decreases with age
Global BOLD component numbers decreased with age in all tested relations, as follows: power law, exponential, and linear (R2 = 0.32, 0.33, 0.28, respectively). Agreement with all three relations is consistent in the present case as they were all monotonic, and the effect size was large. The maximum global BOLD component number was seen for a 9-year-old subject, who expressed 95 BOLD networks. The minimum value was for a 42-year-old subject, who expressed 32 BOLD networks. Figure 2 shows associations of subject age with global BOLD component number. In addition, we confirmed that global BOLD component number was not correlated to the level of subject head motion. The number of BOLD components was also not correlated to the number of non-BOLD components (r = 0.04, p = 0.7). In contrast, the level of subject head motion was correlated to the number of non-BOLD components according to a positive linear relationship, which is shown in Figure 3.
A, Scatterplot of global BOLD component numbers across subjects. Power law fit (blue) is shown for the scaling of component number versus subject age. Significant linear (L) and power law (PL) trends were found (L: R2 = 0.29, p < 0.001; PL: R2 = 0.30, p < 0.001). B, Scatterplot of global BOLD component number versus corresponding subject values of mean participation coefficient, based on functional connectomes of the Deskian–Killany parcellation and Louvain community detection.
Scatterplots of non-BOLD component number across subjects showing significant linear (L) or power law (PL) fits. A, Total BOLD component number versus total non-BOLD component number, no significant relation (r = 0.05, p = 0.7). B, Global BOLD component number versus subject motion, no significant fits (L: R2 = 0.02, p = 0.26; PL: R2 = 0.02, p = 0.24). C, Subject motion versus age, no significant relation (L: R2 = 0.02, p = 0.31; PL: R2 = 0.01, p = 0.58). D, Global non-BOLD dimensionality versus subject motion, significant linear relationship (L: R2 = 0.30, p = 0.001; PL: R2 = 0.30, p = 0.001).
We compared the estimate of BOLD dimensionality from ME-ICA to a separate estimate of dimensionality provided by PPCA of the time series with non-BOLD signals removed. We found a highly significant correlation between the two dimensionality estimates (r = 0.97, p ≪ 10−7). BOLD dimensionality from ME-ICA also correlated with graph theoretic measures of integration derived from subject-level functional connectomes. Across connectivity thresholds corresponding to r = 0.4–0.8, the mean participation coefficient based on thresholded, weighted, undirected functional connectomes was significantly (p < 0.05) correlated with global BOLD dimensionality. The connectivity threshold that led to the most significant such correlation was R > 0.5, leading to correlations between BOLD dimensionality and participation coefficient of Pearson's r = −0.34 (p = 0.013) and Spearman's ρ = −0.38 (p = 0.0057).
Regional BOLD component count and its reduction with age
Maps of BOLD components were generated in units of TE dependence (F-R2*; thresholded at α < 0.01, corrected). The mapping of the overlap of BOLD components, the T2*-COLAP map, is shown for a representative subject in Figure 4A. The TE dependence maps of ICA components on which the T2*-COLAP map is based show spatial patterns comparable with conventional amplitude-based maps from spatial ICA. The components mapped with TE dependence as found in individual subject data included canonical resting-state networks such as default mode, sensorimotor, and frontoparietal networks (Damoiseaux et al., 2006).
Assessment of regional BOLD component overlap versus subject age. A, Procedure: thresholded (p < 0.01, uncorrected) F-maps indicating component voxels with significant TE dependence are warped into standard (MNI) space via normalized cortical surface mesh and nonlinear volumetric warp (see Materials and Methods), then are binarized and summed to produce an overlap map for component-level TE dependence indicating regional BOLD dimensionality. These were termed T2*-COLAP maps. B, Voxels (grouped into clusters post hoc) showing significant Spearman correlation of component overlap with subject age. Clusters represent significance with correction for multiple comparisons (α < 0.01), rendered with a high-contrast color map. C, Scatterplots representing cluster averages for regional BOLD component overlap plotted against subject age. A high magnitude of power law behavior is apparent (note that the initial voxelwise correlation was nonparametric), which varies in parameters depending on brain region. BA, Brodmann's area; Cereb., cerebellum; L., left; R., right; PreSMA, pre-supplementary motor area; Ling., lingual gyrus; PCC, posterior cingulate cortex; MFg, middle frontal gyrus; PoCg., posterior cingulate gyrus; STg, superior temporal gyrus; IFg, inferior frontal gyrus; Pcn, precuneus; PoPut, posterior putamen; SPl, superior parietal lobule; IPl, inferior parietal lobule; MeFg, medial frontal gyrus; MOg, middle occipital gyrus; aIns, anterior insula.
For a representative subject, Figure 4A illustrates the result of thresholding, binarizing, and summing over BOLD component maps to produce a component overlap map. The overlap map conveys the number of BOLD components with a significant linear TE dependence at each voxel. For example, high component overlap is observed for regions associated with the default-mode network, and low component overlap is observed in regions such as the motor cortex.
TE dependence maps and their overlap showed two additional aspects of component maps of functional networks. One is that functional networks, rendered in units of TE dependence, are typically inclusive of widespread brain regions. This is attributed to low-amplitude contributions of components being detected in TE dependence maps that may not have been detected in maps in magnitude units thresholded according to Z-scores. Second is the indication that TE dependence maps of networks can be more comprehensive than magnitude-based maps on the basis of greater overlap across components, which is not seen in the case of magnitude-based maps.
Regional BOLD component reduction with age
The relationship between the number of overlapping BOLD components and subject age at group level was mapped for each voxel using a Spearman correlation. After voxel and cluster thresholds were applied to this correlation map (see Materials and Methods), a distinct set of clusters was found. These clusters indicated an age-related decrease in BOLD component count in gray matter regions. The cluster averages of regional BOLD component overlap versus age showed highly significant decreases of BOLD component count with age across regions (Fig. 4B). Thus, at both voxel and cluster levels, age-related decreases in component count were observed. Importantly, the reduction of component count was not homogeneous across the brain. Figure 4B maps specific regions that showed this reduction, averaged over the voxels of the region, as follows: prefrontal cortex (i.e., bilateral dorsolateral prefrontal cortex, medial prefrontal cortex, frontopolar cortex), bilateral superior parietal cortex (i.e., sensory cortex, precuneus), and bilateral cerebellar hemispheres spanning Crus VI-VIIab. Regions that did not show age-related change in BOLD component overlap included premotor and primary motor cortices, the temporal lobes, insula, thalamus, and inferior frontal cortex.
Linear, exponential, and power law models were used to characterize regional age-related reduction of BOLD component overlap with age. The power law model gave the most statistically significant fit (r = −0.54, p ≪ 0.001). The finding of regional reduction in BOLD component overlap was consistent with the global reduction of BOLD component number, but with even greater statistical significance at the regional level than the global level.
Cross-validation of subject age as predicted by regional variation of BOLD component overlap
Regional values of BOLD component overlap were used to predict subject age using a SVR and cross-validation strategy (Fig. 5). This approach predicted age with an average error of 6.5 ± 0.6 years, with error based on leave-one-out cross-validation. The significance of the classification, based on cross-validation involving 1000 permutations of age values, was p = 0.009. These results suggest that regional BOLD component overlap is a significant predictor of subject age on the order of years.
Predicted versus actual age from training on regional BOLD component overlap in LOO cross-validation with support vector regression. The mean error in estimates of age is 6.5 ± 0.6 years. The significance of classification based on permutation testing (1000 permutations) is p = 0.009.
Seed-based functional connectivity between regional changes in BOLD component overlap
Seed-based connectivity was computed between functional regions that showed age-related changes in the number of overlapping BOLD components. Seed-based connectivity was computed using the complete Fisher transformation including degrees of freedom in terms of the global BOLD component number, a subject-level variable. The group-level connectivity between these regions was represented in a circular network diagram (Fig. 6). Each region of this network showed connectivity to at least one of the other regions with age-related change in component overlap. Bilateral cerebellum and right precuneus were most connected to the other regions (i.e., these were the highest degree nodes).
Circle plot representing the graph of average inter-regional connectivity between regions with decreasing BOLD component overlap versus subject age. Each edge that is shown represents positive connectivity of a seed and target region, such that connectivity is detected in the target region with cluster-level significance of α < 0.01. Self-connections are not shown. The degree of connectivity of each region is shown, with dark coloration representing a relatively low degree of connectivity, and bright red coloration representing a high degree of connectivity. A high degree of connectivity is observed of left and right cerebellar nodes as well as right precuneus. Inf, inferior; Sup, superior; Mid, middle; Gyr, gyrus; Ant, anterior; Med, medial; Lob, lobule; fi (e.g. f2, f3), focus i; Cb, cerebellar lobule; PFop, PFcm, PGa are cytoarchitectonic regions as defined in von Economo and Koskinas, 1926.
Correlation of seed connectivity with age after controlling for total BOLD components
We found that the regions with decreasing regional BOLD dimensionality with age had increasing functional connectivity with age. Across these regions, we found that the number of voxels with significant trends of increasing functional connectivity with age were significantly greater (p = 0.0004) than voxels with decreasing trends of connectivity with age (Fig. 7A). The two-sample t test design of this analysis was robust to false-positive error. The right cerebellum showed particularly high change of functional connectivity with age. The analysis of seed-based functional connectivity of the cerebellum versus subject age is shown in Figure 7, B and C. Age-related connectivity was seen from the cerebellum to a bilateral frontoparietal network, caudate head, and medial thalamus (Fig. 7B). Age-related connectivity increases were greater in the right hemisphere. Cluster averages of the values of connectivity across subjects versus subject age are shown in Figure 7C. Except for the right parietal lobe, which showed decreasing connectivity with age, all plots showed significant age-related increases in connectivity with age, even after controlling for global age-related change in BOLD component number.
A, Bar charts comparing number of voxels with positive versus negative trends (p < 0.05) of change in seed-based functional connectivity with age, with seeds being those regions with changing BOLD component overlap with age. B, Group-level seed-based functional connectivity map of left cerebellum (Lobule VII-Crus II), using central voxel of cerebellar cluster as seed and one-sample t test, thresholded α < 0.01 (cluster corrected). Blue colors indicate mean negative functional connectivity. C, For left cerebellum seed connectivity, scatterplots showing mean connectivity per subject versus subject age, within clusters (α < 0.01) showing age-correlated connectivity change. Linear fits are overlaid, with significance of the Pearson correlation in the legend.
Discussion
In this study, we showed that age modulates the number of components of functional brain activity in the resting state. This association suggests increasing age-related integration of functional brain networks, the first such demonstration in resting-state fMRI data. The number of BOLD components varies with age parametrically, as a power law. This pattern is highly robust. It is expressed with high statistical significance at the global whole-brain level, and with even greater significance locally, in individual brain regions, expressed as a change in the number of overlapping BOLD components as a function of age. The effect is not expressed across all brain regions homogeneously. Frontoparietal and sensory association cortices show the regionally specific effect, while other regions do not. The ME-fMRI approach allowed us to take a data-driven approach to finding those brain regions that showed age-related changes. The estimation of seed-based connectivity among these regions, after controlling for the effect of the number of BOLD components, showed that these regions form a functional network with each other. Among regions with age-related change in BOLD component overlap, the cerebellum was found to be the most highly connected to the others, as well as to have the greatest change in connectivity with age.
The pattern of change in component number versus age is consistent with expectations based on age-related structural brain change. Brain development is a complex process that continues into early adulthood. After a fourfold increase in brain volume from birth to about 5 years of age, gross brain morphology stabilizes, with only an ∼10% increase in total volume up to adulthood. In adulthood, up to the third decade of life, microstructural change is the main mechanism of brain development (Raz, 2004; Lebel et al., 2008). Microstructural changes include increasing myelination of white matter and synaptic pruning. The inter-relationship of brain structure and function leads to an expectation that the functional organization of the brain would also change substantially due to brain development. In this study, we show evidence of such change in terms of functional organization as well as a key statistical characteristic of functional organization called BOLD dimensionality, or component count. Recent work has already shown that network segregation is one of the most robust and reproducible findings in developmental network neuroscience (Fair et al., 2007; Supekar et al., 2009; Dosenbach et al., 2010; Baum et al., 2017). Studies have consistently shown that short-distance connectivity decreases while long-distance connectivity increases with subject age. This pattern in fact concurs with the decrease in BOLD dimensionality with age, as more regionally specific networks that dominate at younger ages integrate into anatomically distributed and functionally distinct networks over neurodevelopment. The present work specially demonstrates that this pattern of the integration in the developing connectome also manifests with a change in the global and regional statistical characteristics of BOLD signal itself.
Virtually all fMRI studies on functional connectivity to date are based on the assumption that the number of signal components in resting-state data is a random effect (van den Heuvel and Hulshoff Pol, 2010). This is evidenced by the current standard of estimating connectivity based on a version of the Fisher transform that drops its SE term. This term would normally be used to control for dataset-level variability in the number of components in the time series data overall. In a prior study, we showed that this assumption is not valid from a methodological standpoint. Specifically, head motion reduces the number of BOLD components by reducing acquisition sensitivity (Kundu et al., 2013). Here we provide the first evidence that the assumption is also invalid on a neurobiological basis because the number of signal components varies systematically with age.
It is a basic principle for how correlation is interpreted, from the canonical form of the Fisher R-Z transform, that an accurate estimate of the statistical significance of correlation depends on a proper factoring of the number of independent components in the data. Given the magnitude of the effect of age-related change in the number of BOLD independent components, up to a factor of three between adolescence and middle age, the consideration of BOLD component counts may be important for interpreting correlation-based connectivity estimates along the lifespan and in comparisons involving disease. The ongoing European Autism Interventions-Multicenter Study (EU-AIMS) study is acquiring multiecho resting-state fMRI in healthy individuals and patients with autism across the age range, and that study will permit further assessment of this effect in health and disease (Murphy and Spooren, 2012).
The present analysis incorporated subject age as a correlate of BOLD component number, which in turn led to the identification of brain regions with similar trajectories of functional change with age. These regions included dorsolateral and dorsomedial prefrontal cortex. The developmental sensitivity of these regions found here is consistent with existing findings of their prolonged developmental trajectory. The dorsal prefrontal cortex mediates the most complex, higher-level aspects of cognition (Petrides, 2000; Koechlin et al., 2003; Koechlin, 2011; Passingham and Wise, 2012). Activity within the dorsolateral prefrontal cortex relates mostly to monitoring and sequencing information in working memory (Petrides, 1995, 2005; Owen et al., 1998; Amiez and Petrides, 2007). It also relates to the hierarchical organization and sequencing of other cognitive operations and to performance monitoring (Duncan and Owen, 2000; Koechlin et al., 2003; Duncan, 2010). Activity in the dorsomedial prefrontal cortex is associated with social cognition both when making judgments about others' mental state or intentions (Amodio and Frith, 2006; Lieberman, 2007; Behrens et al., 2008; Krienen et al., 2010) and while reflecting on one's own self, beliefs, intentions, and actions (Amodio and Frith, 2006; Brass and Haggard, 2007; Desmet et al., 2011). Importantly, these functions are considered to reach maturity only in adulthood. This study also brings into focus age-related changes in the precuneus, which only recently has become prominent in the study of brain development (Dosenbach et al., 2010). The precuneus is a central node of the default-mode network, which is involved in self-referential processing. However, further work is needed to examine the behavioral significance of age-related regional variability in BOLD components.
The results also highlight a potentially critical role of the cerebellum in age-related brain connectivity. The high signal-to-noise ratio (SNR) of the cerebellum after ME-ICA denoising enabled the detection of substantial age-related changes in BOLD dimensionality and functional connectivity in this region. Historically, the cerebellum has been viewed mostly as a motor control region, but more recently there is increased recognition of its role in the regulation of autonomic function and cognition (Reis and Golanov, 1997; Parsons et al., 2000; Craig, 2002; Singer et al., 2004) and in affective processing (Schmahmann and Caplan, 2006). Cerebellar lesions give rise to a constellation of cognitive and affective abnormalities (i.e., the cerebellar cognitive affective syndrome; Schmahmann and Sherman, 1998). From a developmental standpoint, the cerebellum is reported to have low heritability in morphology, indicating a greater role of environment in its neurodevelopment (Gilmore et al., 2010). Functional neuroimaging studies report that regions of the cerebellar cortex most responsive to cognitive demands—namely lobules VI and VII (Stoodley and Schmahmann, 2009; Stoodley et al., 2012)—interact closely with prefrontal and parietal association cortices (Habas et al., 2009; Buckner, 2013). These prior findings agree with the present results that show augmentation in the functional relationship between prefrontal regions and the cerebellum during brain maturation between adolescence and middle adulthood, suggesting the need for further study for this critical brain region (Raz et al., 1992).
The limitations of the present study include a relatively small sample size compared with recent large-scale or multisite neurodevelopmental studies. Larger studies, especially those using ME-fMRI, will be better suited to evaluate the relationship of BOLD component number and gender and cognitive variables, after factoring out the apparently large effect from subject age. To better understand the underlying neurobiological processes involved in the global and regional changes in BOLD component number, imaging of brain microstructure may be necessary. This would involve the additional acquisition of diffusion weighted imaging that is sensitive to processes like myelination or synaptic pruning, such as neurite orientation dispersion and distribution imaging (Zhang et al., 2012). To elucidate whether functional and microstructural changes also correlate to changes in brain metabolites, magnetic resonance spectroscopy (MRS) will be needed, possibly with a focus on regions that show change in BOLD component number. Also important is the consideration of how changes in BOLD component number relate to age-related changes in variability in regional activation during a task, which has already been associated with subject age (Garrett et al., 2011). Quantitative measurements of cerebral blood flow and metabolism alongside ME-ICA may also be informative as to how age-related hemodynamics may relate to BOLD component number. This is possible given a recent MRI sequence that acquires multiecho multiband EPI with simultaneous arterial spin labeling, and has been shown to be compatible with ME-ICA (Cohen et al., 2017). Altogether, future studies could use advanced MRI sequences in imaging larger cohorts with in-depth behavioral data to better characterize the neurobiological changes that age-related change in BOLD component number reflects.
Footnotes
- Correspondence should be addressed to Prantik Kundu, Section on Advanced Functional Neuroimaging, Brain Imaging Center, Icahn School of Medicine at Mount Sinai, New York, NY 10029. prantikk{at}gmail.com or prantik.kundu{at}mssm.edu
