Individual functional parcellation revealed compensation of dynamic limbic network organization in healthy ageing

Abstract Using group‐level functional parcellations and constant‐length sliding window analysis, dynamic functional connectivity studies have revealed network‐specific impairment and compensation in healthy ageing. However, functional parcellation and dynamic time windows vary across individuals; individual‐level ageing‐related brain dynamics are uncertain. Here, we performed individual parcellation and individual‐length sliding window clustering to characterize ageing‐related dynamic network changes. Healthy participants (n = 637, 18–88 years) from the Cambridge Centre for Ageing and Neuroscience dataset were included. An individual seven‐network parcellation, varied from group‐level parcellation, was mapped for each participant. For each network, strong and weak cognitive brain states were revealed by individual‐length sliding window clustering and canonical correlation analysis. The results showed negative linear correlations between age and change ratios of sizes in the default mode, frontoparietal, and salience networks and a positive linear correlation between age and change ratios of size in the limbic network (LN). With increasing age, the occurrence and dwell time of strong states showed inverted U‐shaped patterns or a linear decreasing pattern in most networks but showed a linear increasing pattern in the LN. Overall, this study reveals a compensative increase in emotional networks (i.e., the LN) and a decline in cognitive and primary sensory networks in healthy ageing. These findings may provide insights into network‐specific and individual‐level targeting during neuromodulation in ageing and ageing‐related diseases.


