Functional Connectivity Dynamics of the Resting State across the Human Adult Lifespan

The aging brain undergoes alterations of Structural Connectivity (SC), but also transformations of the flexible interactions between brain regions during cognition, described by Functional Connectivity (FC). FC in the resting state (rs) displays a variety of age-related changes, but most studies do not address the spontaneous reorganization of FC across time. Here, we show by means of functional MRI (fMRI) human brain imaging that aging also profoundly impacts on the dynamics of this rs FC. Analyzing time-dependent correlations between human rs fMRI blood oxygen level dependent (BOLD) time series, we reveal a switching, markedly slowing down with age, between epochs of meta-stable FC and transients of fast network reconfiguration. Furthermore, we identify communities of functional links whose temporal fluctuations become increasingly anti-correlated in elderly subjects. This effect is stronger when the performance is reduced in cognitive screening tasks. Such remodeling of Functional Connectivity Dynamics (FCD) discloses qualitatively novel effects of aging that cannot be captured by variations of SC or static FC. Our new metrics manifest thus a strong biomarking potential, by achieving a quantitative parameterization of FCD variations across the human adult lifespan. Significance statement Aging alters the brain’s anatomical connectivity and inter-regional activity correlations (Functional Connectivity, or FC). Here we show that aging also affects the temporal remodeling of FC networks (Functional Connectivity Dynamics, or FCD) during the resting state. By analyzing human fMRI data, we find that the speed of FCD decreases with age. Furthermore FCD properties predict differences in cognitive performance between same-age subjects. This slowing down of FCD intriguingly parallels the reduction in information processing speed often invoked as a cause of age-related cognitive decline. We believe that FCD analyses will lead to better biomarkers of healthy and pathologic aging, since they avoid misinterpreting temporal as inter-subject variability.


Introduction
As we age, our brain undergoes structural and functional changes. Tractographic studies revealed a tendency toward cortical 'disconnection', as manifested by an agerelated decline in white matter integrity (1,2). These decreases in structural connectivity (SC) are paralleled by disruptions of resting state (rs) functional connectivity (FC). Interregional BOLD correlations are reduced within rs networks (RSNs) such as the Default Mode or the Dorsal Attentional networks (3,4), implicated in attention, memory and executive control functions, which decline in cognitive aging (5,6). However, changes in rs FC are not just reflecting SC changes. For instance, FC between RSNs tends to increase with age (7,8). The correlation between SC and FC may either increase or decrease depending on the considered region and these variations in SC-FC coupling predict age better than SC or FC alone (9). Here, we show that aging not only affects rs FC but also its spontaneous reconfiguration over time, i.e. Functional Connectivity Dynamics (FCD).
Recent studies emphasized the structured temporal variability of rs FC (10-15). If rs FC is dynamic, a wealth of information may be lost by averaging over long imaging sessions. Averaged temporal variability might be mistaken as inter-subject variability. Furthermore, temporal FC variability may carry an inherent meaning, by manifesting ongoing cognition at rest (16) with a direct impact on cognitive performance (17, 18) or by reflecting the sampling of a repertoire of dynamical states (19). From a biomarking perspective, pathological conditions, such as Alzheimer's dementia or schizophrenia (20, 21), or differences in general attributes like developmental stage or gender (22, 23), may alter FCD more than they affect time-averaged FC.
Using a brain imaging data set including N = 85 healthy human adult subjects (between 18 and 80 yrs), we quantified age-related variations in the rate of reconfiguration of FC networks -'speed' of FCD' -and describe the spatio-temporal coordination patterns (meta-connectivity) that constrain the time-variability of FC links. We found that FCD at the whole cortex level markedly slows down with age and that its freedom of reconfiguration is increasingly 'frustrated'. FCD-based markers were also predictive of individual differences in cognitive performance between subjects of similar age, as probed by standard clinical assessments (24). This capacity of tracking cognitive performance hints at the relevance of FCD for actual mental operations, where, speculatively, the slowing of FCD may relate to the reduction in processing speed postulated by theories of cognitive aging (25, 26).

