A total of 227 volunteers (aged 18–76 years, male/female = 99/128) completed the study and were recruited into the study through the local media advertisement in the Stockholm region. All participants were right-handed, and native Swedish speakers with normal or corrected-to-normal vision. They all reported being free of a history of neurological, psychiatric, and cardiovascular diseases. None of the participants reported any use of psychotropic drugs. Each subject signed informed consent before completing the magnetic resonance imaging (MRI) examination protocol. They were financially compensated for their participation. The regional ethics committee approved the study protocol 2014/1982-31/1, which was conducted in line with the declaration of Helsinki.

The MRI data acquisition was conducted on a whole-body 3T clinical MRI scanner (Magnetom Trio, Siemens Medical Solutions, Erlangen, Germany) equipped with a 32-channel phased-array receiving head coil. All data was acquired at Karolinska University Hospital, Huddinge, Stockholm, between noon and 5:00 PM. The MRI data acquisition protocol included the following scanning sessions: (1) 3-plane localizer; (2) Conventional clinical MRI scans including 3D T1-weighted MPRAGE, T2 and FLAIR scans; and (3) A session of 375 s long R-fMRI measurements. The main acquisition parameters for the R-fMRI data included the following: TE/TR 35/2,500 ms, flip angle = 90°, 34 slices of 3.5 mm thick, FOV = 225 mm, matrix size = 76 × 76, data acquisition acceleration with GRAPPA parallel imaging method (iPAT = 2), and 150 dynamic timeframes. The T1-weighted MPRAGE images used for co-registration with functional images were acquired with the following parameters: TR = 1,900 ms, TE = 2.52 ms, FA = 9 degrees, FOV = 256, voxel size 1 mm × 1 mm × 1 mm. The acquisition parameters for the FLAIR image were the following: TE/TR = 89/9,000 ms, flip angle = 130°; inversion time (TI) = 2,500 ms, slice thickness = 4.0 mm, and FOV = 199 × 220 mm. An experienced radiologist inspected both the FLAIR and T1-weighted images for potential signs of neuropathology.

We used foam patting to fix each subject’s head carefully in the head coil to reduce involuntary head motions. During the R-fMRI data acquisition the participants were instructed to focus their sight on a white cross in black background projected on a screen installed in front of their eyes. The subjects were also instructed to not think about anything particular during the R-fMRI session.