| INTRODUCTION
Healthy ageing is typically characterized by subtle declines and slowing in general cognitive abilities, which may further progress into pathological impairment and even dementia. The development of noninvasive neuroimaging methods, such as various magnetic resonance imaging (MRI) sequences, has made it possible to detect structural and functional brain markers in ageing, such as regional cortical thickness (Frangou et al., 2022), white matter integrity (Miller et al., 2016), and functional network organization (Power et al., 2011). Structural brain changes have been extensively studied in the context of healthy ageing and related neurodegenerative disorders (Yan et al., 2018) but are assumed to occur later than functional changes (Jack et al., 2010;Zonneveld et al., 2019). Functional studies, especially studies on functional network organization, could help motivate investigations of or interventions for ageing-related neurodegenerative disorders, such as Alzheimer's disease (AD) (B. Wang et al., 2017;Xu et al., 2021). For example, functional networks can serve as individual intervention targets in healthy older adults and AD patients when using noninvasive neuromodulation to enhance cognitive functions (Nilakantan et al., 2019).
Resting-state functional MRI (rs-fMRI) measures intrinsic lowfrequency (<0.1 Hz) blood oxygenation level-dependent (BOLD) signal activities (Cordes et al., 2001) and has revealed large-scale functional networks, including those linked to high-order cognitive and emo-  (Yeo et al., 2011). These networks are organized to support different brain cognitive functions by segregating and integrating within and between networks H. Y. Zhang et al., 2021) and have revealed neural dedifferentiation in healthy ageing, which indicates decreased independence, decreased segregation of functional networks and inability to specify relevant neural circuits to mediate specialized functional processes (Chan et al., 2014;King et al., 2018). Static functional connectivity (FC) analysis measures statistical temporal correlations of mean BOLD time series signals between distinct brain regions across the entire scanning time and has shown that the ageing brain undergoes complex functional reorganization and compensation (H. Zhang et al., 2017;Zonneveld et al., 2019). Ageing-related reorganization is characterized by weaker within-network connectivity and controversial between-network findings, that is, greater between-network connectivity (Geerligs et al., 2015;Spreng & Turner, 2019;Zonneveld et al., 2019) or weaker between-network connectivity (Varangis et al., 2019;H. Y. Zhang et al., 2021).
Using a "dynamic" analysis, in which time-varying, dynamic FC is measured within a series of overlapping "sliding windows" in the BOLD timeseries data (Preti et al., 2017), studies have found that older adults spend more time in a baseline, weak connectivity state (Tian et al., 2018), indicating loss of dynamics and the inability to adapt to environmental variations (Garrett et al., 2021). In addition, older adults may spend less time in a series of different states, as revealed by the following discrepant findings. The brain state may be characterized by antagonistic activity of the DMN and AN (K. Y. Chen et al., 2022), by high connectivity within the MN and the cognitive control network (Tian et al., 2018), or by many positive connections among subnetworks (Xia et al., 2019). Overall, static between-network connectivity and dynamic functional network state findings are controversial, and we suspect that these discrepancies are due to the use of different analytic methods, especially node or region of interest (ROI) spatial definitions, when constructing individual functional networks.
Like many other complex networks, nodes and edges are two basic elements in brain functional networks. Generally, nodes are typically defined by ROIs in a predefined, group-average and nonoverlapping resting-state network (RSN) parcellation. Edges, also named FCs, are typically defined by Pearson or partial correlation coefficients of averaged BOLD time series in ROIs. However, of particular importance, different RSN parcellations may assign the same voxel or vertex to different networks, especially voxels (or vertices) in subcortical networks (Doucet et al., 2019). In addition, most group-level parcellations are derived from young adult (age <40 years) data (Gordon et al., 2016;Power et al., 2011;Yeo et al., 2011) but are not suitable for participants across the lifespan. Inspired by this question, one study mapped individual functional parcellation and revealed location reconfiguration of functional regions in ageing adults (Geerligs et al., 2017), indicating inaccurate calculation of node signals and FC using young adult brain parcellations. Recent studies have focused on age-appropriate brain parcellation, including an age-appropriate functional parcellation derived from older adults, ages ranging from 55 to 95 years (Doucet et al., 2021), and five cohort-specific parcellations (age range: 20-34, 35-49, 50-64, 65-79, and 80-93 years). However, no matter what the group-level parcellation is, averaging individual BOLD signal data based on group-level network parcellation (Braga & Buckner, 2017) may underestimate certain participantspecific properties and details of functional network architecture, such as cross-individual variations in the shape, size, and position of functional networks (Bijsterbosch et al., 2018). Therefore, it is necessary to map individual-level parcellation and further analyse ageing-related changes based on it.
Recent methodological breakthroughs, in which individual-level parcellation has been mapped using iterative clustering analysis on group-level parcellations (M. L. Li et al., 2019;D. H. Wang et al., 2015), offer the opportunity to characterize dynamic functional network organization changes with age at the individual level. In addition, researchers have noticed the drawback and limitation of fixedlength sliding window analysis and proposed data-driven segmentation of sliding windows (Choe et al., 2017), such as the hidden Markov model (G. M. Zhang et al., 2020), the dynamic conditional correlation model (Lindquist et al., 2014), and activation-informed temporal segmentation (Duda et al., 2021). Considering differences in interindividual parcellation, internetwork organization and interindividual sliding window activities with increasing age, we used individual participant parcellation and individual network dynamic analysis on individual sliding window lengths to advance the understanding of dynamic network organization. In addition, we correlated recurring dynamic brain states with multiple cognitive behavioural measures to assign different states to corresponding cognitive levels. According to the compensation hypothesis (Cabeza et al., 2018;Reuter-Lorenz & Park, 2014), we expected to find inverted "U-shaped" correlations (impairment) between high cognitive level brain states and age in most networks and "U-shaped" correlations (compensation) in specific networks.
The included participants were not diagnosed with diseases that would impact brain functions, such as dementia, AD, Parkinson's disease, multiple sclerosis, stroke, or epilepsy (more exclusion criteria are shown in Shafto et al. (2014)). Two participants were excluded from the analyses due to excessive head motion (details below), resulting in 637 participants (18-88 years old, mean ± SD age = 54.25 ± 18.45 years; 311 males and 326 females) included in the final sample. All participants gave written informed consent, and the Cambridgeshire 2 Research Ethics Committee approved the study.
Participants' exclusion criteria of Cam-CAN and more details are shown in a previous publication (Shafto et al., 2014). In addition, this study included 13 cognitive variables derived from 8 outside-MRI cognitive tasks, including the fluid intelligence task, the hotel task, the picture-picture priming task, the proverb comprehension task, the visual short-term memory task, the choice motor coefficient of variation task, the face recognition task, and the emotion expression recognition task. These cognitive tasks are used to evaluate five cognitive domains, including executive functions, language functions, memory function, motor function, and emotional processing. A brief description of each variable is summarized in Table 1, and full descriptions are given in a previous publication (Shafto et al., 2014).

