Longitudinal functional brain network reconfiguration in healthy aging

Abstract Healthy aging is associated with changes in cognitive performance and functional brain organization. In fact, cross‐sectional studies imply lower modularity and significant heterogeneity in modular architecture across older subjects. Here, we used a longitudinal dataset consisting of four occasions of resting‐state‐fMRI and cognitive testing (spanning 4 years) in 150 healthy older adults. We applied a graph‐theoretic analysis to investigate the time‐evolving modular structure of the whole‐brain network, by maximizing the multilayer modularity across four time points. Global flexibility, which reflects the tendency of brain nodes to switch between modules across time, was significantly higher in healthy elderly than in a temporal null model. Further, global flexibility, as well as network‐specific flexibility of the default mode, frontoparietal control, and somatomotor networks, were significantly associated with age at baseline. These results indicate that older age is related to higher variability in modular organization. The temporal metrics were not associated with simultaneous changes in processing speed or learning performance in the context of memory encoding. Finally, this approach provides global indices for longitudinal change across a given time span and it may contribute to uncovering patterns of modular variability in healthy and clinical aging populations.

Based on such studies, we can conclude that functional connectivity within higher order resting state networks (i.e., DMN) decreases while, between-network connectivity (i.e., DMN and executive control network) increases, resulting in less modular, "dedifferentiated" brains in older age (Chan et al., 2014;Ng, Lo, Lim, Chee, & Zhou, 2016).
In addition, these changes have been associated with less efficient cognitive functioning, highlighting the critical importance of this organizational characteristic of the human brain (Iordan et al., 2018;Gallen et al., 2017;Geerligs et al., 2015;Chan et al., 2014;Müller et al., 2016).
Although these studies have provided important insights into modular reorganization, most of them are, however, based on crosssectional data and, thus, fail to describe true longitudinal change in advanced aging.
Recent work in network neuroscience has given rise to multilayer network modeling suitable to track the evolution of modular reconfiguration throughout time (Braun et al., 2018;de Domenico, 2017). This framework is more powerful in contrast to more traditional approaches as it can provide an aggregate measure of change across multiple time points, such as network flexibility, which characterizes how frequently brain regions switch allegiance from one module to another over time (Bassett et al., 2011). Although still in its infancy, this approach has been used to investigate the dynamic (within-session) modular structure during resting state (Betzel, Satterthwaite, Gold, & Bassett, 2017;Gerraty et al., 2018;Shine, Koyejo, & Poldrack, 2016), task performance (Bassett et al., 2011;Bassett, Yang, Wymbs, & Grafton, 2015;Braun et al., 2015;Schlesinger et al., 2017;Telesford et al., 2017), development (Betzel et al., 2015), and brain disorders (Braun et al., 2016), showing interesting patterns of flexible network reconfiguration dependent on the "task at hand." Indeed, one study has used multilayer networks to investigate how the brain's modular organization evolves across the human lifespan, suggesting that some modules tend to be highly flexible and exhibit substantial reconfiguration throughout adulthood (Betzel et al., 2015).
For the present study, we used a longitudinal dataset comprised of four occasions of resting state functional brain imaging and cognitive testing in healthy older adults. We applied the multilayer modularity model to investigate how the brain's modular structure changes in healthy aging over a span of 4 years. Further, based on the previous studies connecting this functional property to cognitive functioning, we assessed how modular reconfiguration relates to changes in cognitive performance. More specifically, we calculated the functional flexibility, promiscuity, cohesion strength, and disjointedness to study the degree to which brain regions switch between communities as well as the pattern of it, and the recruitment coefficient to study the association of the obtained community structure to well-known resting state networks.
Moreover, we explored the relationship between functional flexibility and (a) processing speed and (b) learning performance in the context of memory encoding. These two cognitive domains were chosen as they have been found to be particularly vulnerable to aging effects (Salthouse, 2010;Schaie, 2005) and have been previously related to brain modularity in healthy elderly (Gallen et al., 2017;Geerligs et al., 2015;Iordan et al., 2018;Ng et al., 2016).
We hypothesized that the functional networks would show high flexibility across the time span of 4 years and with increasing age, suggesting longitudinal modular reconfiguration and instable modular architecture in older adults. Finally, we assumed that greater flexibility is associated with poorer performance in both cognitive domains.