The R-fMRI datasets underwent a preprocessing procedure, which has been described elsewhere in details (Li et al., 2021) and was performed with AFNI (Version Debian-16.2.07∼dfsg.1-3∼nd14.04+1, http://afni.nimh.nih.gov/afni) and FSL¹ programs with a bash wrapper shell (Wang and Li, 2013, 2015). After temporal de-spiking, six-parameter rigid body image registration was performed for motion correction. The average volume for each motion-corrected time series was used to generate a brain mask to minimize the inclusion of the extra-cerebral tissues. Spatial normalization to the standard MNI template was performed using a 12-parameter affine transformation and a mutual-information cost function. During the affine transformation the imaging data were also re-sampled to isotropic resolution using a Gaussian kernel with 4 mm full width at half maximum (FWHM). The co-registered average image volume for the cohort has 28,146 non-zero voxels inside the brain and was used to generate the average brain mask for the preprocessed whole-brain R-fMRI data with 4 mm spatial resolution. Nuisance signal removal was performed by voxel-wise regression using 14 regressors based on the motion correction parameters, average signal of the ventricles and their 1st order derivatives. After baseline trend removal up to the third order polynomial, effective band-pass filtering was performed using low-pass filtering at 0.08 Hz. Local Gaussian smoothing up to FWHM = 4 mm was performed using an eroded gray matter mask (Wang and Li, 2015).

Pearson’s correlation coefficients (CC) were computed between the time courses of all pairs of voxels inside the brain, leading to a whole-brain functional connectivity matrix for each subject. This computation was performed for all voxels located within the brain mask, which was generated by overlapping the registered brains of all participants. This brain mask contained 28,146 voxels and each voxel inside the brain was used as the seed voxel in turn. Therefore, the size of the CC matrix size is 28,146 × 28,146. Each row or column of the CC matrix corresponds to the CC image volume for the seed voxel with the rest of the brain. That is the connectivity map for the seed voxel. As schematically illustrated in Figure 1, based on the CC histogram for each row of the matrix we derived the following two types of threshold-free voxel-wise RFC metrics: the CSI and CDI. As we are interested in systematically investigating all relevant synchronized activities in the whole brain, we quantify the negative and positive portions of the CC histogram separately to avoid information cancelation, sensitivity reduction, and statistical interference. From here on, the subscripts “N” and “P” are used to indicate the negative and positive portions of the RFC metrics, respectively. The metrics without subscripts refer to the mixed measures without distinction of the negative and position correlations.

As shown in Figure 1, the voxel value for the CSI_P and CSI_N are defined as the averages of the positives and negatives in each row of the CC matrix, respectively. That is

Where CC_row refers to a row in the CC matrix. np and nn refer to the number of positive and negative correlation coefficients in a row of the CC matrix, respectively. The voxel values for CDI are defined as the convolution between the CC histogram and a kernel function. That is

The CDI_P and CDI_N correspond to the positive and negative portions of the convolution defined in Eq. (3), respectively. To facilitate statistical comparison it is useful to transform the raw RFC metrics into standard Z-score using the following formula:

Where μ and σ are the mean and standard deviation of the corresponding RFC metrics, respectively. For optimization of the CDI sensitivity, we investigated 6 different kernel functions including

where x⊂ [−1,1] corresponds to the interval of the correlation coefficients. The kernels are also graphically depicted in Figure 2. The kernel should weight the higher correlation coefficients more than the lower ones. The widely used threshold approach can be considered as the case of the square-well kernel function k₆. For illustration, an arbitrary threshold of 0.3 was used here. The CSI metrics can also be considered as a special case of CDI corresponding to a kernel of the sign function.

To investigate if and how the RFC metrics are influenced by heathy aging for the studied cohort we performed voxel-wise linear regression analyses of the CSI and CDI metrics versus the subject’s age, while gender was treated as a covariate by using the AFNI program 3dRegAna to extract the regression parameter β and linear coefficient r. The statistical significance was assessed by using a two-step approach. Firstly, we imposed a voxel-wise threshold p < 0.001 (uncorrected corresponding t-score ≥ 3.34) to form the initial cluster candidates. Secondly, we performed permutation simulations without assuming a particular form of probability distribution for the voxel values in the statistic images to identify the brain ROI out of the initially detected clusters at family-wise error rate p ≤ 0.05. Using the detected ROIs as masks, the mean values of the RFC metrics for each ROI were evaluated and plotted against the subjects’ age. Besides linear regression analysis with age, we performed also verification using two-sample t-test between the young and elderly subgroups. For this, we selected all subjects aged 18–30 years as the young subgroup (n = 124, males/females = 51/73), and all subjects aged 64–76 years as the elderly subgroup (n = 76, males/females = 35/41). To keep sufficient age gap between the young and elderly subgroups the remaining 27 subjects in the age range of 31–63 years old were excluded from the t-test. In the selection of subgroups we attempted to minimize the number of excluded subjects with intermediate ages, maximize the age gap between the subgroups, and keep similar number of subjects and age ranges. It should be emphasized that all 227 subjects were included in the regression analysis.

The CC histogram for each seed voxel in the brain is dependent on its location in the brain (see Supplementary Materials). Figure 3 shows the average CC histogram of the cohort for a seed voxel in the PCC as illustrated by the cross in Figure 1. The histogram is somewhat asymmetric and shifted toward the positive side. This is quite typical at least for voxels within gray matter. Selecting different threshold values along the histogram allows us to examine the RFC networks of different connection strengths associated with the selected seed voxel. As shown in Figure 3, at high negative threshold (Figures 3A,B) we observe the DMN. At low negative threshold, we observe its association with cerebral spinal fluid (CSF) space and white matter (Figure 3D). At moderately high positive threshold, the PCC voxel is not only a part of the DMN, but also connected to most of the cortical gray matter (Figures 3E,F). At high positive threshold, the PCC voxel is associated with the posterior portion of the DMN and the visual cortex (Figure 3C). The visual cortex looks relatively bright at high threshold indicating a relatively high number of voxels are associated with the visual network or voxels within the visual cortex are associated with each other at high threshold criterion. The idea of the QDA framework is to avoid the arbitrary threshold and optimize the contribution of meaningful informatics to the quantitative RFC metrics.

Figure 4 shows an axial slice of CDI_P and CDI_N images for a typical R-fMRI dataset (from a 36 year old male subject). Multiple brain regions depict high CDI_P including the bilateral mPFC, superior and middle temporal gyri (MTG), inferior and superior parietal lobule, precuneus and PCC. These regions have been described as RFC hubs implying their important role in neural signaling and communication across the brain (Buckner et al., 2009; Tomasi and Volkow, 2011). On the other hand, the PCC, insula cortex. White matter and CSF regions have high CDI_N metric. The contrast and intensity variations across each row in Figure 4 demonstrate that selection of the kernel function can optimize the contrast and signal-to-noise ratio of the CDI metrics.

The linear regression results for the CSI, CSI_N, and CSI_P data versus subjects’ age are summarized in Figure 5 and Table 1. The corresponding results for the CDI, CDI_N, and CDI_P are shown in Figure 6 and Table 2. The CSI metric without separation of the negative and positive correlations shows decline of the functional connectivity strength with age in the superior and middle prefrontal gyrus (MFG) and increase of connectivity strength in the precuneus and right inferior parietal lobule (r-IPL). The more specifically defined CSI_N and CSI_P metrics are more sensitive to the adult age effect and the detected brain volumes with significant aging effect are nearly tripled compared with that for the CSI metric. With CSI_N and CSI_P we also observe a more intricate pattern of change with the adult age, which are summarized as follows:

To study the specific connectivity associated with the two ROIs defined by the overlap between the CSI_N and CSI_P metrics, we computed Pearson’s correlation maps for the time courses of the seeds as defined by the overlapping ROIs depicted in Figure 7A. As expected, the associated RFC network for the PCC ROI is obviously the well-known DMN and include 4 negatively correlated brain regions, which are the bilateral IPL and insula cortices (see Figure 7B). On the other hand, the associated RFC network for the insula ROI includes the PCC and bilateral precuneus as the negatively correlated brain regions (Figure 7C) Figure 8. Shows the anti-correlated brain regions between the above 2 RFC networks as obtained by multiplying the two correlation maps with each other and applying negative threshold at CC ≤ (−0.5). The mutually inclusive anti-correlation between the PCC and the right insular cortex explains why both CSI_P and CSI_N metrics in these regions depict declines with the adult age.

Figure 9 shows the ROI average of the CSI_N and CSI_P metrics in the PCC and right insula cortex as a function of the subject’s age. With normal aging, both the CSI_P and CSI_N are reduced in these brain regions (overlap shown in Figure 7A). Therefore, the PCC and right insula are particularly sensitive to the adult age effect. However, the aging effect is barely detectable by the unseparated CSI metric (see Table 1).

As expected, the CDI_P and CDI_N metrics derived by using the different kernels differ in their sensitivity in detecting the adult age effect. Figure 10A shows the detected brain volumes where the CDI_P and CDI_N metrics are significantly associated with the adult age. The sensitivity difference of the kernels is also manifested in the regression parameter β which are detailed in Table 3 and Figure 10B. To compare the similarity of the detected aging effects among the CDI metrics of different kernels, we assessed the joint overlapping brain regions detected by the different CDI_N and CDI_P metrics of different kernels. The observed overall trends of RFC enhancement or decline with age are quite similar. The joint overlapping volumes for the CDI_P and CDI_N metrics of different kernels are 733 and 671 voxels, respectively. Moreover, there is also a reasonable anatomic consistency between the results of the connectivity strength metrics and connectivity density metrics. As detailed in Tables 2, 3, the anatomical locations of the joint overlapping regions for the different CDI_P metrics match those for the 3 largest ROIs identified by the CSI_P results (see Table 1). Similarly, the brain regions of the joint overlapping for the different CDI_N metrics are largely the same as those identified by the CSI_N data (see Table 1). However, it should be noted that the β parameters for the CDI_N and CSI_N have opposite signs even through the trend of change with the adult age is the same. This is because the negative connectivity strength (CSI_N) is negative in nature, while the connectivity density corresponding to the negative correlation (CDI_N) is always positive. Therefore, the enhancement of the negative connectivity strength (CSI_N) with age (for example in the sensorimotor network) corresponds to a negative β, while the connectivity density result corresponds to a positive β value.

Age is an important risk factor for declines of neural cognitive functions and pathology of neurodegenerative diseases. It is also a complex metric difficult to precisely interpret the involved physiology. Healthy individuals of similar age may have quite different vascular and brain-health status. It follows that age is not a single strongest predictor for the RFC in the brain. This is likely to be the reason why the linear regressions of the RFC metrics with the adult age depict substantial scatters and relative low correlation coefficients. The impact of the potential confounds and pre-processing strategies that can mitigate them have been extensively investigated in the published literature (Bluhm et al., 2008; Biswal et al., 2010; Weissman-Fogel et al., 2010; Zuo et al., 2010; Dennis and Thompson, 2014; Alarcon et al., 2015; Zhang C. et al., 2016; Geerligs et al., 2017; Hussein et al., 2020). Here we focus on comparing our findings in the context of documented literature results, particularly the adult age effect in the DMN, dorsal attention network (DAN), sensorimotor network, and subcortical brain regions.

With QDA, we found support for RFC decline with advancing adult age in multiple brain regions of the DMN and DAN, including superior and MFG, PCC, MTG, and IPL. Age-related RFC decrements in the DMN and DAN have previously been reported in numerous R-fMRI studies using ROI and ICA based analysis (Buckner et al., 2008; Damoiseaux et al., 2008; Dennis and Thompson, 2014; Scheinost et al., 2015; Luo et al., 2020). Our findings regarding to the RFC changes in the DMN are overall in agreement with previous reported results (Damoiseaux et al., 2008, 2016; Koch et al., 2010; Ystad et al., 2010; Williams, 2013; Lu et al., 2014; Persson et al., 2014, 2015; Salami et al., 2014). Besides the DMN and DAN, normal aging was associated with RFC increase in the sensorimotor, subcortical network, and para-hippocampal cortex. This has also been reported previously (Persson et al., 2015, 2016; Damoiseaux et al., 2016; Geerligs et al., 2017; Hussein et al., 2020; Luo et al., 2020). We didn’t find significant age-related RFC declines in precuneus and specific sub-regions of the hippocampal cortex as reported in previous studies (Salami et al., 2014; Damoiseaux et al., 2016). Since we assessed the negative and positive correlation separately, this may allow us to detect more intricate age-related RFC changes in the brain. To illustrate this point, we analyzed further the 3 ROIs with significant correlation between the CSI and the subject’s age. As shown in Tables 1, 4 and Figure 11, the detected ROI in the precuneus depicted significant positive linear correlation between CSI and the subject’s age (β = 9.50 × 10^–3, r = 0.459), even though the CSI_P and CSI_N in the same ROI showed only a slight (not significant) increment and decrement with age, respectively, i.e., contribution from a low-significant CSI_P increment and a non-significant CSI_N decrement resulted in a highly significant increment trend in the CSI metric. With the same line of reasoning, we can explain why the MFG ROI detected by the CSI metric is much smaller than that detected by the CSI_P metric, because the decremental trend in the CSI_P metric was partially canceled by the CSI_N contribution. This can also explain why we didn’t detect significant CSI decrement with the adult age in the PCC and R-insula, because both the CSI_P and CSI_N metrics exhibited significant decremental trends with age and their contributions annulled each other. Therefore, it is important to pay attention to the precise definition of the RFC when comparing the results of different studies.

Both CSI_P (r = 0.506, see Table 1) and CDI_P (r = 0.577, see Table 2) showed age-related enhancement Caudate/putamen and the association are quite strong. This finding based on QDA approach are consistent with previous reports (Manza et al., 2015; Rieckmann et al., 2018) from ROI-based studies aimed to investigate the aging effect on specific functional connectivity of in the striatum-cortical system. It is well known that the striatum-somatomotor connection is primarily associated with motor performance, especially the “automatic” performance of already learned movements. It has been reported that posterior putamen and pallidum decreases in connectivity to left somatomotor cortex with age (Manza et al., 2015). This provides a reasonable explanation for the commonly observed motor deficits as in elderly subjects. There is also growing evidence support the notion that striatum-cortical connectivity is potentially important for memory function at older age. Intriguingly, this is related to the enhanced RFC with age in the putamen/caudate. Several studies have reported that increased striatum functional connectivity in older adults typically reflects less negative connectivity between two regions belonging to different networks, and the increased connectivity is often negatively associated with cognitive performance (Rieckmann et al., 2018). Therefore, better understanding the RFC change with adult age in the striatum-cortical system can be potentially useful for assessing motor as well cognitive functions in elderly subjects.

The QDA framework proposed in the study is a voxel-wise and data-driven approach. It has the following two unique features: (1) It can avoid confounding caused by the cancelation of the negative and positive correlations by assessing the negative and positive portions of the CC histogram separately; (2) It derives different RFC metrics based on the connectivity strength and density by utilizing the concept of convolutions with different kernels. The metrics weight all the correlations of a given voxel with the rest of the brain according to the amplitudes of the correlation coefficients and disregard the anatomical distance between the correlation pairs. This permits a comprehensive characterization of the intrinsic activities of each voxel without the use of an arbitrary threshold. The QDA approach can encapsulate the widely used threshold approach as a special case of the square-well kernel function. The widely used degree centrality corresponds precisely this square-well kernel situation which adopts a somewhat arbitrary threshold and every connection above the threshold are weighted equally. Even the CSI metrics can be encapsulated under the convolution concept for a special kernel of the sign function. This provides a unified view for RFC and can facilitate its further optimization. The QDA framework uses the time course of each voxel within the brain as the seed reference to compute voxel-wise whole-brain correlational coefficient matrix. For whole-brain R-fMRI data acquired at 4 mm spatial resolution, the correlational coefficient matrix is in the order of 10⁵ and is currently not practical for direct visualization and statistic assessment. Particularly, when data are acquired with higher spatial resolution, e.g., 2 mm, the matrix size is increased by 8 × 8 times. Therefore, for data reduction, we derived two types of voxel-wise RFC metrics from the correlational coefficient matrix without the need for specifying a particular seed, threshold, model, or mode. As their names indicated, CSI and CDI are aimed to capture the local (voxel) connectivity strength and density with rest of brain, respectively. The QDA metrics can assess the general connectivity with the rest of the brain without specifying a specific path or network. The QDA method does not highlight the specific connectivity changes between selected brain regions. The precise neural correlate of R-fMRI signal is currently not well understood (Hyder and Rothman, 2010). However, it is reasonable to assume that the R-fMRI signal fluctuations indirectly reflect the slow modulations of neuronal activities at rest. Furthermore, the sigmoid function has been widely accepted as the logistic function of neuronal activation instead of a square-well. With current status of knowledge, we cannot identify a convolution kernel to reveal a particular feature of the neurophysiology. However, we can attempt to optimize the kernel to reduce bias and improve sensitivity of the RFC metrics to pathophysiological changes.

The current results based on the QDA framework should be interpreted in the context of some technical and biological limitations. Firstly, at a TR of 2,500 ms, the cardiac and respiratory fluctuation effects might be aliased into the low frequency R-fMR signal fluctuations. The regression up-to the 1st order derivative of the head motions and lowpass filtering could not eliminate the effects of these physiological noises (Muschelli et al., 2014; Pruim et al., 2015; Bright et al., 2017; Parkes et al., 2018). Thus, these aliasing effects could reduce the specificity of the RFC metrics, or even might further confound the detected RFC differences between the young and elderly sub-groups. With the more up-to-date acquisition techniques, such as multi-band simultaneous acquisition of multiple slices and compress-sensing with high under-sampling factor, it is possible to use a shorter TR (e.g., 500 ms) and higher spatial resolution for the data acquisition. Therefore, these physiological effects may be further mitigated.

Secondly, the resting state is associated with spontaneous thoughts and cognitive processing, we cannot exclude the possibility that differences in spontaneous thoughts may exist between the young and elderly subjects (Wu et al., 2007). However, considering the overall consistency of our results with the previous studies, particularly the results from the longitudinal studies (Fjell et al., 2014; Ng et al., 2016; Staffaroni et al., 2018; Li Q. et al., 2020), it is unlikely that these differences have major influence on our findings. These initial findings encourage the future use of QDA as a tool to analyze longitudinal R-fMRI data aimed to develop a comprehensive understanding of age- or pathology-related brain functional changes.

Thirdly, the generalizability, or external validity issue should be considered. This is due to the non-random recruitment procedures and relying on a sample of convenience. The sample size used in this study (N = 227) is moderate, includes unbalanced young and elderly subgroups reflecting the difficulties to recruit elderly healthy subjects. The ages of the participants range from young to old adulthood (reflecting the age of participants in most neuroimaging studies). The age-related RFC differences observed in this study were relatively small but quite robust. However, the results from this cross-sectional study of the cohort cannot distinguish whether the RFC changes in the brain regions are due to gradual changes throughout the adulthood or a more sudden change at later stage in life.

Fourthly, the R-fMRI data were acquired under open-eye condition. Recent studies indicate that opening versus closing eyes at resting-state results in RFC difference between V1 with DMN and SNs (Costumero et al., 2020). This may explain why we did not detect significant RFC change with age in the visual cortex.

As discussed above negative correlation is an important fraction of the CC histogram irrespective of the tissue type and anatomical location of the voxel in question. In the published literature, there is also a rapid growing interest in studying the negative correlations between the voxels (Fox et al., 2009; Weissenbacher et al., 2009; Bianciardi et al., 2011; Schwarz and McGonigle, 2011; Gonzalez-Castillo et al., 2012; Gopinath et al., 2015; Liu et al., 2015; Spreng et al., 2016; Chen et al., 2020). It is clear that the negative portion of the CC histogram is more dominant for voxels in CSF (Gruszecki et al., 2018) and white matter (McColgan et al., 2017; Gore et al., 2019; Wu et al., 2019; Li M. et al., 2020). However, the negative portion cannot be ignored even for voxels in the gray matter (see Supplementary Materials). To avoid confound caused by inappropriate preprocessing pipelines, we have carefully tested and updated our preprocessing pipeline. We did not implement the global signal regression (GSR) which removes the mean signal averaged over the entire brain. GSR removal via linear regression is one of the most controversial procedures in the analysis of R-fMRI data (Fox et al., 2009; Weissenbacher et al., 2009). On one hand, the global mean signal contains variance associated with respiratory, scanner-, and motion-related artifacts. Its removal by GSR can improve various quality control metrics, which enhances the anatomical specificity of RFC networks, and increase the explained behavioral variance. On the other hand, GSR alters the distribution of regional signal correlations in the brain, can induce artefactual anti-correlation patterns, may remove real neural signal, and can distort RFC metrics. The brain masked “global signal” is usually misunderstood, because it is not “global” and its variance contains dominant contributions from different domains of the voxels with temporally coherent signal variation.

To limit the study in a reasonable scope, in the discussion of the adult age effect on RFC we focused on gray matter and did not discuss white matter and CSF related issues. However, it should be pointed out that aging effects in white matter (McColgan et al., 2017; Gore et al., 2019; Wu et al., 2019; Li M. et al., 2020) and CSF (Sakka et al., 2011; Gruszecki et al., 2018) are also worth exploring. There is indeed a rapid growing interest in these arenas in published literature (Sakka et al., 2011; McColgan et al., 2017; Gruszecki et al., 2018; Gore et al., 2019; Wu et al., 2019; Li M. et al., 2020), particularly in the context of the age effect for the glymphatic system.