| Data preprocessing
Rs-fMRI data were preprocessed using FSL (https://fsl.fmrib.ox.ac.uk/ fsl/fslwiki/) with the following steps: discarding the first 10 volumes, correcting head motion by MCFLIRT, slice timing correction, extracting nonbrain tissues with BET, spatial smoothing with full width at half maximum = 6 mm, normalizing intensity, high-pass temporal filtering (cut-off frequency = 0.01 Hz) and registering the rs-fMRI to highresolution T1-weighted structural images. To exclude the head motion effects, head movements exceeding 2 mm or 2 in any direction were discarded (n = 2). The final sample included 637 participants. In addition, the covariate head motion was calculated as the average of the root mean squared realign parameters at all 251 time points using MCFLIRT.
T1-weighted structural images were processed using the FreeSurfer version 6.0.0 software package. Total cortical grey matter volume was calculated as a covariate in statistical analyses (details below). The structural and functional images were aligned using boundary-based registration. Rs-fMRI data were aligned to a spherical coordinate system by sampling from the cortical ribbon in a single interpolation. The rs-fMRI data of each individual were first registered to the FreeSurfer surface template, which consisted of 40,962 vertices in each hemisphere. The smoothed data were then downsampled to a mesh of 2562 vertices in each hemisphere using the mri_surf2surf function in the FreeSurfer software package. individual-level cortical network parcellations were segmented into discrete patches using a clustering algorithm. Third, the patches were matched to 116 cortical ROIs extracted from 18 group-level networks, resulting in individual ROI parcellations. It is important to note that if a patch did not overlap with any cortical ROI and was not near any ROI, the patch would be labelled "unrecognized." Thus, the number of individual cortical ROIs was less than or equal to 116. Considering the intersection of individual ROI parcellation, only 27 homologous ROIs were retained for 100% of participants, and the number of ROIs was too small to analyse individual functional variations. Therefore, we kept ROIs that were defined for 90% of participants. Finally, 88 ROIs were kept for subsequent analysis.

| Individual parcellation and functional network size calculation
To better summarize between-participants functional variations in network size and compare them with previous findings, we further grouped the 18 networks into seven well-studied functional networks, including the AN, DMN, FPN, MN, LN, SN, and VN. Given any of seven networks, the change ratio of the network size (S r ) was quantified as follows: where m indicates the number of vertices in the given network, l k indicates the assigned label of vertex k within individual-level ROI parcellation, and l j indicates the assigned label of vertex j within group-level ROI parcellation. The assigned label was 1 or 0. Finally, the change ratio of the network size was calculated for each network for each participant.

| Clustering-based dynamic functional state construction
For each participant, 251 volumes (time points), with 88 ROIs in each volume, were included in the following analysis ( Figure 1). At each time point, the average of the preprocessed BOLD signal across all vertices in a given ROI was extracted. Considering the different resting-state FC patterns of different networks, we separately applied a sliding window and clustering analysis to each of the seven networks. Given one network and one participant, dynamic FC was calculated using a sliding window method, with varied window length (ranging from 18 TR to 28 TR) and fixed sliding step length (1 TR). The optimum window length (L) at time point t was defined as follows: x tþ a À 1, roi ð Þ ! where x indicates whether the given ROI signal at time point (t +a À 1) was local extrema. If the ROI signal was larger (or lower) than the signal at the former and the latter time points, the ROI signal at the current time point was a local extremum, and x was equal to 1; otherwise, x was equal to 0, a indicates a different window length (range from 18 to 28), and n indicates the number of ROIs within a given network.  (Aggarwal et al., 2001). To improve the stability of clustering, we iterated k-means clustering 100 times.
Since clustering analysis was performed on each network and clustering centres were matched for all networks, k was set to 2 to match F I G U R E 1 Workflow of analysis. (a) Individual dynamic functional connectivity construction. Group-level network templates including seven networks were extracted from a previous publication (Yeo et al., 2011). Individual preprocessed functional magnetic resonance imaging (fMRI) time courses were iteratively clustered based on the group network template to generate individual network parcellation. Then, the network parcellation was separated into the region of interest (ROI) parcellation. ROI signals were extracted to calculate dynamic functional connectivity (FC) using individual sliding windows. (b) Groupwise clustering and statistical analysis. Dynamic FC within each network of all subjects was combined into a feature matrix. Clustering analysis was performed on these seven network feature matrices, resulting in clustering centres for each participant. Using canonical correlation analysis, clustering centres were grouped into the strong state and weak state. Finally, groupwise changes with age were calculated and summarized to determine ageing impairment and compensation. AN, attentional network; DMN, default mode network; FPN, frontoparietal network; LN, limbic network; MN, motor-sensory network; SN, salience network; VN, visual network the clustering results with low deviation. In addition, two clusters were widely reported in previous fMRI studies (Kim et al., 2017;H. Li et al., 2022;Zheng et al., 2022). Finally, clustering algorithms labelled FC matrices into one of the two clusters, and the median FC in the same cluster was calculated as a cluster centroid.
Here, the cluster centroids were defined as brain states.

| Statistical analysis
To exclude the influences of covariates, we first constructed a linear The Kolmogorov-Smirnov test was used to assess normality, and amplitude differences in the two clustering states were tested by the Wilcoxon rank sum test. Since multiple variates were included in amplitudes of each state (7 networks) and cognitive task (13 variables), correlations between amplitudes and cognitive tests were modelled using canonical correlation analysis (CCA) (Tibon et al., 2021), a powerful method to simultaneously examine linear relationships between multiple amplitudes and cognitive tasks derived from each participant. Then, the statistical comparison of correlation coefficients of CCA was performed using the "cocor" package (Diedenhofen & Musch, 2015).
The significance level was set to p < .05, and a false discovery rate correction was applied to correct multiple correlations. The above statistical analyses were performed using MATLAB R2014a, R 4.

| Dynamic functional state changes with ageing
To investigate changes in the two states with ageing, we examined how the temporal properties, including the mean dwell time, and FPN (AIC diff = 1.23) showed significant negative linear correlations with age ( Figure 4a). SN (AIC diff = 8.15) showed significant quadratic correlations with age, and "two-lines" tests were significant (Figure 4b).
For State 2 (high amplitude and cognitive correlations), dwell time of AN (AIC diff = 14.52) and MN (AIC diff = 12.02) showed significant quadratic correlations with age, and "two-lines" tests were significant.
FPN (AIC diff = 7.13) and VN (AIC diff = 11.40) also showed significant quadratic correlations with age, but "two-lines" tests were not significant ( Figure 5) Since the sum of occurrence of two states was constant (= 224 windows), their relationships with age were opposite, and the signifi- quadratic correlation with age: adjusted r 2 = .082, p < .001), and LN (AIC diff = 1.50; quadratic correlation with age: adjusted r 2 = .030, p < .001). Of these, only the correlations of LN first decreased and then increased (Figure 7). GEE showed that the occurrence of the DMN was higher than that of any of the other six networks (all p < .001), and the occurrence of the FPN was lower than that of any of the other six networks (all p < .05). In addition, the occurrence of SN was larger than that of MN (p = .015) and AN (p = .009).
Significant inverted U-shaped patterns were found in AN (upward slope: p = .0213; downward slope: p < .001) and MN (upward slope: showed an inverted U-shaped correlation with age, suggesting the difficulties of maintaining the same brain states in the healthy ageing process. In terms of strong state analysis, DMN occurred and dwelled less with increasing age; FPN and VN began to occur and dwell less in middle age (50-53 years); AN and MN occurred and dwelled more before middle age (48-52 years) followed by a decline. These findings suggest that in the healthy ageing process, the resting brain may allocate insufficient time to the strong state, characterized by a high correlation with cognitive function, which may further lead to inefficient cognitive resource allocation. In addition, we have shown that the LN strong state dwells more with increasing age and occurs more in middle-aged adults (older than 47 years), indicating the compensation of the LN strong state for the impaired process of most functional network states.
In this study, individual cortical parcellation was mapped using established group-level parcellation (Yeo et al., 2011)  is relatively maintained. These findings support the theory that latedeveloping networks (i.e., high-order networks) are sensitive to the ageing process (Douaud et al., 2014;Westlye et al., 2010). Notably, we found a significant positive linear association between age and the change ratio of LN size, suggesting that LN size may consistently expand across the adult lifespan. It has been suggested that changes in the LN are different from other high-order cognitive networks during development and the ageing process. One previous early adult (18-45 years) study showed that the cortical thickness of primary sensory and high-order networks is negatively correlated with age, except for the LN (Bajaj et al., 2017). In addition, one previous study mapped group parcellation of children and compared it to adult parcellation (Yeo's parcellation (Yeo et al., 2011)).
They found that the parcellation of primary sensory networks is similar, but the DMN regions in adult parcellations are partially assigned to the LN in children (Tooley et al., 2022). These findings support our  (Oosterwijk et al., 2012) and plays a significant role in emotional and reward-related processing (Cao et al., 2021;Rolls, 2019). Compared to other higher-order and primary sensory networks, the LN has shown the largest functional network pattern changes during the first 6 years of life (H. T. Chen et al., 2021). An independent component analysis compared FC differences between young (18-33 years) and old (58-85 years) adults and found differences in cognitive networks but not in emotion networks (Nashiro et al., 2017), suggesting cognitive function decline and emotional function maintenance with age. In line with these results, we found that the LN strong state occurred and dwelled more with increasing age, indicating the compensative function of the LN. This is supported by a local FC analysis on the same dataset, which also found decreased local FC within the VN, SMN, and DMN but increased local FC within the basal ganglia network with increasing age (Wen et al., 2020). In addition, one study analysed baseline, 1-, 2-, 3-, and 4-year follow-up rs-fMRI data from old adults (64-87 years) and reported decreased segregation of the DMN, FPN, and SN and an increase in the LN (Malagurski et al., 2020). Overall, compared to the decline in other cognitive networks, such as the DMN, FPN, and SN, the LN may maintain and even compensatively increase with increasing age (Nashiro et al., 2012;Scheibe & Carstensen, 2010).
Several limitations should be considered when interpreting these findings. First, this study reports ageing-related changes in a group of participants, not within individuals, as individual studies require longitudinal follow-up. Second, this rs-fMRI dataset collected 261 volumes (8 min and 40 s) for each participant, and more volumes would result in more precise individual parcellation (D. H. Wang et al., 2015).
However, the relatively less precise parcellation is tolerable compared to group parcellation in ageing studies (Geerligs et al., 2017). Third, we clustered all dynamic FCs into two brain states, as done in previous studies (Kim et al., 2017). However, some studies clustered more brain states, such as three (Tian et al., 2018) and five (Xia et al., 2019).
Perhaps each participant has an individual number of recurrent brain states, but this would make group comparisons more difficult. Thus, more research is needed to describe individual brain states.