Results
Aging slows down Functional Connectivity Dynamics. A widespread manner to extract FC from rs fMRI is to parcellate the brain into N macroscopic regions -see Table S1 and reference (27) -and to compute pairwise correlations between the regionaveraged time series of neural activity a i , based on the entire rs imaging session (tens of minutes). The result of this procedure is an N-times-N FC matrix (Fig. S1A). In order to estimate the time dependency of FC, we adopted a common sliding window approachfollowed e.g. by Allen et al. (12) -, repeating the FC construction separately for each time-window (of a fixed duration τ). By this we generated a stream of FC(t) matrices (Fig. S1B) and studied the similarity between different time-resolved networks. To do so, we introduced a metric of similarity between FC matrices and evaluated Functional Connectivity Dynamics (FCD) matrices (19). The FCD matrix entries FCD ab provide the normalized correlation CC[FC(t a ), FC(t b )] between any two FC(t) networks observed at times t a and t b . Fig. 1A shows FCD matrices for subjects of different ages (see Fig. S2 for different τs). FC(t) matrices were stable during epochs lasting several time windows, as visualized in Fig. 1A by square-shaped red blocks along the main FCD matrix diagonal, denoting high similarity between FC(t) networks within the epoch. Such epochs of transient FC stability -or 'FCD knots'often terminated abruptly, and were intertwined with transients of instability, visualized in Fig. 1A by bluish stripes in the FCD matrix, denoting strong dissimilarity from previously visited FC(t) networks. During these transients -or 'FCD leaps' -FC(t) quickly morphed before stabilizing again into the next FCD knot. Note that the alternation between FCD knots and leaps does not imply the existence of crisply separated FC states (see Discussion). The typical duration of FCD knots depended on subject age. As shown by Fig. 1A, FCD in older subjects seemed to slow down with respect to younger subjects, as revealed by longer lasting FCD knots.
To provide a quantitative description of the slowing down of FCD, we provided a statistical characterization of the rate of stochastic change in time of rs FC networks. We described the dynamics of FC as a stochastic exploration of the space of possible FC configurations. We sampled FC(t) at discrete times t 0 , t 1 = t 0 +τ, … t k = t 0 +kτ, … and evaluated the lengths of the steps traveled in FC space between each of these 'stroboscopic' observations. We measured step lengths as the correlation distance between two consecutive FC(t) networks, i.e. d τ [t i ] = 1 -CC[FC(t i ), FC(t j + τ)]. In this way, FCD step lengths are large when two consecutive FC(t) observations are poorly correlated (as within a FCD leap) and close to zero when two consecutive FC(t) observations are highly correlated (as within a FCD knot). These measured d τ can also be interpreted as describing a speed of FCD reconfiguration, since consecutive FC(t) observations were always separated by the same amount of time τ. We sampled the statistical distributions of d τ for different ranges of τ and different subjects (see Materials and Methods). Distributions of d τ for two representative subjects of different ages are shown in Fig. 2B. For all considered τ-s and subjects (see also group averages in Fig. S3A and different τ-s in Fig. S2B), these distributions displayed a peak at a value d typ , which we call the typical FCD speed. As shown in Fig. 1C, d typ significantly decreased with age (-0.74 < R 50 = -0.58 < 0.44, 99% bootstrap c.i.), thus confirming our intuition.
For the analyses underlying Figs. 1B-C, we used a range of time scales (pooled windows sizes from 20 to 8 s), which is unusual and relatively very fast when dealing with BOLD signals. This allowed however extracting longer d τ sequences, yielding an improved estimation of d τ distributions at the single subject level (see Discussion). The same trend was found as well for larger τ-s, but age correlations were weaker (Fig. S3B).
The sampled d τ distributions had a second property. Their left tail had a slow decay (red distributions in Fig. 1B) and was significantly fatter than chance expectations (grey distributions in Fig. 1B), evaluated under the null hypothesis of no sequential correlations between consecutive steps in the FC(t) stream (removed by random shuffling of the order within the stream). Such over-representation of short step lengths d τ was consistently observed across all subjects and τ-s (Fig. S3C) and is due to the fact that FC(t) networks transiently stabilized within FCD knots. Their speed of reconfiguration was thus reduced to only increase again when a FCD knot 'melted' into a FCD leap. The graphical cartoons in Fig. 1D visualize our interpretation of the phenomenon. The FC sequences in experimental data, alternating between short steps within FCD knots and longer steps during FCD leaps (Fig. 1D, top), bear analogies with anomalous stochastic processes showing persistence -i.e., long-lasting sequential auto-correlations of increments -such as Lévy flights (28, 29). The resulting FCD trajectories contrasted with the null hypothesis of vanishing sequential correlations between step lengths (Fig. 1D, bottom) in which the trajectories would spread homogeneously over the same volume in FC space.
We rigorously confirmed that FCD instantiate an anomalous stochastic process by performing a Detrended Fluctuation Analysis (DFA) of the sequences of d τ , supplemented by a (Bayesian) model-selection step (30). This robust procedure quantifies the strength of auto-correlations in a sequence by detecting a power-law scalingdescribed by a scaling exponent α DFA -in the divergence of a quantity ! , probing the strength of the fluctuations of d τ at different scales of observation k (see Supplementary  Methods). A value of α DFA = 0.5 corresponds to a Gaussian white noise process, in which the expected fluctuation after N uncorrelated steps grows as √N. In contrast, larger values 0.5 < α DFA ≤ 1 correspond to so-called fractional Gaussian noise which is anomalously persistent (α DFA = 1.0 corresponding to a Lévy walk). Model comparison allows for discarding subjects for which a genuine power-law scaling is not present (30). Fig. 1E shows DFA log-log plots for two representative subjects (based on a FCD analysis with a short τ = 12 s). The slopes of the straight lines fitted to log ! provide an estimate of α DFA . For all window-sizes up to τ < 29 s, at least 70% of the tested subjects showed a genuine power-law scaling of ! , according to (Bayesian) model selection (for τ = 12 s, N = 46 out of 49 subjects). For all of them, α DFA was significantly larger than 0.5 (p < 0.05, bootstrap c.i.). The median α DFA was closer to 1.0 for elderly than for younger subjects (Fig. 1F). We also found a significant positive correlation between α DFA and age at the single subject level (Fig. S4A), robustly up to window sizes of τ ~ 29 s. (Figs. S4B-C). For longer τ-s -and correspondingly shorter d τ sequences -, correlations were no longer significant.
A slower FCD does not necessarily mean that FC topology becomes more stable. On the contrary, when considering the actual modular structure of FC(t) networks, we observed that its flexibility over time increased with age ( Fig. S5), as also reported by Schlesinger et al. (31). Smaller FCD changes in elderly subjects are thus more effective in mixing FC communities (see Supporting Results).
Aging frustrates Functional Connectivity Dynamics. We converted the FC(t) stream of Fig. S1B into a collection of time series describing the time-dependency of individual FC pairwise couplings (Fig. S1C). By extending the FC matrix construction from regional nodes to inter-regional links, we built a N(N-1)-times-N(N-1) matrix of correlations between the time-dependent strengths of N(N-1) FC links (N 2 pairs of regions, minus the self-loops). The resulting Meta-Connectivity (MC) matrix described inter-link covariance, similar to FC describing inter-node covariance. MC matrices are particularly suitable for inter-group comparisons. In fact, they can be averaged over multiple subjects within a homogeneous group (expected to share common inter-link correlations), while this cannot be done for FCD matrices (portraying specific realizations of FC fluctuations, different for each subject in a group).
MC captures higher-order correlations between triplets or quadruplets of brain regions, beyond pairwise FC. Indeed, estimating the strength of a meta-link MC ij,kl between two functional links FC ij and FC kl requires monitoring the coordinated activity of minimum three brain region (for MC trimers, when the two links share a common root region, i = k) or, generally, four (for MC tetramers, when the two links do not share any vertex). Thus MC analysis tracks the static high-order correlation structures (32, 33) shaping the dynamics of second-order pairwise FC. Fig. 2A shows average MC matrices for different age groups, based (here and in the following) on sliding FCD windows of τ = 29 s. Both these MC matrices exhibited a modular structure, characterized by communities of temporally co-varying FC links. For each individual MC module, Fig. 2B displays a cortical map of the associated trimer weights MC(i), i.e. of the sum of the weights of trimer meta-links MC ij,il belonging to a given MC module and being rooted in a region i. These maps (averaged over the two hemispheres, given the absence of large asymmetries) allow to visualize meta-hub regions with elevated MC(i), around which the different modules are organized. While a conventional FC hub region is a network node whose activity fluctuations highly correlate with those of adjacent nodes, a meta-hub region is characterized by a high degree of correlation between the fluctuations of its incident FC links. MC modules and the associated meta-hub regions characteristically distribute across the whole rostrocaudal range (Fig. 2B).
These modules do not overlap with any standard RSN (34). We would like to remark that focusing on trimer meta-connectivity exaggeratedly emphasizes the localized nature of the MC modules. By considering all meta-links together -the majority of which were tetramers -, we found that every MC module 'touched' nearly every region through at least one meta-link (Fig. S6A). Hence, functional meta-networks formed pervasive higher-order correlation scaffolds that coordinate the temporal fluctuations of brain-wide FC networks.
A striking effect of aging on our MC results is the surfacing of negative meta-links, apparent from the darker blue hue of the elderly group MC matrix ( Fig. 2A, right). The results in Fig. 2C revealed that, although the median strength of individual MC weights decreased for elderly subjects (p < 0.001, Kruskal-Wallis), the bulk of their distributions largely overlapped (blue 95% c.i. in Fig. 2C). Nevertheless, the elderly group's distribution developed a tail of negative MC outlier entries (dashed red ranges in Fig. 2C). Negative meta-links were found within and between all MC modules ( Fig. S6B-C). At the single subject level, we found that neg-MC -defined as the average strength of all negative MC entries -, anti-correlated significantly with age (Fig. 2D). Thus, the freedom for links of a given MC module to follow its own 'fluctuation choreography' was frustrated by the rise of inter-module conflicts.
Functional Connectivity Dynamics predict cognitive performance. We were also interested in whether FCD alterations correlate with the age-related decline of cognitive and motor performance, which varies greatly between subjects of a same age group, and particularly among elderly subjects (35). We first related FCD to cognitive performance in a subset of old subjects between 60 and 70 yrs tested via the Montreal Cognitive Assessment (MoCA). MoCA is a standard cognitive screening tool for Mild Cognitive Impairment (24), which recapitulates into a single global score the partial performances achieved in subtests for visuo-spatial abilities, language, orientation, attention, working memory and other executive functions.
The global MoCa score significantly correlated with the typical FCD speed d typ (the faster the FCD, the better the score, Fig. 3A) and with neg-MC (the more negative the MC, the worse the score, Fig. 3B). We then considered performance in a simple unimanual visuo-motor coordination task, adapted from (36). In this task, performance was measured, following (37), as the average frequency-locking -quantified by a coefficient 0 < Φ < 1 -between the rotations of two circles presented on a screen, a visual cue and a second, subject-generated, cyclic force (squeezing an air-filled rubber ball at certain frequency, see Supporting Materials and Methods). Once again, the correlations of Φ with d typ (Fig. 3C) and with neg-MC (Fig. 3D) were significant (see Table 1).
Importantly, correlations between FCD-related metrics (d typ or neg-MC) and performance indicators (MoCA or Φ) continued to be significant even when regressing out the common declining trend with age (see partial correlations, CC(·, ·| Age), in Table  1). Thus, neg-MFC variations at a common age were indicative of different cognitive performance levels.
Among the partial MoCA scores, only the working-memory subtask score correlated significantly with neg-MC or d typ (e.g. CC(neg-MC, MoCA-wm | Age) = 0.65, p < 0.01, bootstrap). DFA could not be performed on the subset of subjects undergoing cognitive and motor performance assessment due to the shorter length of their rs fMRI scans.