| Participants
Longitudinal resting-state fMRI (rs-fMRI) data were taken from the Longitudinal Healthy Aging Brain Database Project (LHAB; Switzerland)-an ongoing project conducted at the University of Zurich (Zöllig et al., 2011). We used data from the first four measurement occasions (baseline, 1-year follow-up, 2-year follow-up, 4-year follow-up). The baseline dataset included 232 participants (age at baseline: M = 70.8, range = 64-87; females: 114). At each measurement occasion, participants completed an extensive battery of neuropsychological and psychometric cognitive tests and underwent brain imaging. The brain imaging session was conducted in close temporal proximity to the behavioral assessments (difference between behavioral and MRI assessments in days [M ± SD]: baseline: 2.2 ± 5.2, 1-year follow-up: 2.6 ± 5.2, 2-year follow-up: 4.3 ± 13.0, 4-year follow-up: 4.6 ± 9.3).
Inclusion criteria for study participation at baseline were age ≥64, right-handedness, fluent German language proficiency, a score of ≥26 on the Mini Mental State Examination (Folstein, Folstein, & McHugh, 1975), no self-reported neurological disease of the central nervous system and no contraindications to MRI. The study was approved by theethical committee of the canton of Zurich. Participation was voluntary and all participants gave written informed consent in accordance with the declaration of Helsinki.
For the present analysis, we only included participants with complete rs-fMRI data (four measurement occasions). This was necessary because four temporal windows enable a more reliable maximization of the multilayer modularity function. At 4-year follow-up, the dataset still comprised 74.6% of the baseline sample (n = 173), of which 86.7% (n = 150, age at baseline: M = 69.8, range = 64-83; females: 71) had complete data for rs-fMRI.
To estimate whether attrition was selective, we compared the full sample at baseline with participants that had rs-fMRI data from all four measurement occasions. The total selectivity was computed by standardizing the difference between the mean in the baseline sample and the sample with no missing data, on the SD of the baseline sample in the variable of interest (Lindenberger, Singer, & Baltes, 2002).
The size of the selectivity index was interpreted with reference to an effect size. As it can be seen in Supplementary Table S1, the total selectivity was negligible for all measures (i.e., none of the measures exceeded the cutoff of 0.20 for a weak effect according to Cohen (1988)), indicating that the participants with no missing rs-fMRI data did not significantly differ from the full baseline sample in terms of age at baseline, education, initial cognitive ability, or physical and mental health.

| Neuropsychological assessment
All participants completed a neuropsychological test battery assessing multiple cognitive domains at each measurement occasion. For the present analysis, data from the domains "processing speed" and "learning/memory encoding" were used. Individual scores were standardized to T scores (M = 50, SD = 10) with respect to baseline and averaged across subtests to calculate the domain-average composite scores.
Processing speed was assessed using four psychometric paper-  (Reitan & Wolfson, 2004) (the scores were reversed so that the higher scores equaled better performance); and (d) number of correct responses (within 2 min) on the LPS14, a subtest from the Leistungsprüfsystem 50+ (LPS), a German intelligence test developed to measure Thurstone's (1938) primary mental abilities (Horn, 1983).
Learning/Memory encoding was defined using the (a) number of correctly reproduced abstract designs at first five trials of the DCS figural memory test (Diagnosticum für Cerebralschädigung; Weidlich & Lamberti, 2001), and (b) total correct responses over five immediate free recall trials from the Verbal Learning and Memory Test (Helmstaedter & Durwen, 1990), a German equivalent of the Rey Auditory Verbal Learning Test.
Further preprocessing steps were performed using the CONN toolbox (Whitfield-Gabrieli & Nieto-Castanon, 2012). The nuisance regressors were defined according to the 36-parameter model : six motion parameters, signals estimated from CSF and WM, global signal, their derivatives, quadratic terms, and squares of derivatives were regressed out from functional data separately for each run. The rs-fMRI data was temporally bandpass filtered in the 0.01-0.1 Hz frequency range. We applied simultaneous filtering/nuisance regression, because it was shown to reduce correlation between time-series fluctuations and motion (Hallquist, Hwang, & Luna, 2013). Global signal regression was performed in accordance with previous studies on aging (Chan et al., 2014;Ng et al., 2016), as this has been shown to be effective in the reduction of the effects of physiological signals and head motion (Lydon-Staley, Ciric, Satterthwaite, & Bassett, 2019).

