The dataset was acquired from the NKI/Rockland Sample database (http://fcon_1000.projects.nitrc.org/indi/pro/nki.html), including the imaging data recorded from 207 participants aged 4–85 years old. All participants underwent semi-structured diagnostic psychiatric interviews and completed a battery of psychiatric, cognitive, and behavioral assessments. All participants had no history of taking medicine, psychiatric illness, or any brain surgery that might affect their central neural system; they had no serious unstable illnesses requiring medication or hospitalization, including hepatic, renal, gastroenterology, respiratory, cardiovascular (including ischemic heart disease), endocrinologic, neurologic, immunologic, or hematologic disease; were not pregnant, breastfeeding, or had other contraindications to magnetic resonance imaging (MRI) at the time of the study. All participants signed a written informed consent form according to the Helsinki Declaration and satisfied the criteria for MRI. The Ethics Committee of the Department of Psychology at Renmin University of China authorized the current analysis procedure.

All participants were scanned using a 3-Tesla magnetic resonance scanner (B15, Siemens AG, Munich, Germany). They were required to remain quiet, relaxed, and awake during the ten-minute resting state scanning session with their eyes closed and heads fixed using comfortable sponge fixators. Functional images were acquired using a single-shot T2-weighted gradient-echo-planar imaging (EPI) sequence. The scanning parameters were as follows: repetition time = 2500 ms, echo time = 30 ms, flip angle = 90°, matrix size = 64 $\times$ 64, field of view (FOV) = 216 mm $\times$ 216 mm, slice number = 38 (axial interleave acquisition), slice thickness = 3.0 mm, voxel size = 3.0 mm $\times$ 3.0 mm $\times$ 3.0 mm.

Raw fMRI data were preprocessed using the DPABI toolbox [27] in MATLAB 2013a. The First ten volumes were removed to remove potentially unstable signals. The image preprocessing procedures include slice-timing, realignment to correct head motion, spatial normalizing of the individuals’ images to a standardized spatial Montreal Neurological Institute (MNI) space [28], resampling to the voxel size of 3.0 mm $\times$ 3.0 mm $\times$ 3.0 mm, and spatial smoothing using an 8 mm full-width at half maximum Gaussian kernel to reduce spatial noise [29]. Several optimized denoising procedures were used to reduce spurious temporal signals further. We (1) implemented independent components analysis to remove nuisance components [30] using the GIFT v4.0 toolbox (http://mialab.mrn.org/software/gift); (2) scrubbing procedure to detect and correct the time points with a frame-wise-displacement larger than 0.5 mm [31] and replace them with temporal interpolation; (3) regressing out the global noise, white matter signal, cerebrospinal fluid signal, and the head motion parameters as nuisances; we performed global signal removal because the common noise embedded in the global signal may cause severe contamination by false-positive correlation in traditional FC analysis and severely degrade the estimation of TI and GC; (4) detrending to remove signal drift; and (5) bandpass (0.01–0.08 Hz) filtering to acquire a stable low-fluctuating blood oxygenation level dependent (BOLD) signal [10, 32].

For the total independence (TI) and GC (Granger causality) estimations, we used the hemodynamic response function (HRF) toolbox (http://users.ugent.be/ dmarinaz/code.html) to perform deconvolution on the preprocessed fMRI data and retrieve the signal changes at the neural level to alleviate potential confounds caused by HRF inhomogeneity [33, 34].

We use the Neurosynth online meta-analysis tool (https://www.neurosynth.org) to define the ROIs of four functional networks closely related to cognitive control, including CON, DAN, CEN, and DMN, and calculate the connectivity within CON and between CON and the rest four networks to elucidate how CON communicate with the cognitive-control-related networks (For the convenience of description, we summarized these networks as cognitive control-related functional networks (CCRNs) in the current study). Specifically, the results of the meta-analysis regarding the keyword “attention” in 1831 studies provided a well-proposed spatial pattern (uniformity test, p 0.01, False discovery rate (FDR) corrected) of the task-positive regions of CON, DAN, and CEN. On the other hand, the results regarding the keyword “default mode” in 777 studies outlined the classical spatial pattern (uniformity test, p 0.01, FDR corrected) of the DMN, which is thought to be a source of internal interference to goal-directed cognitive control [35]. CON ROIs included dACC and bilateral anterior insula (AI) [36]. DAN ROIs included bilateral frontal eye field (FEF) and bilateral intraparietal sulcus (IPS) [37]. CEN ROIs included bilateral precentral gyrus/frontal middle gyrus/subfrontal gyrus/dorsolateral prefrontal (PreCentral/MFG/IFG/DLPFC, collectively referred to as DLPFC in the current study [14]. DMN regions included the medial prefrontal cortex (MPFC), posterior cingulate cortex (PCC), bilateral inferior parietal lobules (IPL), and lateral temporal lobe (LTC) [38]. Please see Fig. 1 for the spatial locations of the ROIs.

Fig. 1.

The cognitive-control-related networks were defined based on a meta-analysis. CON, including dACC, rAI, lAI; FPN, including rFEF, lFEF, rIPS, lIPS; DMN, including MPFC, PCC, rIPL, lIPL, rLTC and lLTC; CEN, including bilateral PreCentral/MFG/IFG/DLPFC.

The seed time series representing each CON ROI’s activity was obtained by averaging all voxel time series in that CON ROI (dACC, right AI, and left AI, respectively). Then we calculate CC, TI, and GC between the CON ROI seed time series and each voxel of all four CCRNs (CON itself, DAN, CEN, and DMN) to yield a map that delineates the interaction between the CON ROI and other regions of CCRNs.

A common approach to demonstrate the topography is providing the thresholded T-map of an age group using a one-sample t-test. However, the sample size affects the p-value threshold selection. For age groups with different numbers of participants, setting one p-value threshold for all three groups is not viable for comparing their topography. Therefore, we applied a ranked mapping approach to intuitively demonstrate the difference in connectivity maps between the three age groups. Specifically, we ranked the connectivity index values across all participants for each voxel. Then we summarized that in each group, the portion of participants fell into the upper half (stronger half) of all participants. For each group, a voxel was highlighted if the portion of strong participants it represented was larger than 50%. Namely, the map illustrated whether a voxel in a given age group is strong relative to the distribution of all participants.

We calculated the cross-correlation [39, 40] between time series to measure the FC between CON and the other four networks. The Cross-correlation (CC) is the most commonly used FC measurement describing the consistent collaboration or antagonism relationship between regions/networks. In the current study, the raw CC indices were Fisher-transformed to project the values on [-$\mathrm{\infty }$, $\mathrm{\infty }$]. CC measures the point-to-point linear relationship between two time series but neglects their non-zero lag temporal interdependence [32].

Unlike the traditional CC, total interdependence (TI) measures total temporal interdependence and accounts for the non-zero-lag relationship existing across different lags between time series (e.g., between x(i) and y(i+n)) [32]. TI was defined by Gelfand and Yaglom [41] as:

(1)$$T{I}_{x,y}=-\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\ln \left(1-C_{xy}^2(\lambda )\right)𝑑\lambda$$

where Cxy($\lambda$) is the coherence between the two random processes, x and y, at frequency f = $\lambda$/2$\pi$.

For a given sampling frequency fs, the equation can be converted into an implementable form as:

(2)$$T{I}_{x,y}=-\frac{2}{f_s}\sum _{i=1}^{N-1}\ln (1-C_{xy}{}^2(i\mathrm{\Delta }f))\mathrm{\Delta }f$$

where $\mathrm{\Delta }$f is the frequency resolution of the spectrum of the time series.

Therefore, by calculating TI, we can more intuitively and quantitatively prove the complete interdependence between two time series not captured by CC.

A directed network model would be meaningful to describe the information flow between CON and the other networks. GC analysis is an effective connectivity approach useful in directed network modelling. Let X and Y be two time series, if using both X’s and Y’s past, one can predict Y’s future better than merely using Y’s past, then we say that there is GC influence from X to Y [42]. The estimation of GC is realized using the multi-variable auto-regressive (MVAR) model [43, 44, 45].

Specifically, let the two stationary time series be denoted by X${}_t$ and Y${}_t$, each is independently represented by the following univariate autoregressive (AR) models:

(3)$$\begin{array}{c}\hfill X_t=\sum _{j=1}^{\mathrm{\infty }}{a}_{1j}{X}_{t-j}+{\epsilon }_{1t},\mathrm{var}\left({\epsilon }_{1t}\right)={\mathrm{\Sigma }}_1\hfill \\ \hfill Y_t=\sum _{j=1}^{\mathrm{\infty }}{d}_{1j}{Y}_{t-j}+{\eta }_{1t},\mathrm{var}\left({\eta }_{1t}\right)={\mathrm{\Gamma }}_1\hfill \end{array}$$

and jointly represented as the following bivariate AR model:

(4)$X_t={\displaystyle \sum _{j=1}^{\mathrm{\infty }}}{a}_{2j}{X}_{t-j}+{\displaystyle \sum _{j=1}^{\mathrm{\infty }}}{b}_{2j}{Y}_{t-j}+{\epsilon }_{2t}$
$Y_t={\displaystyle \sum _{j=1}^{\mathrm{\infty }}}{c}_{2j}{X}_{t-j}+{\displaystyle \sum _{j=1}^{\mathrm{\infty }}}{d}_{2j}{Y}_{t-j}+{\eta }_{2t}$

where $\mathrm{\Gamma }$ ${}_1$ denotes the variance of the residue of the univariate AR model fitting and $\mathrm{\Gamma }$ ${}_2$ the covariance of the residue matrix of the bivariate (or multivariate) AR model fitting. Then GC of X$\to$Y can be defined as:

(5)$$F_{X\to Y}=\ln \frac{{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_2}$$

According to the bayesian information criterion (BIC) and akaike information criterion (AIC) model order estimations and our previous studies [46, 47, 48, 49], the order of the autoregressive model (AR) model was set to 2 in this study. More details of the mathematical realizations were provided in our previous study [45]. GC measures directed information flow and has been widely used in directed network modeling based on multi-unit recordings, electroencephalogram (EEG), and fMRI, providing insightful results for complex interactions among large-scale networks.

We could generate connectivity maps of dACC-CCRNs, rAI-CCRNs, and lAI-CCRNs for each participant using the above connectivity measurements. Then for each connectivity index and each CON seed ROI (say region S), we calculate the correlation coefficient between each CCRNs voxel value and the age across participants to yield a map delineating the age-varying pattern of the connectivity between S and CCRNs. We start with performing the connectivity-age correlation across all participants. Further, from a lifespan perspective, it is expected that many biological indices vary in a nonmonotonic fashion against age, and taking all participants in one shot may overlook the detailed age-related pattern. It is reasonable to divide the participants according to their natural life stages and perform the correlation test in each stage. Therefore, we divided the 207 participants into three groups: Group [4–20 years old] (N = 56), Group [21–59 years old] (N = 119), and Group [60–85 years old] (N = 32). The correlation analysis in Group [4–20 years old] helped to assess how the CON-CCRNs interaction develops before adulthood, the analysis in Group [21–59 years old] helped to assess the changing of CON-CCRNs interaction during young and middle-aged adulthood, and the analysis in Group [60–85 years old] helped to assess the changing of CON-CCRNs interaction during the elder adulthood.

The connectivity-age correlation test included the following steps. First, we performed a z-transform to normalize CC, TI, and GC maps. Second, for each CCRNs voxel, the correlation between each indicator and the age of participants was calculated using an age-ranked grouping method [50, 51, 52, 53]. Specifically, (1) participants in the current age group were arranged in an age-descending order. (2) Then the analysis moved from the youngest to the eldest participant and with an increment of 1 year, and at the ith step of age x(i), calculate the mean value of the connectivity indices I(i) and mean age xmean(i) of the participants fall into the span between x(i) and x(i)+3 years. (3) The correlation coefficients test between I and xmean was performed to yield the correlation coefficient R, significance value p, and the coefficient of determination R${}^2$.

In order to explore the difference across the three groups, we performed one-way ANOVA based on F-test for CC, TI, and GC mapping results, respectively. Statistical significance for each ANOVA was assessed at the p 0.05 level and corrected using False discovery rate (FDR) multiple comparison correction. Post hoc analysis based on the Scheffe Test was applied within the CCRN mask to reveal regions with significant changes in multiple indicators.

First, we use the ranked mapping approach to demonstrate the spatial difference of the connectivity map between the age groups. Noted that the classical FC was measured using cross-correlation defined on [–1, 1], negative CC values representing anti-correlated activity could rank low, but that did not mean a weak interaction. For example, CON and DMN show prominent anti-correlation during the resting state; whether the negative correlation between CON and DMN BOLDs represents a strong but antagonist interaction or merely extremely weak coupling remains debatable. The aforementioned ranked mapping approach does not fit the complicated issue of negative FC. Therefore, we only applied the ranked mapping approach to TI and GC (Fig. 2). The result showed that the TI and GC output in Group [4–20 years old] had weaker connectivity within CON, with lIPS, and with DMN, compared to the other two groups. For GC input, CON seemed to receive more influence from DMN and IPS but less from bilateral AI.

Fig. 2.

A summary of ranked seed map. The colored regions delineate the CCRN regions that showed relatively strong connectivity with CON seed ROIs (red, seeded by dACC; green, seeded by rAI; blue, seeded by lAI).

Second, we performed ANOVA to explore the difference across the groups. No significant difference was detected across the three groups (p 0.05, FDR correction).

Third, we performed the connectivity-age correlation across all participants. No connection indicators significantly correlate with age at a lifetime span were found (R${}^2$ $\ge$0.6, p 0.001).

The result showed that most of the significant connection indicators detected were positively correlated with age in this group. The detailed spatial map of the CCRNs regions showed significant (R${}^2$ $\ge$0.6, p 0.001, FDR corrected) correlation/anti-correlation, as shown in Figs. 3,4, and with detailed information listed in Table 1. The results demonstrated several prominent patterns of the development of CON-CCRNs interactions. First, within the CON, enhanced CC/TI among dACC and bilateral AI and enhanced GC of dACC$\to$lAI/rAI and rAI$\to$dACC were detected, indicating increased communication within the CON from childhood to adulthood. Second, most of the indicators detected between the CON and DAN were found enhanced, such as the CC of lAI-rIPS/lIPS and rAI-rIPS/rFEF, TI of dACC-rFEF, lAI-rIPS/lIPS, and rAI-rIPS, and GC of lIPS$\to$lAI. Only a few indicators regarding the CON and the left DAN nodes, including the CC of dACC-lIPS, TI of rAI-lEFE, and GC of rAI$\to$lIPS, showed attenuation. Third, all significant connectivity indicators detected between the CON and CEN, including the CC of dACC/rAI/lAI-rDLPFC and lAI-lDLPFC, TI of dACC-rDLPFC, and GC of dACC$\to$rDLPFC/lDLFPC and lAI-lDLPFC, showed enhancement against age. Fourth, intense enhanced connection indicators found between the CON and DMN, including the CC of dACC-PCC/lIPL/lLTC, rAI-PCC/lIPL/rIPL, and lAI-PCC/lLTC/lIPL/rIPL/DMPFC/VMPFC/ACC, and TI of dACC-PCC/lIPL, rAI-PCC/lLTC/lIPL/rIPL, lAI-PCC/lLTC/lIPL/rIPL. Only TI: lAI/VMPFC/ACC and GC VMPFC/ACC$\to$lAI and rAI$\to$PCC showed attenuation against age. Though we set up a significant threshold of R${}^2$ $\ge$0.6, most of the enhanced indicators we detected in Group [4–20 years old] had R${}^2$ $\ge$0.8.

Fig. 3.

CCRNs connectivity changed against age in Group [4–20 years old]. The colored CCRN regions show significant correlations (R${}^{\mathrm{2}}$ $\ge$0.6, p 0.001, cluster size $>$20 voxels) between the connectivity indices with age in Group [4–20 years old]. The color schema is the same as that in Fig. 2.

Fig. 4.

Examples of the scatter plot of the connectivity index correlated with age in Group [4–20 years old].

Unlike the results of Group [4–20 years old], Group [21–59 years old] only showed a minimal number of connection indicators that showed a significant negative correlation with age. That is, within CON, the TI of lAI-dACC/rAI, and between CON and DMN, the GC of PCC$\to$rAI were found attenuated against age (0.6 $\le$R${}^{\mathrm{2}}$$\le$ 0.8, p 0.001, FDR corrected, see Fig. 5 and Table 2 for more details). The correlation is weaker than that of most indicators detected in Group [4–20 years old]. No significant positive correlations were detected between the connection indicators and age.

Fig. 5.

CCRNs connectivity changed against age in Group [21–59 years old]. (A) CCRNs regions show significant correlations (R${}^{\mathrm{2}}$ $\ge$0.6, p 0.001, cluster size $>$20 voxels) between the connectivity indices with age in Group [21–59 years old]. (B) demonstrates examples of the scatter plot of the connectivity index correlated with age. The color schema is the same as that in Fig. 3.

The results of Group [60–85 years old] (Fig. 6, Table 3) showed an inversed pattern compared to that of Group [4–20 years old]. Most of the connection indicators detected (R${}^{\mathrm{2}}$$\ge$0.6, p 0.001, FDR corrected) showed attenuation against age, distributed within CON and between CON and DMN. Within CON, the CC and TI of lAI/rAI-dACC were found significantly attenuated. Between CON and DMN, the CC of dACC/rAI-PCC, dACC/lAI-rIPL, lAI-VMPFC/ACC, and dACC-DMPFC/VMPFC/ACC, TI of dACC/AI-PCC, and GC of lAI$\to$lIPL were found significantly attenuated. Only the TI of lAI-rIPS was found to be enhanced.

Fig. 6.

CCRNs connectivity changed against age in Group [60–85 years old]. (A) CCRNs regions show significant correlations (R${}^{\mathrm{2}}$ $\ge$0.6, p 0.001, cluster size $>$20 voxels) between the connectivity indices with age in Group [60–85 years old]. (B) demonstrates examples of the scatter plot of the connectivity index correlated with age. The color schema is the same as that in Fig. 3.

Although the ranked mapping result showed some trend of change across different age groups, the ANOVA results did not reveal significant change between the age groups. Four points may explain why no significant results were disclosed using ANOVA: (1) The Group [4–20 years old] and Group [59–85 years old] belong to the two populations that were experiencing significant changes in their brain functional connectivity. For example, the Group [4–20 years old] covered participants from preschoolers, schoolchildren, preteens, teenagers, and even early young adults. It was not surprising that the variance of the connectivity indicators contributed by the developmental effect in this group was large. Similarly, the Group [59–85 years old] might also have a large variance contributed by aging. (2) The individual difference in functional connectivity also contributed much to the variance besides the age in all three groups. (3) The variance caused by the two effects mentioned above blurred the boundary of two neighboring groups, making detecting group differences difficult. (4) If the change was not monotonic, it was possible that few significant changes could be detected across groups. For example, if the connectivity increased with age in Group [4–20 years old] but decreased with age in Group [59–85 years old], the connectivity might show no change across the two groups. In summary, for the dataset analyzed in the current study, compared with correlation analysis, ANOVA might not capture the detailed alteration pattern across the groups. The null results of ANOVA indicated that dividing the participants into age groups and analyzing how CON-CCRN interaction changes with age within each Group is necessary.

Although the results of including all participants in the connectivity-age correlation analysis were insignificant, they also showed that it was necessary to divide the participants into different age stages, calculate the correlation within each stage, and summarize the changing pattern of a lifetime. A few points are worth discussing here:

(1) Including participants of all age groups, the number of samples increased, and the p-value might decrease. However, when examining the correlation with age, both R and R${}^2$ showed a significant decline, which showed that on a lifelong scale, the changes of most connectivity indicators were probably not linear, and each age stage has its change pattern, which demands the analysis strategy of examining within age groups. Our results also explained this point. For example, Group [21–59 years old] had the largest sample number (N = 119) and the largest age span, but the number of brain regions with changed connectivity significantly correlated with age was minimum, indicating that the connection between CON and CCRN regions is stable, and generally does not change significantly with age. When analyzing the whole cohort, participants aged 21–59 might be dominant in the sample size and age span. Therefore, it was likely that the details of network connection change of the other two age groups would be lost.

(2) Existing studies have shown that the trend of brain connection changes in the three age groups is very inconsistent. The meaning of using linear correlation to explore possible general change patterns at a lifetime scale from a data-driven perspective needs to be clarified.

(3) Although the results were insignificant regarding the criteria we applied in the current study, providing the so-called “null results” were still helpful. In the current study, it is reasonable to start with a correlation analysis across the whole cohort. The “negative results” of the whole cohort correlation analysis, like the “negative results” of ANOVA, provided insight into the reason for performing correlation analysis within each age group and enhanced the motivation of the analysis strategy applied in the current study.

Two of the main results of the current study were that the connections between the CON regions (dACC and bilateral AI) showed significant enhancement against age before adulthood but showed attenuation during aging in adulthood. The former suggested increased integration between nodes within CON, indicating a typical pattern of CON development at the age of 4–20 years, which is consistent with the previous study [54]. According to previous research [55], this integration may be related to the central regulatory role of CON among the large-scale functional networks [4, 12, 13, 56, 57]. The result also showed that the connectivity enhancement of dACC-lAI/rAI is more prominent than that of lAI-rAI, indicating a core position of dACC in the development of the hierarchy of CON node, which was in line with our previous findings that dACC is on a higher level than AI in control [52, 53]. On the other hand, in Group [4–20 years old], the enhancement of dACC-rAI connections seemed to be more significant than those of the other connections in CON. For example, the R${}^2$ of the correlation between rAI$\to$dACC and age reached 0.96, the largest R${}^2$ among the connections within CON. This result supports the view of the Triple Network Regulation Model that the rAI is central to the dynamic regulation of functional networks [13, 58]. It is worth noting that the ventral AI area showed significantly enhanced FC. In contrast, dorsal AI showed enhanced GC. This discrepant pattern was consistent with previous findings that ventral AI is more associated with emotional processing, and the dorsal side is more associated with higher-order cognitive processes [58, 59, 60]. These results suggest that the dorsal and ventral of AI undergo functional differentiation and integration into different subnetworks during development.

In both Groups with ages over 20, we observed attenuation of TI of dACC-lAI and rAI-lAI. The attenuation is mild and only showed in TI indicators when participants are younger than 60 (Group [21–59 years old]). For participants over 60 (Group [60–85 years old]), the attenuation became more significant in both CC and TI indicators. The mild attenuation from early adulthood to late middle age may reflect an early sign of the decline of CON. The more severe attenuation in participants over 60 reflects the consequence of the early decline. It is worth noting that the connection indicators between dACC and rAI, though prominently enhanced before 20 years old, did not show significant change during aging. This result indicated that the interaction between dACC and rAI might be a binding signature of developing the cognitive control functions in younger participants and maintaining them in the older participants under normal conditions. The varying pattern of dACC-rAI interaction implies that its abnormal alteration may be associated with a mental illness involving impaired cognitive control [61]. This needs further investigation beyond the normal population.

Task-positive networks such as CON, DAN, and CEN play a crucial role in the performance of cognitive tasks. Dosenbach et al. [56] suggested that the CON served as the regulatory core of the task-positive network, monitoring the CEN and DAN, which were involved in specific task trials. The present study found that the coupling between CON and DAN/CEN increased with age during development, which is consistent with some previous findings which demonstrated a gradual increase in FC tween CON and the DAN/CEN [23, 62] from childhood to adulthood. This enhancement has been proposed to be critical for cognitive control maturity [24]. In addition, GC indicators further revealed the stronger influence of rAI$\to$lDLPFC, dACC$\to$r/lDLPFC, and lIPS$\to$lAI but weakened the influence of rAI$\to$lIPS. Sridharan et al. [13] argued that the rAI, as the main site of stimulus reception/response, is first activated upon receipt of external stimuli and subsequently generates priming to activate the dACC, lAI to trigger hierarchical control. DLPFC and IPS are involved in top-down control of cognitive control signals at the sub-trial level of the task. The present study only found GC influence of lAI/dACC$\to$DLPFC but not the rAI$\to$DLPFC, and enhanced interaction between the lAI and IPS but a weakened interaction between the right AI and IPS during development, which may indirectly support the hypothesis that right AI regarding initiation.

Taken together, the enhanced coupling, total interdependence, and directional GC influence between CON and DAN/CEN before adulthood delineated a procedure of developmental maturation characterized by the integration of the so-called task-positive networks, including CON, DAN, and CEN.

In Group [21–59 years old], only a small number of connection indicators correlated with age. The scatter plot of most connection indices did not show a significant pattern of increasing or decreasing or nonlinear trends against age, despite this Group containing the largest number of participants (n = 119) in the current study. The result indicated that most of the CON-CCRNs connections do not change significantly or systematically against age and enter a stable stage after adolescence. The negative results of most connections may be a sign of matured CON-CCRNs integration and are in line with the proposed trend of cognitive control in humans [55, 63, 64]. On the other hand, disentangled CON-CCRNs interaction may impair emotional and cognitive functions. For example, previous research has reported altered connections of CON-CCRNs in psychiatric disorders such as depression [65], anxiety [66], and Alzheimer’s disease [67]. This stable stage could also be critical to maintaining normal functions in adulthood.

The enhancement during the development and the attenuation of CON-DMN connection indicators were the most noticeable results we observed in the current study (Fig. 7). DMN is well-known as a task-negative network constantly deactivated when an individual performs cognitively demanding tasks and shows more active signal change when the task is absent [38, 68, 69, 70, 71]. DMN is thought to be associated with mind wandering and free thinking and the source of internal interference to the external goal-directed task performance [38, 52]. The regulation of DMN during a task is important for many behaviors that demand cognitive control, such as attention [52, 72], working memory [73, 74], meditation [69], etc. The attenuated deactivation of DMN during cognitive tasks has been reported constantly in elders, patients with mild cognitive impairment (MCI) [75], patients with Alzheimer’s Disease [76], Children with attention deficit, and other developmental problems [77, 78]. Previous studies have suggested that effective CON-DMN communication is crucial for successfully regulating DMN activity during cognitive control [52, 79]. The antagonistic relationship (measured using anti-correlation [40]) between CON and DMN activity is a sign of effective communication both during the task and during the resting state.

In the current study, the CC indicator showed enhanced anti-correlation between CON-DMN in Group [4–20 years old], indicating increased communication from childhood to adulthood. Consistently, the TI of CON-DMN also showed enhancement against age in Group [4–20 years old], which confirmed the increased CON-DMN communication. This developing pattern of CON-DMN interaction may suggest a balance of internal and information processing and be related to the functional development and improvement of individual functional networks to cope with more complex cognitive tasks when children and adolescents grow up [76, 80].

On the other hand, the results in Group [60–85 years old] showed significant attenuation of the CON-DMN interaction in elder participants. For example, the anti-correlation measured in CC between many CON and DMN nodes was significantly weakened against age. Consistently, for dACC and rAI, two major nodes of CON, and PCC, a representative hub of DMN, the TI of dACC/rAI-PCC were also significantly weakened when participants’ age became older. The result showed a pattern contrary to that in the children and young group. This attenuation in CON-DMN interaction may underlie the neural basis of the decline in cognitive function during aging [81]. The attenuation was not observed in the connection indicators of CON-DAN/CEN, indicating that cognitive decline during aging is mainly related to the disrupted communication between CON and DMN instead of those between the task-positive CCRNs [81]. The attenuation was not observed in the connection indicators of CON-DAN/CEN, indicating that cognitive decline during aging is mainly related to the disrupted communication between CON and DMN instead of those between the task-positive CCRNs.

Although there are many studies on the change of functional connectivity with age, only a small number of studies [22] focused on how the connectivity of CON and CCRNs changes in a lifetime and support the proposed core control network theory [9]. However, the previous research has the following limitations. (1) It was only based on a single FC indicator using cross-correlation and could not provide directional information. (2) It performed correlation analysis across participants of all ages and only examined the relationship of linear, monotonic changes, ignoring that the change of CON-CCRN communication in a lifetime could be nonmonotonic.

Therefore, as an essential core cognitive control network, the changing communication pattern between CON and other CCRNs needs to be clarified. This study can be regarded as a follow-up study of the previous work examining the interaction within CON and between CON and other CCRNs, delineating how the interactions changes in different stages of life by dividing the samples into three different age groups and performing analysis in each group. From the perspective of analytical methods, the current study applied a comprehensive analysis based on multiple connection indicators, including classic FC and time-dependent TI and GC, where GC may also provide directional information on the interaction. The approaches applied in the current study help overcome the above-mentioned limitations of the previous study, providing further details of CON-CCRN communication change. The current analysis of the three different age stages clearly demonstrated that CON-CCRNs interaction change in a nonmonotonic fashion during the lifetime and the change pattern of each age stage has its own characteristics.

The current study divided the participant into three age groups. The selection of the separating point of the age may affect the result of the connectivity-age joint analysis. The definition of “young”, “middle-aged”, and “old” may vary in different contexts. The maturing of the brain also showed brain-regional differences and individual differences. Therefore, finding absolute separating points fitting for all conditions is nontrivial and unnecessary. The age of 20 is commonly used to separate mature individuals from immature in daily life, and age 60 is an empirical point that people start recognizing as “old”. We tried using the age of 18, 22, 55, or 65, but the results did not change much. Most of the patterns were preserved, indicating that the results were not sensitive to the selection of the empirical separating point.

In this study, we used global signal regression to alleviate the contamination of common noise embedded in the global signal. However, it is worth noting that global signal removal is a common practice in resting-state and task fMRI research, and there are both for and against it. Some people think it may confound the results if the global signal is strongly modulated by the task or other experimental factors [82, 83, 84]. Others suggest that it is conducive to observing localized neuronal effects that can be obscured by various global noise factors [85, 86, 87, 88, 89]. We proposed whether global signal regression should be performed depending on specific situations, such as different datasets and methodological considerations.

The current study has some limitations. First, it only conducted correlation analysis across participants of different ages. The results may be affected by individual differences across the participants. Because of the limitation of the dataset, each participant accounted for a single time point (the age when scanned), and it is hard to acquire continuous life-long longitudinal MRI data. We could not examine how the CON-CCRNs interaction develops and changes within an individual on a life-time level. Second, the current study focused on the linear relationship between connectivity and age; we adopted a classic and intuitive correlation method but could not explore the potential non-linear, complex relationship between connectivity and age. More sophisticated methods for detecting the non-linear trend may be needed in future studies better to delineate potential complex changes in the connection indicators. The current study only focused on the interaction between CON and several well-defined CCRNs. However, it is worth noting that using a data-driven approach across parcels defined by an atlas may provide insightful results that show overall brain connectivity development.