Discussion
The intrinsic time-variability of rs FC may carry useful, generally disregarded information about ongoing neural dynamical processes and their age-related alterations (10-23). Here we characterized the coordinated evolution of rs FC, revealing correlations of FCD and MC with subject's age and also with inter-subject variability in task performance.
Many studies in 'chronnectomics' have attempted isolating FC states that stereotypically reoccur across time and subjects (13, 14). However, the retrieved states and their number may be affected by the clustering algorithm used and by the availability of data. Furthermore, it may be difficult to statistically prove that FC is switching between states that are really distinct (38, 39). Here we bypassed FC state extraction, by focusing on sequential correlations in continuous and stochastic FC(t) streams. We found that FCD trajectories 'explore' the space of FC configurations at a temporally inhomogeneous speed (Fig. 1D). Such anomalous random process may occur even within a single FC state, in such a way that the distance between FCD knots could become compatible with a scenario of overall FC stationarity (38).
The alternation between FC knots and leaps may reflect the underlying complex dynamics of cortical circuits, with system's trajectories alternatively exploring faster and slower manifolds in a complex system's phase space (19, 36). From a functional perspective, walks of the Levy type have been associated to optimal search (37), e.g., for food in an ecological environment. It is tempting to consider anomalous stochastic FCD an efficient 'foraging' for cognitive resources by searching for FC patterns adapted to the current information sharing and transfer demands (42).
Theories of cognitive aging have advanced that a cause for declining performance would be the insufficient access to cognitive resources due to a reduced speed of processing (25, 26). Cognitive aging has been associated to deficits in disengaging from active brain functional states (43, 44). In fact, aging affects FCD by reducing the reconfiguration speed of FC networks. Slowing down of cognition hence parallels slowing of FCD. Having said that, we do believe that more experiments should be designed to probe speed of processing or task switching to move beyond mere conjectures. Speculatively, certain modifications of FCD through aging may also be interpreted as compensatory mechanisms (45), e.g., the fact that elderly subjects' FCD leads to a more flexible FC modular structure, despite being slower (Fig. S5).
Aging has a widespread impact on the strength of SC and FC (1-4, 7) and on the temporal variance σ(FC) of individual FC links (44). As presented in the Supporting Results ( Fig. S7) age alterations of SC, FC, σ(FC), and MC are quite independent and yield non redundant age information. As a result, mixing different categories of predictors can lead to superior precision and sensitivity in age prediction and age-group discrimination (Fig. S8). The combined tracking of SC-, FC-and MC-based features may thus lead to superior biomarkers of healthy and pathological aging.
The age dependency of FCD may be the expression of artifacts, such as modifications of the neurovascular coupling itself (47) or other physiological processes of non-neural origin (48). Furthermore, stronger age-related effects on BOLD FCD are observed for short temporal windows (cf. analogous results in (31)), which are 'borderline' for the time resolution of fMRI. Future studies may simultaneously look at FCD based on signals that reflect neural activations more directly and have a better time resolution, such as EEG (49). Our findings of correlations between FCD and cognitive performance (Fig. 3, Table  1) add nevertheless to independent studies (17, 18) hinting at an actual link between rs BOLD FCD at fast temporal resolution and task-relevant neural information processing.
Mean-field whole-brain computational models (50) provide a welcome avenue for assessing the non-artifactual nature of FCD -virtual brains do not have blood -, while allowing for reverse-engineering its dynamic underpinnings. Generic whole-brain models are already able to qualitatively reproduce FCD (19), but future simulations might be fitted to individual subjects via automated pipelines (51). Models embedding SC typical of different age classes may reproduce the slowing down of FCD as an emergent byproduct of SC 'disconnection' itself. Or, more likely, they may show that this disconnection must be compensated to account for observations by a drift of the global 'dynamic working point' of operation of cortical networks, which could be possibly induced by altered neuromodulation (52).