| Network definition
The regions of interest (ROIs), used to build the network, were selected from the atlas of Schaefer and colleagues , corresponding to 200 cortical regions (i.e., ROIs) classified into seven well-known resting state networks according to the The time-series were averaged over the voxels in each ROI and correlated between each pair of nodes using the Pearson's correlation analysis. Each correlation coefficient was then Fisher's r-to-z transformed. Due to ambiguity regarding the meaning of negative correlations in the context of global signal regression (Murphy & Fox, 2017), negative z values were excluded from the analysis in accordance to previous studies (Chan et al., 2014(Chan et al., , 2018. This resulted in subject-specific 200 × 200 correlation matrices with diagonal and negative values set to zero. Therefore, the main network analysis was performed on positive weighted networks.
In addition, supplementary analysis has been done on positive weighted networks containing only statistically significant connections at an FDR-adjusted significance level of p < .01. This has been done to ensure that our results are robust in both unthresholded and thresholded matrices based on the significance level of node connections.

| Multilayer modularity
Modular structure in a network indicates that nodes in a module are more interconnected with one another then with the rest of the network. Modules within complex networks are identified using community detection algorithms such as the maximization of the modularity quality function (i.e., Q; Newman & Girvan, 2004). These algorithms are applicable to traditional single layer networks (e.g., within-session static functional connectome) or multilayer networks in which nodes are connected across layers (Kivelä et al., 2014). These layers may correspond to different modalities (e.g., structural and functional connections between brain regions) or different temporal instants at which the network was observed (i.e., brain connectivity across a span of several years). So, in addition to being able to calculate single window modularity scores, the multilayer framework is interesting for longitudinal studies because it can produce a "global index" which reflects the changes in network organization across a defined period (i.e., 4 years). In the present study, we define multilayer networks where each layer is the functional connectivity matrix at a given time point, resulting in a four-layer temporal network for each of the subjects, separately.
Thus, to investigate the changes in network organization across time, the multilayer modularity was optimized (Blondel, Guillaume, Lambiotte, & Lefebvre, 2008;Mucha, Richardson, Macon, Porter, & Onnela, 2010) as follows: where 1 is the number of layers in the multilayer network, A ijl is the functional connectivity matrix, P ijl is the corresponding null model matrix (i.e., Newman-Girvan null model) defined as the where m is the average edge weight in the matrix, g il gives the community assignment of node i in layer l, and g jr gives the community assignment of node j in layer r. The γ is the structural resolution parameter which defines the weight of intralayer connections and thus the number of obtained modules, while ω is the temporal resolution parameter which sets the weight of the interlayer edges that link each node i to itself across layers. When the value of γ is small, the maximization of Q produces relatively large communities, while large values result in more communities with smaller number of nodes. Given that the value of ω defines the consistency of multilayer modules, large values relative to intralayer edges result in communities that are more similar to one another across layers. We set these parameters to frequently used default values of γ = 1, ω = 1 Telesford et al., 2017). Additional analyses with varying parameters around these default values can be found in the supplementary material. As the multilayer community detection algorithm is stochastic, the modu- The algorithm was used with the randomization option "moverandw" instead of the default "move" option, as this has been shown to mitigate some undesirable behavior for multilayer modularity with ordinal coupling (Bazzi et al., 2016).

| Multilayer metrics
In order to describe the temporal variability of community (i.e., modular) structure, we calculated the node flexibility score which represents the number of times that a node switches communities over time, normalized by the total possible number of switches (Bassett et al., 2011). We obtained the global flexibility scores by averaging over all brain regions included in the analysis.
To better understand the underlying mechanism of the network flexibility over time, we calculated additional three measures: node promiscuity (Papadopoulos, Puckett, Daniels, & Bassett, 2016), node cohesion strength, and node disjointedness (Telesford et al., 2017).
In order to investigate the tendency of brain regions to change allegiance between limited or multiple (i.e., temporal integration) communities, we calculated the node promiscuity which reflects the fraction of all networks in which the node participates at least once, across all network layers (Papadopoulos et al., 2016). The global promiscuity was defined as the average promiscuity over all nodes.
Node cohesion strength and node disjointedness quantify the node changes based on mutual versus independent changes, respectively (Telesford et al., 2017). and lower values imply infrequent switches of communities of a given node independently of other nodes. The calculated metrics were averaged over all brain regions in order to obtain the global indices of disjoint and cohesive changes of a brain network across time.
To understand, how the present dataset corresponds to known functional networks, we compared the community partitions obtained on this data to the seven predefined resting state networks (Schlesinger et al., 2017): FPCN, DMN, DAN, SVAN, LIMB, SM, and VIS.
First, the cooccurrence of brain regions was summarized with a module allegiance matrix, where the ijth element present the percentage of time points in which both region i and region j belong to the same community (Mattar, Cole, Thompson-Schill, & Bassett, 2015).
Thus, regions that are consistently grouped together in the same network, across time, have high allegiance values.
Then, using the module allegiance matrix, we computed the dynamic recruitment of a network, that estimated the probability that its regions cooccur in predefined networks with regions from the same network across time points (Mattar et al., 2015). The recruitment coefficient was averaged over all brain regions in order to obtain the global dynamic recruitment.
In addition, all of the calculated multilayer metrics were summarized across predefined networks (Schlesinger et al., 2017); to further explore the significance of the obtained results in the context of wellknown resting state networks.
As previously stated, the multilayer community detection algorithm is stochastic and therefore the obtained measures were averaged across 100 optimizations of the modularity quality function.
To provide better intuitive understanding of these metrics, Telesford et al. (2017)

| Null models
To determine the statistical significance of the temporal evolution of functional brain networks and to test against the null hypothesis that there is no smooth reconfiguration between consecutive time points, we compared the real functional networks to a temporal null model.
This null model was created by shuffling the time layers in the multilayer network uniformly at random across time (Chai, Mattar, Blank, Fedorenko, & Bassett, 2016;Sizemore & Bassett et al., 2017). Importantly, the temporal null model constructed in such a way, preserves connectivity within a network layer but eliminates the dependencies between layers over time.
Further, we computed a nodal null model to contrast the network recruitment coefficient obtained on our data against the null hypothesis that roles of the regions in the network are identical (i.e., no functional subnetworks). Hence, we ensure that the community structure we calculate is not random but instead captures the modular organization of functional connectivity. This null model, similar to the configuration model for static graphs, was constructed by randomly rewiring edges occurring at the same point (Sizemore & Bassett et al., 2018).
Therefore, we computed 50 null models (temporal and nodal) for each subject's multilayer functional matrix and then optimized the multilayer modularity quality function 100 times on each of these null model networks.
The real networks were compared to null networks using Welch's two-sample t test implemented in R (v. 3.5.2; The R Project for Statistical Computing; http://www.R-project.org/).
The significance level of p-values was adjusted to p < .05/n, where n represents the number of multilayer measures that were compared between the observed and null networks.
Further, we calculated the Cohen's d index of effect size using an R-based package lsr (v. 0.5), and interpreted it as follows: d ≥ 0.2 was considered a "small" effect size, d ≥ 0.5 represented a "medium" effect size and d ≥ 0.8 a "large" effect size (Cohen, 1988).
2.9 | Statistical analysis 2.9.1 | Brain modular reconfiguration and age at baseline To test the hypothesis that the brain modular reconfiguration is related to aging, we performed multiple linear regression analysis , in which the outcome was a particular multilayer measure and the predictor was age at baseline (grand-mean-centered variable). Gender (female = 1, male = 0) and education (on a scale from 1 to 3; 1 = high school with or without vocational education, 2 = higher education entrance qualification, business school or university of applied sciences, or 3 = university degree) were entered as nuisance covariates into the model, as previous studies have related these variables to brain's functional network organization (Chan et al., 2018). Further, we included motion as an additional nuisance covariate, defined as the average framewise displacement (FD) across four measurement occasions, as head motion has been shown to have an important impact on brain network topol- The effect sizes (i.e., partial η 2 ) were calculated using the lmSupport package (v. 2.9.13) in R (v. 3.5.2). The significance level of p-values was adjusted to p < .05/n, where n represent the number of models tested.
In addition, to investigate the relationship between the metrics, Pearson's correlation analysis was used to calculate the association between the modularity index (Q), global recruitment, and global flexibility. The analysis was performed using an R-based package called psycho (v. 0.4.0.). In addition, partial correlation analysis was conducted to test the relationship between these metrics while controlling for the effects of age at baseline using the ppcor (v. 1.1) package in R.

| Longitudinal change in cognitive performance
We performed linear mixed effects (LME) analysis (lme4 package The results were visualized using the ggplot (v. 2-3.0.0) package in R.

| Brain-cognition association
To investigate the "change-change" brain-cognition association, we The individual rate of change in cognitive performance, defined as the subject-specific slope of the regression line between time and the cognitive scores, was derived from the LME models described in the previous section (Ng et al., 2016(Ng et al., , 2018. The results were visualized using the ggplot (v. 2-3.0.0) package in R.

| Longitudinal reconfiguration of functional modules
We used the modularity maximization algorithm to define the temporal communities in the multilayer functional connectivity matrix (distributions of the number of modules are plotted in Supplementary Figure S1).
First, we calculated the dynamic recruitment coefficient, which quantifies the probability that nodes of a network are consistently assigned to the same module across different time layers. Then, we compared the observed global recruitment to that of a nodal null model, to ensure that the community structure we obtained is not random but instead captures the organization of well-known resting state networks. The module allegiance matrix in Figure 1a provides a summary representation of how brain regions and networks are dynamically engaged across four time points.
The observed global recruitment coefficient was significantly higher than in the nodal null network (t(298) = 71.5, p < .0001, d = 8.26) (Figure 1b). This result suggests that the regions tend to be recruited to their own networks and are less integrated with other networks across time, suggesting that the community structure we obtained is not random but represents meaningful modular organization. Self-recruitment scores for each of the networks can be found in Supplementary Figure S2.
Further, we calculated several temporal measures in order to characterize the patterns of modular change across time, which we then compared to a temporal null model to determine whether the obtained values were higher or lower than expected.
We investigated the longitudinal network reconfiguration by calculating the flexibility score, which indicates the number of times nodes switch their community assignment across temporal layers, and global promiscuity, which is defined as the fraction of all communities in which the node participates at least once, across all network layers.
The observed global flexibility (t(298) = 59.2, p < .0001, d = 6.84), and global promiscuity (t(298) = 69.6, p < .0001, d = 8.04), were significantly higher than in the temporal null model suggesting that the functional brain displayed more change than expected in a temporal null network (Figure 2a,b).
Further, to better understand the underlying mechanism of brain flexibility over time, we calculated another two measures: global cohesion strength and disjointedness. Finally, Supplementary Figure S4 shows that our results were robust in both unthresholded and thresholded matrices.

| Negative correlation between modularity and flexibility
Next, we wanted to assess the relationship between flexibility and the quality of modular decomposition. In order to do so, we used the modularity index Q, describing the quality of modular structure across p < .0001) remained significant after controlling for age at baseline.

| The role of age at baseline in modular reconfiguration
We were interested in assessing if multilayer measures were related to age at baseline in our sample of elderly subjects (Figure 4).
There were no significant age effects on global promiscuity (b = .0008, p = .277, partial η 2 = 0.0081), suggesting that there were no age differences in the number of networks the nodes switch between (Table 1).
We did not observe any significant association between the global cohesion strength (b = .0845, p = .306, partial η 2 = 0.0072) or disjointedness (b = .0001, p = .182, partial η 2 = 0.0123) and age at baseline, implying that the extent of mutual or independent changes in community structure (over the 4-year time interval), was not linked to age at baseline.
F I G U R E 1 (a) The module allegiance matrix represents the probability that two brain regions are part of the same community across the 4-year interval. The brain regions are ordered according to the predefined network they belong to. Higher values along the diagonal of the matrix suggest that networks from the Schaefer atlas tend to be recruited together in the same communities across the 4-year interval. (b) Comparison of global recruitment of real networks to a nodal null model. The center red lines represent the mean, and the light red bars and light blue bars represent 95% confidence interval and SD, respectively. This figure was generated using notBoxPlot (https://github.com/raacampbell/ notBoxPlot) These results suggest that although there are age differences in the extent of global flexibility, the mechanisms of this change were not significantly related to age and thus homogeneous across the sample.
Next, the global recruitment was significantly negatively associated with age at baseline (b = −.0035, p = .006, partial η 2 = 0.0501), suggesting a lower probability of a given region to be  grouped in the same community as other regions of its own community across time in older age (Table 1). This result points to higher instability in modular organization and network reconfiguration with advancing age.
Motion and education were not significantly related to any of the metrics.
Importantly, the associations between age at baseline and flexibility survived multiple comparison correction (p < .05, corrected for five models). As with the global measures, we ran multiple regression models to assess the relationship between network-specific flexibility and age at baseline (Table 2). Interestingly, the flexibility of DMN (b = .0076, p < .001, partial η 2 = 0.0824); FPCN (b = .0043, p = .014, partial η 2 = 0.0407); and SM (b = .0038, p = .047, partial η 2 = 0.0269) was higher with older age, suggesting that older subjects have higher flexibility in these networks ( Figure 6).

Results
Linear models with remaining networks did not show any significant association between age at baseline and network flexibility (Supplementary Table S2).
Of all of the tested models, only the association between age and DMN flexibility survived multiple comparison correction (p < .05, corrected for seven networks).
Finally, we wanted to explore the mechanisms of flexibility for the networks that had significant effects of age at baseline; therefore, we computed linear regression models with network-specific promiscuity and recruitment for the DMN, FPCN, and SM networks.
The results showed higher DMN promiscuity (b = .0041, p = .002, partial η 2 = 0.0627) with older age (Figure 7a, Supplementary   Table S3). This can be interpreted as a higher tendency of brain regions belonging to the DMN to segregate from this network and connect to all other networks across the whole brain.
Not surprisingly, the self-recruitment of the DMN (b = −.0063, p = .006, partial η 2 = 0.0501) and SM (b = −.0085, p < .001, partial η2 = 0.0851) was lower with higher age (Figure 7b,c, Supplementary Table S3), suggesting that the network flexibility across the 4-year interval was related to lower probability of brain regions from these networks to be categorized into same communities as other regions from the corresponding networks.
In Supplementary Figure S6, in which the module allegiance matrix was calculated for the youngest 10% and for the oldest 10% of the sample, we can see that modules, and especially the DMN, tend to be more consistently recruited across time points in younger in comparison to older subjects.
All statistically significant effects survived multiple comparison corrections (p < .05 corrected for six tests).

| The role of modular reconfiguration in cognitive performance
First, we performed LME analysis to investigate longitudinal change in cognitive functioning (Supplementary Figure S7). We found a statistically significant decline in processing speed (b = −.5276, p < .001) and learning/memory encoding (b = −.5001, p = .003). Older age was associated with lower performance in both domains (Table 3). Finally, there was a significant interaction between age at baseline and time (b = −.0512, p = .029), suggesting that older participants had a more significant decline in processing speed across the 4-year time interval.
Also, participants with better education had higher scores in learning/memory encoding (b = 2.1896, p = .002). In addition, there was a significant retest effect on processing speed (b = 1.7722, p < .001), but not on learning/memory encoding (b = .3933, p = .470).
F I G U R E 5 Visualization of summary statistics for networkspecific flexibility across the 4-year interval. The center red lines represent the mean, and the light red bars and light blue bars represent 95% confidence interval and SD, respectively. This figure was generated using notBoxPlot (https://github.com/raacampbell/ notBoxPlot) T A B L E 2 Association between network-specific flexibility and age at baseline in multiple regression models. Statistically significant effects (p < .05) appear in bold F I G U R E 6 Network-specific flexibility in association with age at baseline. The results indicate greater modular flexibility of the frontoparietal control, default mode, and somatomotor networks across time in older participants All statistically significant effects (except the interaction effect on processing speed) survived multiple comparison corrections (p < .05 corrected for two domains).
Next, we conducted Pearson's correlation analysis to assess the association between flexibility and change in cognitive performance.
We did not find a significant relationship between change in F I G U R E 7 Network-specific promiscuity and recruitment in association with age at baseline. The results indicate greater modular promiscuity of the default mode and lower self-recruitment of the default mode and somatomotor networks across time in older participants Finally, the brain-cognition association was also tested for network-specific flexibility, however, there were no significant relationships found between any of the networks (i.e., FPCN, DMN, DAN, SVAN, LIMB, SM, VIS) and change in processing speed or learning/ memory encoding (see Supplementary material for further details).

| DISCUSSION
In the present study, we investigated the temporal change of the brain's modular structure in healthy aging. We applied the multilayer model, which resulted in several measures summarizing the characteristics of longitudinal network reconfiguration in healthy elderly.
We showed significantly higher variability and a greater range of modular "dynamics" over the course of 4 years in older subjects' functional networks in comparison to temporal null networks. Further, flexibility was significantly associated with age at baseline, with older participants having higher global and network-specific flexibility, particularly evident in the FPCN, DMN, and SVAN.
We also observed a decrease in the global network recruitment with older age, indicating community structure reorganization in which some brain regions are inconsistently assigned to their modules across different time points. This was most evident for the DMN which had lower time-dependent self-recruitment in older participants.
Over the same 4-year time interval, we observed a significant decrease in processing speed and learning/memory encoding (a finding which is well in line with previous cognitive aging studies (Ng et al., 2016;Salthouse, 2010;Staffaroni et al., 2018). However, the decline in cognitive performance was not related to the multilayer brain dynamics, implying an absence of simultaneous (i.e., changechange) relations between changes in functional brain network organization and cognition measures.
The multilayer modularity approach provides several advantages in comparison to more traditional methods for community detection.
Similar to "single-layer" modularity, it does not require a predefined set of networks, it is completely data-driven (de Domenico, 2017) and applicable to an individual level (Shine et al., 2016). However, in contrast to "single-layer" modularity, it partitions all temporal layers simultaneously, maintaining a consistent set of modules across all layers, thus ensuring the same definition of networks across all time points (Mucha et al., 2010).
To our knowledge, this is the first study with the application of multilayer community detection on longitudinally acquired data in healthy elderly. Nonetheless, our findings are in line with previous studies, suggesting unstable network architecture and substantial functional reconfiguration with aging.
We also found a highly significant negative correlation between global flexibility and recruitment, which implied that higher network flexibility is related to a more random nature of brain dynamics in which functional modules are not persistently recruited across time.
Modular brain networks exhibit a fine balance of dense withinnetwork connections and sparse connections between regions in networks with different processing roles (Meunier et al., 2010).
In our recent study, including the present sample (but also encompassing participants with missing data at some time points), we explored changes in the functional segregation of resting state networks across time . We showed a decrease over a 4-year interval in the functional segregation of associative networks, including the default mode, FPCN and SVAN networks. Thus, it is possible that a loss of within-network integrity might have resulted in nodes grouping in fewer and larger modules, and at the same time losing functional segregation and seeing more nodal movement between modules across the 4-year time interval (Schlesinger et al., 2017).
These findings are in line with other research that suggested increased modular variability or heterogeneity within higher order cortices in healthy elderly, indicating that the brain reconfigures during the aging process and varying cognitive demands (Peraza et al., 2018;Schlesinger et al., 2017). Furthermore, another study showed less similarity of network partitions in older healthy subjects, both as a group and across time, compared to younger participants, which implied reduced stability of network organization with aging (Iordan et al., 2018).
In the present study, we also calculated node promiscuity, cohesion strength, and disjointedness, and compared these metrics to temporal null models in order to better understand the underlying mechanism of brain flexibility over time. The observed cohesion strength was significantly higher than in the null network, while the disjointedness was significantly lower, implying that with aging subgroups of nodes cohesively reorganize into new modules instead of individual brain nodes switching communities independently of other nodes. Nevertheless, this reconfiguration pattern is related to less specific modules with more fluid connectivity between them and could possibly indicate compensatory reconfiguration of functional networks due to declining cognitive performance or impaired recruitment mechanisms (Sala-Llonch, Bartrés-Faz, & Junqué, 2015). Further, cohesion strength was not significantly associated with age at baseline, which means that the mechanism of longitudinal change in functional configuration was more uniform across the sample.
Although global promiscuity was significantly higher than in the temporal null model it was not significantly associated with age, further reinforcing the notion of a more uniform mechanism of global flexibility across the age span found in our sample.
Importantly, we found a significant relationship between flexibility and age, more specifically, older subjects tended to have higher flexibility in several predefined resting state networks, such as the FPCN, DMN, and SM networks. All networks exhibited significant levels of longitudinal flexibility, but these networks seem to present more variable change across the age range present in our sample.
However, it should be noted that only the DMN flexibility-age association survived multiple comparisons correction, so the overall results should be interpreted with caution.
Moreover, given that there were some network-specific effects on functional flexibility, we wanted to explore if the mechanism of change-promiscuity and recruitment-was diverse across age groups, despite the absence of age effects on the global average of these metrics.
Our findings pointed to higher promiscuity coupled with lower self-recruitment of the DMN in older age, which can be interpreted as brain regions of this network being inconsistently recruited across time and less segregated from regions belonging to other communities. Lower self-recruitment was also found for the SM network, pointing to a lower functional specialization of these networks across the 4-year interval in healthy elderly.
The most consistent finding across studies on aging is that older adults have lower functional integrity in the DMN, compared to younger adults (Chong et al., 2019;Damoiseaux, 2017). Moreover, this network has been shown to be highly vulnerable to aging-associated diseases such as Alzheimer's disease or cerebrovascular disease (Chong et al., 2017;Crossley et al., 2014;Kim et al., 2016).
Importantly, the relationship between aging-related functional changes within these networks and changes in cognitive performance is not yet fully understood. In our study, we observed significantly lower processing speed and learning/memory encoding in older participants with a decline over the 4-year time interval. However, we did not find any significant association between this reduction in cognitive performance and the multilayer measures, contrary to our expectation.
Although there are no previous studies investigating longitudinal network flexibility in the context of cognitive performance, some crosssectional results did indicate relevant associations between modular properties and cognition in older adults (Gallen et al., 2017;Geerligs et al., 2015;Iordan et al., 2018). This lack of simultaneous relation may have been driven by some compensatory mechanism, where healthy aging individuals delay the effects of functional reconfiguration for a certain time and thus maintain cognitive ability (Reuter-Lorenz & Park, 2014). Our study did not focus on identifying the causal role of temporal variability in cognitive functioning, but future studies should investigate lagged coupled changes in order to test if functional changes precede cognitive changes in healthy aging.

| Methodological considerations and limitations
Although the methods applied here provide an interesting framework for investigating the longitudinal functional reconfiguration in the elderly, this approach has several methodological considerations that should be taken into account.
First, it is well-known that the choice of nodes can significantly influence the calculation of network properties (Fornito, Zalesky, & Bullmore, 2016). We defined our nodes according to a functional atlas comprising seven resting state networks Yeo et al., 2011), also commonly used in studies on healthy aging. Second, the number of networks in the original atlas approximated the number of modules obtained in our analysis (even though there was high interindividual variability), thus allowing more straightforward comparability between our data-driven modules and predefined networks.
Apart from that, we only included participants with complete data (all time points) in order to maximize the overall number of temporal layers. Selective exclusion of subjects with incomplete data might introduce some bias into this type of analysis, although there was no selective attrition in our sample. As a consequence, the future research should explore different strategies for handling missing data in the context of multilayer modularity.
Finally, although multilayer metrics present a convenient method for summarizing change information across time points, this "summarization" might also obscure different rates of change existing between specific time windows, which could be more easily explored within a different statistical framework.

| CONCLUSION
This study, for the first time, illustrates substantial functional network reconfiguration in healthy aging across a 4-year time interval. In particular, the whole brain network flexibility, which reflects the tendency of brain nodes to switch between modules was significantly higher in healthy elderly than in a temporal null model and with increasing age. The modular temporal variability was not related to simultaneous changes in cognitive performance, however, further studies should include more cognitive domains or investigate lagged changes to better understand the temporal implication of the multilayer modular reconfiguration. Finally, this approach provides simple intuitive indices for overall longitudinal changes across a desired time span and it can be useful for uncovering patterns of modular variability in healthy and clinical aging populations.

ACKNOWLEDGMENTS
The current analysis incorporates data from the Longitudinal Healthy Aging Brain (LHAB) database project carried out at the University of Zurich (UZH). The following researchers at the UZH were involved in the design, set-up, maintenance, and support of the LHAB database: Anne Eschen, Lutz Jäncke, Franz Liem, Mike Martin, Susan Mérillat, Christina Röcke, and Jacqueline Zöllig.

CONFLICT OF INTERESTS
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
Susan Mérillat, Franziskus Liem, and Lutz Jäncke contributed to the design, set-up, maintenance, and support of the Longitudinal Healthy Aging Brain (LHAB) database. Brigitta Malagurski performed the data analysis and wrote the first draft of the manuscript. All authors discussed the results, contributed to manuscript revision, read, and approved the submitted version.

DATA AVAILABILITY STATEMENT
The data for this manuscript are not publicly available. Since data collection was started in 2011, when public data sharing and open science were not yet widely discussed, the used consent does not allow for the public sharing of the data. We are currently working on a solution for this matter. At the moment, data can only be accessed via collaborations with the URPP "Dynamics of Healthy Aging".