Material and Methods
See Supporting Material and Methods for full detail.
Experimental data. Healthy adult subjects between 18 and 80 yrs were voluntarily recruited at Charité University Medicine Berlin to undergo rs fMRI (N = 49 subjects, 20 min; plus additional N = 36 subjects, twice 5 min; TR = 1.94 s) and DSI-based SC reconstruction, through pipelines as in (49). A Desikan partition (25) including N = 69 regions was adopted (Table S1). In addition we evaluated cognitive performance for a subset of subjects via a MoCA cognitive assessment (N = 21) (34) or via a unimanual visuo-motor task (N = 36) (32), evaluating for each subject the average performance over ten task repetitions. FCD analysis. According to a common approach (12,13), inter-regional rs BOLD correlation matrices FC(t) were estimated by sliding a temporal window of fixed duration τ. We calculated the similarity between any two FC(t 1 ) and FC(t 2 ) matrices observed at different times ( Fig. 2A) via inter-matrix Pearson correlation, following (19). Histograms of FCD step lengths were approximated by measuring the correlation distances d between FC networks in consecutive non overlapping windows and then -for the sake of a better sampling -by pooling together, d τ observations for multiple nearby window sizes τ, as justified by the similarity between their single-τ distributions. Window sizes between 280 and 8 s were pooled into scale-groups containing an equal overall number of observations (e.g., τ = 20 to 8 s for Fig. 1B-C or τ = 60 to 20 s and τ = 240 to 60 s for Fig. S3B). Histogram confidence intervals were evaluated according to standard binomial proportion statistics (Agresti-Coull interval). Histograms for the null hypothesis (absence of sequential correlations) were obtained by repeating the same construction above, but for random re-orderings of the empirical FC(t) sequences (500 permutations per window size). Standard DFA combined with a non-parametric model selection step was performed over sequences of d τ for different choices of window size, following (30) and using the associated toolbox.
MC analysis. Out of the stream of FC(t) matrices, we extracted N(N-1) time series of pairwise FC couplings. A normalized covariance matrix, or MC, between these time series was constructed. In all MC analyses presented in this study, we chose τ = 15 TR (~29 s) with a window step of one TR. MC matrices as well as neg-MC -defined as the average over all the negative MC matrix entries -were extracted for each subject. MC matrices were then averaged over age groups. MC communities were based on the young adult group-averaged MC matrix (18 to 25 yrs). Although we extracted communities of co-varying links, MC matrices were ordinary adjacency matrices of graphs having as nodes FC links. We could thus use conventional community-detection algorithms, such as the Louvain method.    Figure S1). (A) FCD matrices for subjects of increasing age (window size τ = 60 s, see Fig. S2 for other τ-s). Blocks of large inter-network correlation indicate epochs in which FC is stable (FCD knots), separated by transients of faster FC reconfiguration (FCD leaps). Visual inspection suggested that FCD knots lasted longer with increasing age. This is confirmed by computing the rates of FC reconfiguration d τ , or FCD speed. (B) Distributions of FCD speed, shown here for two representative subjects (log-log scale, pooled window sizes 8 s ≤ τ < 20 s) displayed a peak at a value d typ (typical FCD speed) and a fat left tail, reflecting an increased probability with respect to chance level to observe short FCD steps (95% confidence intervals are shaded: red, empirical; gray, chance level). (C) The FCD speed d typ decreased with age (see Fig. S3 for larger τ-s). The FC space was seemingly explored through an anomalous random process in which short steps were followed by short steps with large probability (sequential correlations), leading to clustered trajectories (panel D, top). This contrasts with a standard random process, visiting precisely the same FC configurations but without long-   . Partial correlations in which the common declining trend with age was regressed out were also significant (see Table 1).

Supplementary Results
The slowing down of FCD is paralleled by an enhanced temporal 'liquidity' of FC modules. FCD metrics quantify global variations of FC topology across time but do not capture whether these modifications induce or not a reorganization of graph-theoretical properties of FC(t) networks. It might be that local changes in the network sum up into a relatively large step length in terms of correlation distance, but leave global network properties such as modularity and community structure intact.
We studied how the modular structure of FC(t) networks varied across time, as an effect of ongoing rs FCD. We first separately partitioned each of the FC(t) networks into non-overlapping graph modules, using a standard Louvain algorithm. Since the modules obtained were not identical across the different time-windows, we had to relabel them and group them according to the similarity between time-slices -through a clustering of clusters procedure (see Supporting Materials and Methods) -in order to maximize the window-to-window overlap between matched modules. This allowed us for tracking when brain regions were changing its modular class, i.e. transiting toward a different community than the one they belonged to at a previous time.
A dynamic community analysis of the FC(t) stream is shown in Fig. S5. In Fig. S5A (top) the time-dependent modular membership of different brain regions is color-coded, for a representative subject. In the illustrated example, the sampled FC(t) networks could contain up to five different modules, and the number of detected modules fluctuated in time. At every time step we quantified the average 'liquid fraction' Λ(t) within a module, i.e. the expected fraction of regions within each module transiting to a different module at the next time step (see Supporting Methods). We found that the average number of detected modules decreased with age (Fig. S5B) and that the average 'liquid fraction' within a module increased with age (Fig. S5C).
The first result is in apparent contrast to Schlesinger et al. (31), since they find an increasing number of modules with age. This contrast may be due to the different clustering approach adopted in the two studies. While we performed clustering of nodes within each time slice (clustering in 'space') and subsequently re-clustered these clusters across time-slices in a second step (clustering in 'time'), Schlesinger et al. used a singlestep multi-slice clustering approach (clustering in 'space-time'). If modules strongly change in time, then the 'space-time' clustering will find a larger number of possible modules. On the contrary, a blurred inter-module distinction will yield a smaller number of 'space-only' modules. By construction, our two-steps -'space then time' -approach will never retrieve more modules than the largest number of 'space-only' modules that can be observed within individual time-slices (see Supporting Methods).
Both our study and (31) agree however in the important fact that the 'liquidity' of FC(t) modules was enhanced with increasing age. It appears that ongoing rs FCD has a stronger disrupting effect of the modular structure of FC(t), despite its slowing down. In other words, although FCD-induced changes are smaller, they go better 'on target', possibly facilitating inter-module information sharing, a phenomenon that could be speculatively interpreted as a novel type of compensatory mechanism in aging (45).

Relative age-predictive value of MC features.
We were primarily interested in the effects of aging on the temporal reconfiguration of FC networks. Therefore, we focused on correlations between age and FCD features. Many other features correlate with age, such as SC (cf. Fig. S1D), time-averaged FC (cf. Fig. S1A), and even the variance of individual FC links σ(FC) (cf. Fig. S1E). We may ask whether FCD analysis really provides additional information about the aging process, or the age-related information carried by FCD is redundant with the one already obtained by more conventional feature categories.
To identify possible redundancies, we studied all possible pairs of brain regions and for each pair of regions i and j and considered: the strength SC ij of the structural link between these pairs; the time-averaged strength FC ij of the inter-regional activity correlation; the variance over time σ(FC ij ) of this inter-regional activity correlation; and, finally, the integrated strength MC(ij) of all the meta-links involving the time-dependent functional link FC ij (t). We determined correlations with subject age for all these features for every pair of regions. In Fig. S7A we show scatter plots of these different correlations in which every dot corresponds to a distinct pair of brain regions. Plotting agecorrelations with one type of features vs. age-correlations with a different type of features helps to visualize dependencies between them (i.e. redundancy), whenever present. However, the scatter clouds shown in Fig. S7A were quite unstructured, indicating that it is difficult to anticipate whether, e.g., the strength of structural or functional links SC ij and FC ij also correlates with age when the integrated metaconnectivity strength MC(ij) does. All possible scenarios were present. For instance, there were MC-SC inter-dependencies: when MC(ij) was anti-correlated significantly with age, sometimes SC ij was anti-correlated, sometimes it was correlated (clouds of green dots in Fig. S7A) and some other times it did not correlate at all (clouds of red dots).
In Fig. S7B, we performed a similar analysis but considered node-level rather than link-level features: integrated SC strengths SC(i), given by the sum of all structural links SC ij incident to region i; integrated FC strengths FC(i), given by the sum of all functional links FC ij incident to region i; and meta-connectivity strengths MC(i), given by the sum of the strengths of all the trimer meta-links sharing region i as a root (cf. Fig. 2B). The corresponding scatter clouds resulting from the measured age-correlations were once again fairly unstructured.
A lack of low-dimensional inter-dependency between different types of features suggests that they provide complementary and non-redundant age-related information. Therefore, synergies between features may be useful for age prediction. To explore this eventuality, we performed age prediction through a simple linear regression scheme, based on different possible mix of link features, SC ij , FC ij , σ(FC ij ) or MC(ij) (See Supporting Materials and Methods). The cross-validated prediction errors for different predictor mix are shown in Fig. S8A. We first considered a linear age prediction from a single link feature at the time (one SC ij , or one FC ij , etc.). For every predictor category, the single link (ij) achieving the best cross-validated age-prediction performance is reported in Fig. S8B and the associated prediction performance in Fig. S8A. The median single-predictor performances achieved by FC ij and MC(ij) predictors were both better than the ones achieved by SC ij and σ(FC ij ) predictors (p < 0.05, Kruskal-Wallis with multiple comparisons correction). However FC ij and MC(ij) predictors were not significantly different from SC ij and σ(FC ij ) predictors.
Subsequently, we considered quadruplets of predictors of a same homogeneous type (i.e. four SC ij 's, four FC ij 's, etc.). Again, we identified the best quadruplet among all possible ones for each predicting feature category (see Fig. S8C). In this case, quadruplets of MC(ij) performed better than any other homogeneous-type quadruplet (p < 0.01, Kruskal-Wallis with multiple comparisons correction). Finally, we considered all possible quadruplets mixing the four types of predictors (i.e. composed by one SC ij , one FC ij , one σ(FC ij ) and one MC(ij) features), the best among which is reported in Fig. S8D. In our sample, this best quadruplet performed better than the best MC(ij)-only quadruplet (p < 0.05, Kruskal-Wallis with multiple comparisons correction).
To avoid a bias towards specific 'best' predictors (which may be highly specific to our own dataset, despite the care we took in cross-validating the performances), we studied the classification performance at the level of the whole ensembles of predictor quadruplets. We illustrate this briefly in Fig. S8E, even if a full analysis of the compared biomarking potential of different features falls beyond the scope of the present study. The figure depicts a precision-sensitivity analysis (see Supporting Materials and Methods) of the capabilities of every predictor quadruplet to discriminate our subjects into two different categories (age lower or higher than the median dataset age). All the distributions of precision-sensitivity pairs were unimodal and concentrated around their mean. SC ij quadruplets tended to achieve better precision in discrimination, but MC(ij) quadruplets tended to be more sensitive. Interestingly, the peak of the precisionsensitivity distribution for mixed-type quadruplets was displaced toward both the best precision peak of SC ij quadruplets and the best sensitivity peak of MC(ij) quadruplets, hinting again at synergies between different predictor types, at the entire ensemble level and not only limited to the best single quadruplet.
Even smaller errors in age prediction and precision-sensitivity pairs in age class discrimination could have been achieved by using quadruplets of general MC ij,kl strengths rather than just integrated strengths MC(ij). However, we chose not to implement this here to guarantee a fair comparison with SC, FC, and σ(FC) features. While there are O(N 2 ) independent links -as well as MC(ij) integrated strengths-, there are O(N 4 ) possible independent meta-links, which poses a serious multiple comparison problem when assessing relative performances. These problems can be avoided by considering only integrated meta-link strengths MC(ij).

Experimental procedures
Subjects. Overall N = 85 healthy adult subjects (N = 53 females, N = 32 males) were voluntarily recruited at Charité -Universitätsmedizin Berlin to participate in rs fMRI and DSI scans and, for a subset of them, also in a motor study. The first subset of N = 49 subjects ('rs-only') had ages uniformly distributed over the 18-80y range. The second set of N = 36 subjects ('rs+tasks') was further split into a first (N = 15, 20-25 yrs) and a second (N = 21, 59-70 yrs) age groups. All subjects had no self-reported neurological, psychiatric or somatic conditions. For the 'rs+tasks' subset, healthy cognitive function was furthermore assessed with the Montreal Cognitive Assessment (MoCA) (24). For all the analyses of Figs. 1-2 and S1-8, in which cognitive performance was not relevant, we merged the two subsets of subjects. We distinguished for inter-group comparisons between a 'Young group' composed of subjects in the first age-quartile of the 'rs-only' subset and the first age group of the 'rs+tasks' subset (overall N = 28, 18-25 yrs, median age = 22.5 yrs), and an 'Elderly group' composed of subjects in the fourth age-quartile of the 'rs-only' subset and the second age group of the 'rs+tasks' subset (overall N = 33, 59-80 yrs, median age = 65 yrs).
In addition to general exclusion criteria for participation in an MRI experiment, subjects with a self-reported musical background were excluded, as musical training may affect the performance of the rhythmic motor task. Left-handed subjects, identified using the Edinburgh Handedness Inventory, were also excluded. Subjects were informed of the procedure of the study and basics of fMRI acquisition, and written consent was obtained prior to data collection. The study was performed in accordance with the local medical ethics committee protocol. Here, FEAT (fMRI Expert Analysis Tool) first-level analysis from the FMRIB (Functional MRI of the brain) software was used. Motion correction was performed using EPI field-map distortion correction, BET brain extraction, and high-pass filtering (100s) to correct for baseline signal drift, MCFLIRT to correct for head movement across the trial. We checked that additional head movement correction by regressing out the six FSL head motion parameters did not significantly modify correlations with age of FCD and MC markers extracted from our BOLD time series. Functional data was registered to individual high-resolution T1-weighted images using linear FLIRT, followed by nonlinear FNIRT registration to Montreal Neurological Institute MNI152 standard space. Voxel-level BOLD time series were reduced to 68 different brain region-averaged time series, according to a Desikan parcellation (27). See table S1 for the regions under consideration. We did not perform slice-timing correction, smoothing, or normalization of BOLD intensities to a mean.
Resting state. During resting state scans, subjects were to remain awake and reduce head movement. Head cushions were used to minimize head movement, and earplugs were provided. Scans for the 'rs-only' and the 'rs+task' subsets of subjects had different durations. For the 'rs-only' subset, 20m of uninterrupted rs scan were performed. For the 'rs+task' subset, five minutes of rs were collected before the task acquisition (see later), and then further five minutes after the task. Visuo-motor coordination task (frequency-locking Φ score). The N = 36 subjects in the 'rs+tasks' subset also succeeded in performing a visuo-motor coordination task while in the scanner. The task followed a unimanual paradigm, which was adapted from a bimanual paradigm introduced in (36). During the task, subjects were told to lay still inside the scanner with an air-filled rubber ball in their right hand. A screen, animated using a custom-made LabView program, was projected in the scanner. To reduce eye movement, subjects were instructed to fix their gaze at a cross, displayed in the middle of the screen between two rotating disks. The left disk served as visual cue, rotating at a computer-generated speed, while the subject's squeezing of the ball controlled the speed of the right disk. The goal was to make the subject-generated rotating disk align in (counter) rotation with the computer-generated rotating disk, which was done by squeezing the rubber ball in a 4:3 frequency to the visual cue. For perfect performance, the two disks would rotate in synchrony. Because the computer-generated disk rotated at a 4:3 frequency to the subject-generated circle, subjects had to squeeze the ball at 1.35 cycles per second to match the 1.8 cycles per second of the computer-generated disk in order to achieve synchrony Behavioural measures were collected (1 performance score per trial) based on the frequency locking of the two rotating circles. If the two disks rotated perfectly insynchrony (i.e. subject was able to match the frequency of bulb-squeezing to the computer generated cue), the performance score would be 1. Not frequency-locked rotations of the two disks would result in a performance score of 0. More specifically, the frequency locking of the computer-generated circle and the subject-generated disk was quantified by the normalized spectral overlap between the power spectra of the two forces, P x and P y , as described in detail in (37): with ρ = 4/3, corresponding to the target frequency ratio between the two rotating disks. Behavioural performance was expected to improve across trials as subjects learned the task, but in this study focused just on the average performance over the ten trials, ignoring learning.
MoCA cognitive assessment. MoCA assessment was performed by N = 21 elderly subjects of the 'rs+tasks' subset. The MoCA includes multiple sub-tasks probing different cognitive domains such as: short-term memory and delayed recall; visuo-spatial abilities; phonemic fluency, verbal abstraction and naming; sustained attention and concentration; working memory; executive control in task switching; spatio-temporal orientation. The test was administered in a German version (downloadable from http://www.mocatest.org). The maximum global score 'MoCA' achievable is of 30 points, up to 5 of which are contributed from the partial score 'MoCA-wm' from the working memory ('Erinnerung') task. Participants were considered in good/healthy mental state, when achieving scores higher than 25. All details can be found in (24).

Network analysis
Extraction of time dependent Functional Connectivity. As, e.g., in (12), time-dependent Functional Connectivity matrices FC(t) were estimated by sliding a temporal window of fixed duration τ and by evaluating zero-lag Pearson correlations between rs BOLD time series a i (t) from different brain regions: where temporal averages are taken over the restricted time interval [t, t + τ]. All entries were retained in the matrix, independently of whether the correlation values were significant or not or without fixing any threshold (i.e., we treated FC ij entries as descriptive features operationally defined by the above formula).

Evaluation of Functional Connectivity Dynamics (FCD) matrices.
We introduced a notion of similarity between any two FC(t 1 ) and FC(t 2 ) matrices following (19), based on the Pearson correlation between the entries of their upper-triangular parts: FCD matrices depend thus on the window-size τ adopted when extracting the FC(t) stream.
Evaluation of FCD step lengths and their distribution. Correlation distance provide a natural (bounded) metric for the space of FC(t) matrices. We can measure the distance traveled in FC space when stepping from one matrix FC(t) to the next matrix in a sequence by evaluating: Since the time step τ is maintained constant, these step lengths can be interpreted as rates of FC reconfiguration across time, and, thus, FCD speed.
It was difficult to estimate histograms of d τ at the single subject level based on a single window-size τ because the number of sampled points within a single rs fMRI session was too small, especially for the longer window-sizes. Observing that the single-subject histograms that we obtained were similar for close window sizes, we chose to group observations obtained for different scales into three size-groups: 'slow' window-sizes (from 4m to 1m), 'normal' window-sizes (from 1m to 20s) and 'fast' window-sizes (20s to 8s). The exact choice of included τ within each size-group was done to have τ homogeneously distributed over the size-group range and to have, after pooling, the same number of d τ observations available for building each of the three histograms (there are hence larger ranges for the 'slow' and 'normal' size-groups, since the longer the τ the fewer non-overlapping windows could be allocated within the fixed length of the imaging session). Histograms for the three size-groups averaged over all subjects can be seen in Fig. S3A, together with a group-averaged histogram for ultra-fast window-sizes (9s to 6s), which still has a 'well-behaved' shape. Single-subject histograms for the 'fast' size-group are shown in Fig. 1B. Confidence intervals for the histograms were evaluated according to a standard Agresti-Coull binomial proportion approximation.
The typical FCD speed d typ was determined as the mode of a single-subject windowsizes pooled d τ histogram and was as such dependent on the chosen size-group. When there were multiple d τ values associated to the same maximum count then d typ was given by the average of these multiple d τ values.
Histograms under the null hypothesis of lack of sequential correlations were computed as above, apart from the following change. Before computing the d τ values, the order of the observed FC(t) matrices in the observed stream was randomly permuted. In this way the FC mean and variance were preserved for each possible link, but the eventual sequential correlations were disrupted by construction. The randomization of order was performed 500 times and independently for each of the window sizes pooled to build the overall null hypothesis sample. With respect to null hypothesis expectations, single subject histograms very frequently showed: an under-representation of long step lengths; an over-representation of short step lengths; or both. We studied how frequent these two types of discordancies from chance-level expectations were by computing d τ histogram made of 20 d τ bins for all subjects and all window-sizes τ (no longer pooled). We ignored the actual values of the center of the twenty d τ bins (which were different for each of the τ-s) and only considered their ordinal rank -always 1 to 20 for all τ-s -, in order to systematically compare across different window-sizes τ in which bins (relatively short or relatively long step lengths) an over-representation or an under-representation of counts were significantly occurring. The results of this analysis are shown in Fig. S3C.
Detrended Fluctuation Analysis (DFA). The Detrended fluctutation analysis (DFA) allows us detecting intrinsic statistical self-similarity embedded in a seemingly nonstationary time series. It is particularly adapted to the study of time series that display long-range persistence, and it is in this sense similar to other techniques, such as Hurst exponent analysis, requiring however the stationarity of the analyzed signal. See (29) for a review. In order to capture self-similarity and auto-correlations among increments in a time series, DFA infers a self-similarity coefficient by comparing the detrended mean square fluctuations of the integrated signal over a range of observation scales in a log-log plot. If the log-log plot has an extended linear section, (i.e. if the scaling relation is a genuine power-law over a reasonably broad and continuous range of scales, see later for the meaning of 'genuine'), it means that fluctuations 'look the same' across different temporal scales, i.e. we have statistically the same fluctuations if we scale the intensity of the signal respecting the DFA exponent.
We performed here DFA over the time series of d τ 's to detect sequential autocorrelations within them. First the sequence d τ (t 1 ), d τ (t 2 ), … d τ (t L ) was converted into an unbounded process: Prior to construing outcome values, however, it is mandatory to verify that a linear power law scaling actually exists. If it was not the case indeed the output value α DFA could not be interpreted as a scaling exponent, but just as yet another meaningless number.
Following (30), we tested the hypothesis of power-law scaling using a (Bayesian) model comparison approach. This allowed identifying the subjects for which the DFA log-log plot was better fit by a straight line than by any other tested alternative model. Only these subjects with a proper linear section in the DFA log-log plot were retained for the following steps of DFA exponent extraction and analysis of correlations with age.
In order to test the hypothesis of power law against alternative models, we evaluated the density of fluctuations over the consecutive segments, i.e. the density of ! beyond its mean value ! using a kernel source density estimator. Based on this probability density, one can estimate the log-likelihood for a certain model to generate fluctuations of a given strength (on a log-scale) as a function of log . To perform model selection, the toolbox then computes the corrected Akaike Information criterion for each one of the tested models: where p is the number of free parameters to fit in the model and the number of points used to estimated the density of ! . Note that this model selection criterion automatically embeds a penalization for models with larger number of parameters, thus protecting against over-fitting. The model yielding the lowest AICc was selected as the relatively best one, and in the case this was the linear one the corresponding ℒ !"# -fitting parameter was considered as !"# .
To improve fitting and model selection we increased the number of available datapoints by extracting more than just one d τ sequence from each given subject. Additional d τ sequences were constructed based on streams of stepped FCD windows with the same length τ, but with slightly shifted beginnings. In detail, in addition to the standard sequence d τ (t 1 ), d τ (t 2 ), … d τ (t L ), we also extracted also a second sequence d τ (t 1 + δt), d τ (t 2 + δt), … d τ (t L-1 + δt), a third sequence d τ (t 1 + 2δt), d τ (t 2 + 2δt), … d τ (t L-1 + 2δt), … and, finally, a last sequence d τ (t 1 + qδt), d τ (t 2 + qδt), … d τ (t L-1 + qδt) with qδt < τ and (q+1)δt ≥ τ. We subsequently performed independent DFAs for each one of these shifted sequences and merged the associated clouds of k vs DFA(k) points prior to model selection.
We used a range of 5 < k < 20, to discard data chunk sizes that were too short or long data chunk sizes yielding an overall number M of chunks that was too small.
Note that our DFA exponents still depend on the FCD window-size τ, which also provides the actual unit of measure for DFA scales. The k values of d τ 's included in each DFA chunk are evaluated based on (k+1)τ s of BOLD rs activity. Therefore, the robustness of DFA results must still be verified across different FCD window-sizes, as shown in Fig. S4. We based the analysis shown in Fig. 1E-F with a window-size of τ = 12 s (and hence estimated the scaling exponent α DFA over time-scales between 60 s and 240 s). In Fig. S4B the significance of age-to-α DFA correlation at the single-subject level is shown also for different τ-s via permutation testing. The significance of intergroup differences between α DFA values for different τ-s was assessed via the U-Mann Whitney non-parametric two-sided test (Fig. S4C). All tests were performed with different group sizes at every considered τ value, since different subjects could be considered or not to show a genuine power-law scaling for different τ-s.
A genuine power-law scaling in the DFA of subjects within the 'rs+task' subgroup could be established only for very few subjects and not always the same for different τ-s. This was probably due to the shorter length of rs BOLD time-series acquired for the 'rs+task' subset of subjects. Therefore, we limited DFA to the 'rs only' subset of subjects and correlations of α DFA with cognitive and motor performance could not be tested.
FC modular flexibility analysis. To study how the modular structure of FC(t) matrices changes across time as an effect of ongoing FCD, we performed a modular analysis of the graph structure of each FC(t) matrix separately.
We used the Louvain algorithm to extract the modules of each of the networks FC(t a ) separately. This involves an iterative scheme in which modules obtained at a given iteration step were used as the initial condition for the next step until when the obtained modularity goal function converged. The procedure was repeated 500 times and the module partition associated to the best score was maintained. We denote the obtained modules as M 1 (t a ), M 2 (t a ), …, M K(a) (t a ) where K(a) is the number of FC modules extracted at time t a . Since modules were generated independently for each considered timewindow, there was no guarantee that a module M k (t a ) was related to a module M k (t b ) with which it -contingently -shared the same module index. For instance, it might be that the module at time M 1 (t a ) has a larger overlap with the module M 4 (t b ) than with the module M 1 (t b ). We therefore relabeled the module indices in order to guarantee the larger overall similarity across time between all modules sharing a common label. This relabeling enabled us to meaningfully track the displacement of brain regions between different modules. To align the retrieved modules, we determined inter-module distance as the number of elements in the set-theoretical symmetric difference between them: where: This was followed by a modified K-means clustering of the modules: i) choosing K, the maximum number of modules as the maximum K(a) observed across time in the rs session for each given subject; ii) by using the above definition of inter-module distance as metric for the clustering; and, iii) by imposing that at each step of the K-means iterative algorithm (including the final one after convergence) every temporary cluster always included one and only one module deriving from each of the time-windows. After convergence of the clustering algorithm, all modules within each of the clusters of modules were reassigned a common label corresponding to the index of the K-means cluster. For instance, a given cluster of modules Q could contain modules Q = {M 4 (t 1 ), M 2 (t 2 ), M 5 (t 3 ), … M 3 (t max -δt), M 1 (t max )}. All these modules were then relabeled as {M Q (t 1 ), M Q (t 2 ), M Q (t 3 ), … M Q (t max -δt), M Q (t max )} after our matching via re-clustering procedure. For the analyses underlying Fig. S5 we adopted δt = 1.94 s (corresponding to one TR). The module labels indicated in the top panel of Fig. S5A correspond to the labels Q after re-clustering of modules. The middle plot gives the evolution over time of K(a) for a specific subject.
Next, we estimated the fractions of brain regions changing module when moving from a time t to the following time t + δt. The number of region-switching events affecting a given module Q was given by: We evaluated the time-dependent 'liquid fraction' of regions as the following quantity (averaged over the modules present at time t): For a specific subject this quantity is plotted in Fig. S5A's bottom panel. Figs. S5B and S5C show correlations between subject age and the averages over time of the number of modules K(t) and of the liquid fraction Λ(t), respectively.
MC analysis. Out of the stream of FC(t) matrices, we extracted M = N(N-1)/2 time series of pairwise FC couplings given by the entries FC ij (t) for all pairs of regions i and j with i < j ≤ N. We used a window-size of τ = 15 TR (~29 s) with a window step of one TR (1.94 s). The entries of the meta-connectivity matrix MC were given by the conventional Pearson correlation values: We compiled these correlation entries into a format that allowed us to easily identify the pair of links involved into each meta-link and the participating brain regions. That is, we built MC matrices of (N 2 -N)-times-(N 2 -N) size, where different rows corresponded to different directed pairs of regions -i.e. both the pair (i,j) and the pair (j,i) were includedand only links corresponding to self-loops -i.e. of the type (i,i) -were excluded. The price to pay for this choice of clarity was to introduce redundancy among MC matrix entries, since: For the sake of computational efficiency, however, we computed only P = M(M-1)/2 independent Pearson correlations between pairs of possible FC link time series FC ij (t) and FC kl (t), with i < j, k < l, i ≤ k, j < l and copied this inter-link correlation value into the eight degenerate MC matrix entries.
Besides the strengths MC ij,kl of individual meta-links, we also computed integrated meta-connectivity strengths for individual links This quantity, depending uniquely on the two indices i and j, and not anymore on the four indices i, j, k and l, measures the total strength of all inter-link interactions involving a given link (ij). We could compute an analogous integrated strength for each brain region i : We limited in this case the sum to 'trimer' meta-links only, having the considered region i as 'root'. We call 'trimers' these special meta-links MC ikil which give the strength of interaction between two functional links (ik) and (il) that share the common root region i. Trimer meta-links involve therefore a total number of three different brain regions. Trimers are thus distinct from general tetramer meta-links, which involve two links touching a total of four different brain regions.
In the adopted format, the MC matrix can be considered as the adjacency matrix of a graph whose nodes are FC links. Thus communities of temporally co-varying FC links can simply be evaluated by extracting node-based communities out of the MC matrix. To extract MC modules, we adopted, as for the modular structure of FC(t) matrices, a standard Louvain algorithm, and selected the best modularity solution out of 500 independent runs. We extracted the modular structure in Fig. 2A out of the group average of single subject MC matrices taken over the 'Young' group of subjects (see above the 'Subjects' subsection for information about this group's composition).
We found that negative strength meta-links are particularly informative about age. Hence, we defined where θ(x) denotes the Heaviside function which is equal to 1 if x ≥ 0 and equal to 0 otherwise. In this way, neg-MC gives the sum of all the negative entries of the MC matrix.

Age prediction
Prediction via linear regression. We took into account a large number of different possible features as potential predictors of age. These predictors included the strengths of structural connectivity links SC ij or of time-averaged functional connectivity links FC ij , their variance over time σ(FC ij ) and their integrated meta-connectivity strength MC(ij). We also considered node-based features, such as the total strengths of a region under structural connectivity and functional connectivity, i.e. SC(i) and FC(i), and regional trimer meta-connectivity strengths MC(i). As shown in Fig. S7, for every possible brain region i and link (ij), we tested whether these quantities correlated with subject age and whether correlations were significant (via bootstrap with replacement, p < 0.05).
We performed age-prediction based on different sets of features. As predictor variables we used: individual SC ij strengths; individual FC ij strengths; individual σ(FC ij ); individual MC(ij); or quadruplets of SC ij 's; quadruplets of FC ij 's; quadruplets of σ(FC ij )'s; quadruplets of MC(ij)'s; and, finally, mixed quadruplets composed of SC ij , one FC ij , one σ(FC ij ) and one MC(ij) features. For each of these types of predictor sets, we tested all the possible individual features or feature quadruplets, keeping track of their performance in age-predictions and then isolated the best individual link or the best quadruplet of links for each of the predictor sets types.
For each chosen set of predictor features, prediction performance was evaluated following a cross-validation procedure. First, we randomly split the subjects into two agebalanced halves by including within each of the halves subjects belonging to all four age quartiles in equal number. One half was designated as the 'training' set; the second half was designated as the 'testing' set. Then, we fitted a linear relation between age and the chosen predictor set (standard least-square procedure) based on the subjects in the training set, and evaluated this linear relation over the subjects in the testing set to extrapolate their age. Subsequently, we compared the predicted age with the actual age, and evaluated absolute errors of prediction |age real -age predicted | and averaged them over the testing set to obtain the mean absolute error (MAE). Distributions of the expected MAE were then built by replicating 1000 times the random splitting of the overall sample into training and testing subsets. The boxplot in Fig. S8A shows the median, interquartile range and 95% confidence interval for the MAE distributions associated with the best individual link or quadruplet of links for each of the predictor set types. The best predictor sets achieving these performances are reported in Figs. S8B-D.
Discrimination into age-classes. Next to the best features for each of the predictor types, we studied the general performance of all the different features in a simpler task, in which only the age-class of a subject had to be predicted, rather than her/his actual age. Subjects were split into two classes, age 'lower than' and 'higher than median'. We assumed as target class for evaluating true or false positives the second one of the older subjects. Again subjects were separated into training and testing sets. A linear regression was then performed over the training set only. Ages for the testing set were predicted using these fitted linear relations. Finally, an age class was inferred depending on whether the predicted age was lower or higher than the actual whole sample median, used as a discrimination threshold. We evaluated performance metrics in terms of sensitivity (a.k.a. recall): and precision: where TP denotes the number of true positives, FP the false positives and FN the false negatives in age-class discrimination over the testing set. Distributions of Sensitivity and Precision were finally built for each predictor quadruplet, by repeating the random split into training and testing subsets 1000 times. Fig. S8E depicts the joint histograms of the achieved median Precision and Sensitivity pairs sampled over all considered quadruplets of predictors (4x SC ij , 4x FC ij , 4x σ(FC ij ), 4x MC(ij) and four mixed types).