State-switching and high-order spatiotemporal organization of dynamic functional connectivity are disrupted by Alzheimer’s disease

Alzheimer’s disease (AD) is the most common neurodegenerative illness, with an estimated prevalence of 10–30% in people older than 65 years (Hou et al., 2019; Masters et al., 2015). Yet, despite substantial research, we are far from fully understanding the mechanisms that link pathophysiology to cognitive impairments. Neurodegeneration in AD has been traditionally associated with the extracellular accumulation of insoluble amyloid-β₄₂ (Aβ) neuritic plaques (Glenner & Wong, 1984; Lemere et al., 1996) along with the intracellular accumulation of abnormally phosphorylated tau (pTau) (Aβ and Tau), which constitute the neurofibrillary tangles (Spires-Jones & Hyman, 2014). These processes yield to widespread neuronal death, synaptic loss, and atrophy (Bateman et al., 2012), with a progression of structural damages not occurring uniformly throughout the brain (Braak & Braak, 1991). However, the progression of neurodegenerative processes does not correlate linearly with the severity of cognitive impairment possibly due to a “cognitive reserve” accrued through education, cognitive training, and a healthy lifestyle (Rentz et al., 2010; Snowdon, 2003). Furthermore, the severity of cognitive impairment symptoms in a patient can fluctuate substantially within the same day, faster than the time scales of neurodegeneration (Palop et al., 2006). Together, these findings suggest that AD involves alterations of neural dynamics and that these dynamical changes may be the mechanistic substrate leading to functional impairment or preservation.

As molecular and structural changes alone do not fully account for cognitive impairment, alternative studies based on functional connectivity (FC) analyses have sought to fill the gap. In particular, resting-state FC (Fox & Raichle, 2007) quantifies brain-wide correlations of BOLD signals, capturing interactions between regions. In this context it has been suggested that structural alterations in AD lead to FC changes (Dennis & Thompson, 2014), and that the early manifestation of Aβ toxicity preceding overt atrophy can be detected using resting-state functional magnetic resonance imaging (rsfMRI) (Hedden et al., 2009; Sheline et al., 2010a, 2010b; Mormino et al., 2011). Changes in FC in AD include reduced connectivity within the default mode network (Fleisher et al., 2009; Greicius et al., 2004; Jones et al., 2011; Petrella et al., 2011; Rombouts et al., 2005; Sorg et al., 2007; Wang et al., 2006, 2007; Zhang et al., 2009, 2010), in a spatially nonuniform fashion (Damoiseaux et al., 2012). Besides Aβ, the deposition of pTau affects FC as well (Franzmeier et al., 2022). Furthermore, additional FC alterations have been reported, leading to functional disconnection between hemispheres (Shi et al., 2020; Wang et al., 2015) and a reduction of small-world topology (Brier et al., 2014; Sanz-Arigita et al., 2010; Stam et al., 2007, 2009; Supekar et al., 2008).

More recently, investigations of FC in AD have been extended to encompass time-varying, rather than time-averaged FC. Indeed, rsfMRI networks undergo a continuous reconfiguration of their weighed topology, and the statistical structure of spontaneous network reconfiguration carries information potentially useful to discriminate cohorts (Calhoun et al., 2014; Hutchison et al., 2013; Preti et al., 2017). The flexibility of dynamic functional connectivity (dFC) has been shown to correlate with cognitive performance (Bassett et al., 2011; Battaglia et al., 2020; Braun et al., 2015; Jia et al., 2014; Lombardo et al., 2020; Shine et al., 2016). In this view, ongoing variability of FC networks is not noise but rather an actual resource subserving computation. The capacity to actively maintain a spatiotemporally organized yet variable dFC would confer the system resilience to cope with variable cognitive and environmental conditions (Lombardo et al., 2020). Hence, the preservation of a “healthy” structured dFC variability may provide a form of functional compensation and a likely neural substrate for “cognitive reserve” (cf. also other studies linking mental training with enhanced dFC variability, e.g., Premi et al., 2020). Conversely, dynamic FC-based metrics thus promise to better characterize the impact of AD pathology.

A number of studies have quantified dFC changes in healthy aging (Battaglia et al., 2020; Davison et al., 2016; Hutchison & Morton, 2015; Lavanga et al., 2022; Petkoski et al., 2023; Qin et al., 2015; Viviano et al., 2017) and in conditions such as schizophrenia (Damaraju et al., 2014; Sakoğlu et al., 2010), epilepsy (Liao et al., 2014; Liu et al., 2017), and Parkinson’s disease (Fiorenzato et al., 2019; Kim et al., 2017). In AD, probabilities of temporal transitions between alternative FC states have been shown to be altered (Jones et al., 2011; Fu et al., 2019; Gu et al., 2020; Schumacher et al., 2019). Moreover, machine learning applications have achieved greater accuracy in differentiating between healthy control and aMCI or AD subjects when trained with dFC-based rather than static FC metrics (Chen et al., 2016, 2017; de Vos et al., 2018; Wee et al., 2016). Although the contributions of these studies are promising, they are largely descriptive and do not propose an explicit theory of why dFC changes lead to functional consequences. Furthermore, the plethora of methods for dFC quantification (Hutchison et al., 2013; Preti et al., 2017)—from extracting discrete FC states (Allen et al., 2014; Thompson & Fransson, 2016) to continuously time-resolved approaches (Battaglia et al., 2020; Lindquist et al., 2014)—hinder the convergence of results.

Here, we start from a theoretical tenet: efficient cognition requires spatiotemporally organized FC variability, which is neither trivial, nor random, but complex. This assumption is based on empirical evidence. Fluctuations in dFC are not a mere unstructured “Drunkard’s walk”: more highly structured dFC trajectories are observed in individuals with higher performance on general cognition domains (Battaglia et al., 2020; Lavanga et al., 2022). Furthermore, individual FC links do not fluctuate independently but with network reconfigurations governed by higher order coordination patterns, manifest by nontrivial interlink covariance patterns (Davison et al., 2015; Faskowitz et al., 2020; Petkoski et al., 2023); “back-bones,” partially scaffolding dFC (Braun et al., 2015); and dFC flowing under the influence of competing “meta-hubs” (Lombardo et al., 2020). Reiterating, our hypothesis suggests that spatiotemporal structure of dFC between order and randomness allows for rich computation to emerge from the systems’ collective activity (cf. Crutchfield, 2012). Correspondingly, we predict that individuals with higher cognitive performance should display an enhanced organization of dFC compared to those with impaired cognition (aMCI or AD) in which, conversely, a loss of dFC spatiotemporal organization should be evident.

Here we analyze resting-state fMRI data acquired from individuals with better-than-normal or normal cognitive performance—supernormal (SNC) and normal controls (NC)—and those clinically diagnosed with amnestic mild cognitive impairment (aMCI) or AD. We first characterized dFC across groups using two complementary methods. First, we use a state-based dFC analysis paradigm, in which we assume the existence of a small set of possible discrete FC configurations and quantify dwell times in different states and the temporal stability of different FC network links along state-switching transitions (Thompson & Fransson, 2016). Second, we use a state-free dFC analysis paradigm, where FC networks are described as continually morphing in time. Through these complementary but convergent approaches, as described in the following, we find that the fluctuations of different links show different degrees of mutual interdependence across the considered groups, shifting from a “liquid-like” dFC (flexible but constrained) for SNC and NC toward a “gas-like” dFC (uncorrelated and disordered) for patient groups. We also show that these changes in dFC coordination cannot be fully accounted by changes occurring at the level of ordinary pairwise FC, but stem from the weakening of genuine higher order interactions observed especially for regions that are among the first to be physiopathologically affected by AD.

We considered an fMRI dataset including resting-state sessions from subjects with varying degrees of cognitive skills. As our interest focuses not only on disease but also on healthy cognition, healthy controls were subclassified in two groups (SNC and NC) based primarily on composite memory Z scores to define the SNC and NC groups. That is, SNC had a higher performance in the composite memory scores (Z > 1.5) and at least a Z > 0.7 in all other cognitive domains (attention, language, visuospatial, and executive; see Materials and Methods for more details). Healthy control subjects between NC and SNC or below NC were not considered in the study. As shown in Figure 1A, from 73 subjects, 10 were classified as supernormal controls (SNC), 16 as normal controls (NC), 34 as amnesic mild cognitive impairment (aMCI), and 13 as Alzheimer’s disease (AD). Across the four clinical groups, there were no significant differences in age and sex.

Based on rsfMRI time series from these cohorts, we then computed (and compared across groups) a variety of static and dynamic functional connectivity (FC and dFC) metrics, extracted with complementary approaches, assuming or not the existence of discrete FC states in time (Figure 1B). Importantly, as detailed below, we did not uniquely consider pairwise interactions between two brain regions at a time, but also considered more complex coordination patterns between larger groups of regions. Classic FC links express the existence of a correlation between the BOLD fluctuations of two brain regions and are represented as a link between two regional nodes: we refer hence to them as dimers, since they are computed out of two parts. In classical FC analyses, dimers are static, as their strength is averaged over the duration of complete resting-state sessions. In dFC analyses, however, dimer strengths fluctuate in time. We can thus also compute correlations between different dimers. Estimating these “correlations between correlations” requires jointly monitoring the BOLD fluctuations of three (Figure 1C, top) or four (Figure 1C, bottom) regions, hence the names of trimers and tetramers—collections of three or four parts, respectively—used in the following.

We chose to focus in this study on dFC within a network of limbic brain regions of particular interest (Figure 1D). The rationale was twofold: first, the regions included in the chosen limbic subnetwork are highly interconnected brain regions that degenerate early in the disease process (Arnold et al., 1991; Braak & Braak, 1991); second, previous modeling work confirmed their central role in shaping the evolution of FC alterations comparing healthy controls to aMCI or AD stages (Zimmermann et al., 2018b).

In order to assess FC changes along time, we started with a state-based dFC approach, called the PBM and first introduced by Thompson and Fransson (2016). In this framework, different instantaneous images of brain-wide BOLD activation are first clustered via an unsupervised procedure into K states, and state-specific FC matrices FC^(λ) are constructed by evaluating BOLD correlations limited to time frames assigned to a given state cluster (λ = 1 … K, see Materials and Methods for details). Figure 2A shows the weighed adjacency matrices FC^(λ) (obtained as centroids of their respective cluster) for each of four different states of dFC, called S-graphlets by Thompson and Fransson (2016). An alternative graph representation of these templates is shown in Figure S1A. The optimal number of K = 4 was determined based on a statistical elbow criterion (Figure S1B) and confirmed post hoc by the consistency of our results.

Based on these four dFC states, we obtained the spatial profile of neural activation across regions (Figure 2A). The spatial organization of the observed neural activation profiles naturally suggests, in this study, to group the regions in two subsets, characterized by having an activity level transiently higher or lower than their average level. We defined Zone I as the subset of ventral limbic regions including amygdala, temporal pole (superior and medial), hippocampus, and parahippocampal gyrus (limbic system and posterior cingulate cortex). Zone II included the cingulate gyrus (anterior, medial, and posterior). In States 1 and 2, Zone II (dorsal regions) and Zone I (ventral regions) were respectively active above average level (high activation states). In contrast, in States 3 and 4, Zone II and Zone I regions were respectively active below average levels (low activation states).

Furthermore, these four states were noted based on the topology of their FC^(λ) networks and the level of internal synchronization within Zone I. Quantitatively, connection weights between regions within Zone I tended to be stronger for States 2 and 4 than for States 1 and 3 (average within Zone I FC weights = 0.23 ± 0.16 for States 1 and 3 vs. = 0.29 ± 0.18 for States 2 and 4). Hence, States 2 and 4 displayed higher internal synchrony, in contrast to States 1 and 3. Then we computed local and global efficiency metrics (Achard & Bullmore, 2007; Latora & Marchiori, 2001) for the four FC^(λ) networks. Global efficiency quantifies how well communication pathways can be established between any two nodes in a weighed network. Local efficiency quantifies the robustness of communication and the possibility to find alternative routes if local connectivity is disrupted. We found that the high sync States 2 and 4 have a lower global efficiency (Figure 2B; Mann-Whitney U test, p < 0.001) but a greater local efficiency (Figure 2B, Mann-Whitney U test, p ∼ 0.023), reflecting a denser within-zone but a weakened between-zone connectivity (average between Zone I and Zone II FC weights = 0.026 ± 0.069 for States 1 and 3 vs. = −0.013 ± 0.071 for States 2 and 4).

Thus, in short, the overall four states that we find are obtained as combinations of two qualitatively different network topologies and two possible levels of activation, so that each topology can exist in low and high activity versions.

We quantified the stability of dFC both by the longer or shorter duration of transient epochs within a given state (average dwell time, Figure 2C) and by the overall time fraction spent within a state (average state census, Figure 2D). As shown in Figure 2C, group differences were identified in the mean dwell time of low sync states, with longer dwell time for the two control groups (∼3.6 TR = 7.4 s, for SNC at one extreme) and shorter for the MCI and AD groups (∼2.8 TR = 5.7 s, for AD at the other extreme). However, the mean dwell time of high sync states were not different.

Analogously, Figure 2D shows that the relative fraction of time spent in low sync states decreased in aMCI and AD compared to healthy controls (ranging from 62% for AD to 72% for SNC).

In summary, low sync and globally efficient dFC states were less frequent and more transient in aMCI and AD, suggesting a reduction of their overall stability.

The next step, also following Thompson and Fransson (2016), was to map a state-based dFC temporal network to each subject’s resting-state acquisition. To do so, we constructed a sequence of network time frame FC(t) set to be equal to the FC^(λ) graph specific for the state λ visited at time t (Figure 3A; see Materials and Methods for details). Thompson and Fransson (2016) called such a temporal network a T-graphlet.

In this approach, each link can assume up to four possible strength values, corresponding to its strengths in the FC^(λ) associated to each of the four states. Hence, any variability of dFC dimers reflects exclusively state-switching dynamics. Figure 3B shows the time course for a representative fluctuating dFC dimer. The temporal organization of link fluctuations (periodic or bursty) can be highlighted by a binarization procedure, where a link is set to 1 if its instantaneous strength is above the threshold θ, or to 0 otherwise (see Materials and Methods). The result of this procedure is shown in Figure 3C for a few representative links and a specific choice of threshold. A link whose strength remains steadily above (below) threshold will result as constantly—or tonically—“active” (“inactive”). In contrast, a link whose fluctuating strength crosses the threshold through the different dFC-state frames will undergo several activation and inactivation events at specific threshold crossing times. Yet, there can be various types of intermittency, with different temporal statistical properties. The durations of different link activation and inactivation epochs could all be roughly similar, resulting in a more periodic type of intermittency (blue color link activation rasters in Figure 3C). Alternatively, they could be more variable, stochastically alternating between shorter and longer activation epochs (red color rasters in Figure 3C). The degree of temporal regularity in link activation and deactivation dynamics (link burstiness, mean and total activation, memory) can be evaluated, link-by-link, by the quantification of a burstiness coefficient (β). We also define the mean duration of a link’s transient activation events as mean activation (μ), and the total fraction of time in which a link is active relative to imaging session duration, total active time fraction (τ). The burstiness coefficient is bounded in the range −1 ≤ β ≤ 1, with β < 0, corresponding to near-tonic or periodic link activation dynamics, β = 0, corresponding to Poisson (random-like) link activation dynamics, and β > 0 corresponding to time-clustered (bursty) events of link activation. Mean activation times μ are bounded to the length of time series. Total active time fraction is also bounded, 0 ≤ τ ≤ 1.

In this approach, three numbers β (burstiness coefficient), μ (mean activation), and τ (total active fraction) fully characterize the binarized dynamics of a link (for a given choice of the strength threshold θ). These metrics were evaluated for the two categories of dFC dimers: intrazone (between two regions within either Zone I or II) and interzone (between one region in Zone I and one region in Zone II). Our results show that these two categories have distinct distributions of β, μ, and τ, first exemplified in NC subjects (Figure 3D). Whereas interzone dFC dimers are closer to a Poisson-like intermittency (β = −0.229 ± 0.020, median ± MAD), intrazone dimers, present a tonic activation time course (β = −0.890 ± 0.027, median ± MAD). In addition, interzone dimers are also less active (τ = 0.312 ± 0.099 for interzone vs. τ = 0.855 ± 0.027 for intrazone dimers) and activate for shorter transient times (μ = 34.926 ± 4.439 for interzone vs. μ = 178.995 ± 7.378 for intrazone dimers). These results suggest a smaller average strength of interzone time-averaged FC than for intrazone FC. Using NC subjects as reference group, we measure indeed an average FC(t) strength = 0.083 ± 0.135 for interzone and of 0.564 ± 0.155 for intrazone dimers (average ± SD). Similar differences were found for all groups (Table S1). The relative differences in β, μ and τ between intra- and interzone dimers are maintained over the entire range of possible thresholds θ (Figure S1C for bustiness coefficient). Interzone dimers also displayed more burstiness, were more transient and less active than intrazone dimers in all groups.

To achieve a robust and more precise comparison of β, μ, and τ distributions between the cohorts (Figure 3E), we computed percent changes of the three indicators in SNC, aMCI, and AD groups relatively to normal controls. The advantage of relative comparisons is that they can be collated for different threshold values θ, resulting in a threshold-independent analysis. We found that, moving from NC to aMCI and AD subjects, many dFC dimer links tended to have larger burstiness values. In contrast, moving from NC to SNC subjects, dFC dimers tended to be more tonic. These trends of β were smaller yet significant for intrazone FC dimers (Figure 3E), compared to interzone dimers, reaching +1.869 ± 1.663% for aMCI patients, +9.071 ± 3.001% for AD patients and −6.404 ± 1.938% for SNC subjects (Figure 3E) that had larger values.

These results reinforce the notion of a significant reduction of interzone time-averaged FC along the clinical spectrum (cf. Table S1). More importantly, and beyond this reduction of average strength, our results point to a degradation of the temporal regularity of FC fluctuations. While the total active time fraction τ of interzone dFC dimers decreased by less than −2% from NC subjects to aMCI and AD patients (Figure 3E; and even increased for intrazone dimers), the burstiness of interzone links increased over 10%, showing a real alteration in the temporal statistics of link activation, well beyond the trivial decrease necessarily induced by the observed reduction of average strength.

We also observed a significant decrease of the mean activation time μ (Figure 3E), for both intrazone and interzone dFC dimers (−1.275 ± 0.227% for aMCI and AD subjects compared to NCs). For SNC relative to NC, however interzone link burstiness decreased and their activation time increased (+0.613 ± 0.161% for SNCs).

Goh and Barabási (2008) also defined another metric related to burstiness, the memory coefficient. This coefficient λ (see Materials and Methods for exact definition) becomes significantly positive when autocorrelation exists in the duration of consecutive link activation events, that is, when long- (short-) lasting activation events tend to be followed by activation events which also are long (short). Computing λ, we found a weak median autocorrelation in all four groups, for both intra- and interzone links. Values (see Supporting Information Table S2) were small but still significant given the large number of activation events. Furthermore, memory was decreasing across the four groups from SNC to AD, providing yet another indication of increased disorder.

In summary, the temporal dynamics of dFC dimers between regions in different zones are altered along the SNC-AD spectrum from tonic and periodic in SNC to more intermittent in aMCI and AD subjects. Together with the finding of altered dwell times and transition dynamics between dFC states (Figure 2C and D), our state-based dFC analyses based on the PBM approach suggest that changes toward AD involve a degradation of global integration and an increased disorderliness of dynamic functional interactions between zones.

The PBM approach to dFC analyses reduces the description of FC network reconfiguration to the tracking of discrete state-switching events. Alternatively, sliding-window approaches evaluate the evolution of FC links as a continuous reconfiguration along time. As shown in Figure 4A, all dFC dimers FC(t₁) can be evaluated in a time-resolved manner restricting their estimation to BOLD signal time series within a window centered at time t₁. The window is then shifted at a slightly increased time t₁ + δt, providing an updated set of values FC(t₁ + δt). The result is a collection of smoothly varying continuous time series FC(t) for each possible dFC dimer (Allen et al., 2014; Battaglia et al., 2020).

As in the case of node activity time series, it is possible to study covariance between the temporal evolutions of different dimers. The case in which their fluctuations are not independent—or, in other words, that the dimers are “entangled”—will be signaled by significantly positive or negative correlations between dimers. These correlations can be represented graphically by trimer and tetramer diagrams in which the two entangled dimers are linked by a spring (Figure 4A, top left; we will omit in the following to draw this spring, for the sake of a clearer visualization). The stronger the correlation between the fluctuations of different dFC dimers in a trimer or tetramer, the stronger will be their “entanglement” (i.e., metaphorically, the stiffness of the spring).

These strengths of entanglement between FC dimers can be compiled into a meta-connectivity matrix (MC; Figure 4A). The notion of MC (Lombardo et al., 2020) is strongly related to the edge-centric FC discussed by Faskowitz et al. (2020). The key difference is that MC is obtained by using a short smoothing window in the estimation of the stream of FC(t) matrices, while edge-centric connectivity captures coincidences between instantaneous fluctuations. The denoising brought by the smoothing window allows an easier extraction of the modular structure of MC, with respect to edge-centric FC (cf. Lombardo et al., 2020), but the two concepts are otherwise equivalent. The choice of window size (here 5 TRs, Materials and Methods) was motivated by the fact that the state-based PBM method suggested that ∼90% of epochs within a coherent state lasted less than 5 TRs (Figure S2A), indicating a fast intrinsic time scale of link fluctuation. Furthermore, we can observe post hoc that the use of larger (or smaller) windows would not improve the capability to separate our groups based on MC values (Figure S2B).

Group-averaged MC matrices are shown in Figure 4B for the four groups. Their modular structure is evident at simple visual inspection. A module in the MC matrix—also called dFC module or meta-module (Lombardo et al., 2020)—corresponds to a set of cofluctuating dynamic FC links, that is, to FC subnetworks whose overall strength waxes and wanes transiently along the resting state in an internally synchronous manner. The existence of nonuniform MC matrices indicates that the flow of dFC reconfiguration is not mere noise but rather, it is organized by specific arrangements of “springs between the links.” In other words, fluctuations of FC dimers are entangled in complex patterns reflecting higher order correlations (nonvanishing trimers and tetramers) between the coordinated activation of multiple regions.

MC matrices can also be represented as graphs, in which MC nodes correspond to different FC links and MC links appear due to the entanglement between the FC links. An example graph embedding is shown in Figure 4C for the MC matrix of the NC group. Graph vertices are color coded depending on the type of associated FC link (i.e., start and end zones of the links; cf. FC diagram with matching colors at the top right of Figure 4C). Notably, the different dFC modules, visible as blocks in the MC matrices of Figure 4B and as uniform-color node communities in the graph of Figure 4C, are composed of FC dimers with internally homogeneous start and ending zones.

A standard graph-theoretical notion useful when commenting about dimer arrangements into trimers and tetramers is the one of incidence: a link is incident to a node (or a node incident to a link), if the link is attached to the node (the notion of incidence complements the more familiar one of adjacency, where two nodes are said to be adjacent if connected by a link). Equipped with this terminology, we call root the common region incident to both the dimers within a trimer, while the other two regions form the leaves of the trimer. We can then describe the first three dFC modules (1, 2, and 3) of the MC matrix as including mutually entangled FC dimers originating in either one of the Zone II cingulate regions and terminating in Zone I. The entanglement of FC dimers gives thus rise to strong interzone trimers with “roots” in Zone II and “leaves” reaching out to Zone I regions. The two other dFC modules 4 and 5 include dimers within Zone I and Zone II, respectively, forming strong within-zone trimers or tetramers. Entanglement is thus particularly strong between dimers within a same zone and between interzone dimers incident on a common root region (in Zone II).

Although the MC graph is highly modular, it is not split into disconnected components and some entanglement exists also between dimers located in different dFC modules. Intermodule connections in the MC graph can arise, for example, due to the existence of trimers with a root in Zone I (entangling dimers across dFC modules 1, 2, and 3) or interzone tetramers (entangling dimers across dFC modules 4 and 5). In other words, MC reveals some degree of global, widespread entanglement between FC dimers, beyond modular entanglement. The strength of such global entanglement is quantified by the so-called average participation coefficient of the MC matrix, a graph-theoretical quantity measuring intermodule coupling (Guimerà & Amaral, 2005; see Materials and Methods).

The distribution of MC participation coefficients for each group are shown in Figure 4D. We found that the participation coefficients decreased significantly (Figure 4D, left; Mann-Whitney U test, p < 0.001) from SNC to AD, while overall modularity did not vary significantly (Figure 4D, right). These results suggest that, in patients, coordination structure between fluctuations of FC dimers is impoverished: global entanglement is disrupted, making dimer fluctuations in different modules more random and mutually independent.

Before entering a more detailed and regional specific account of changes to dFC organization observed at the regional level along the SNC-to-AD spectrum, it is important to stress that trimer and tetramer analyses are not redundant with the dimer-based analyses. Indeed, some studies have suggested that correlation between edges (captured by higher order trimer and tetramer in a MC matrix) could just be an automatic byproduct of existing lower order dimer interactions (Novelli & Razi, 2022). This can be easily understood through some examples. Let us consider, for instance, two strong dimers FC_ri and FC_rj sharing a common root region r. If a third strong dimer FC_ij also exists—closing the triangle of edges (ri), (rj), (ij)—then it is not surprising that a strong trimer MC_ri,rj is also detected: indeed, the fluctuations of the two leaf regions i and j are coordinated through a transverse dimer interaction, that is, the strength of the trimer would be the byproduct of a triangular motif of dimers and would thus be a redundant consequence of them. Analogously, we may consider the case of a square motif of dimers FC_ij, FC_jk, FC_kl, and FC_li, which could also give rise to strong tetramers because of the presence of one or more pairs of strong dimers. In other words, the detection of strong trimer and tetramer entries within the MC (or other forms of edge-centric FC) is not a sufficient condition for the existence of genuine high-order interactions (Battiston et al., 2020) that cannot be explained as stemming from motif arrangements of lower order pairwise interactions. On the contrary, the existence of genuinely high-order interactions could be established by detecting trimer or tetramer couplings between the dimers in a motif, stronger than the dimers themselves involved in the motif. The question that then arises is, what is the structure of MC that we observe in our data?

To investigate the genuine or spurious nature of trimer and tetramer interactions, we systematically studied the interrelations between MC and FC entries. First, we define the dimer strength FC_r = Σ_i FC_ri of a region r as the sum of the strengths of all the dimers incident to it. Analogously, we introduced the (root-pinned) trimer strength MC_r = Σ_ij MC_ri,rj of a region r as the sum of the strengths of all the trimers of which r is the root. Conceptually, whereas FC_r measures the average coordinating influence that the region r exerts on its adjacent nodes, MC_r can be understood as quantifying the coordinating influence that r exerts on its incident links. As shown by Figure 5A, the correlations between dimer and trimer strengths of a region are weak and not significant, both at the global (black lines) and within each group (bundles of colored lines) levels, and for both within-zone and interzone trimers and dimers strengths. Of note, the average strength of between-zone trimers and dimers strengths had a larger variance across groups, hence the positively slanted shape of the global point cloud when confounding all groups, despite negative trends within each group. Although weak, within-subject correlations between FC_r and MC_r were negative, suggesting that some regions can be “meta-hubs” (Lombardo et al., 2020) but not “hubs,” that is, they can be the center of an entangled star subgraph of incident dimers, even if these dimers are individually weak and unable to systematically synchronize the fluctuations of adjacent nodes. Such meta-hubs could not have been identified through ordinary pairwise FC analyses only and manifest thus the existence of a real high-order multiregional coordination.

We then moved to consider how many trimers cannot be considered as a manifestation of underlying triangular motifs of dimers. We defined a trimer rooted in a region r to be genuine if MC_ri,rj > FC_ij, that is, if the observed trimer strength cannot be fully explained by a strong synchronization between the leaves. We then measured the observed fractions of genuine trimers. As shown by Figure 5B, substantial fractions of genuine trimers could be found for all trimer types: genuine fractions amounted to 32 ± 7% for within-zone trimers (root and both leaves in a same zone) and increased to 43 ± 13% for interzone trimers with leaves in two different zones, or 58 ± 9% for interzone trimers with the root in a different zone than the leaves. Especially for interzone trimers, many trimers could not be trivially explained by the existence of triangles of dimers.

Considering tetramers, we found larger redundancy with dimers. We defined the tetramer strength MC_ij = Σ_kl MC_ij,kl of a link (ij) as its total entanglement with other links. Figure S3A shows that a significant positive correlation existed between the dimer strength FC_ij of a link (ij) and its tetramer strength. That is, the stronger links were also the most entangled. Interestingly, several tetramers could still be considered genuine. We defined a tetramer genuine when MC_ij,kl > FC_ij, that is, when the two composing dimers were strongly correlated, despite (at least one of) the dimers being individually weak. Under this definition, Figure S3B shows that up to 55 ± 10% of tetramers composed of interzone dimers were genuine.

We conclude that, in general, the information conveyed by trimer and tetramer analyses is not completely redundant with the one conveyed by dimers, as many trimer and tetramer metrics cannot be explained solely in terms of dimers and thus express actual higher order correlations.

After defining various metrics to quantify the involvement of specific regions and links into pairwise and higher order interactions, as previously described, we then studied how dimer, trimer, and tetramer strengths varied across the four cohorts in our study.

First, we found that for both dimer and trimer interactions, the stronger effects were found considering interzone interactions. Figure 6A reports group differences for interzone dimers and Figure 6B for interzone trimers (mixed-zone or same-zone leaves are not treated separately). Results for within-zone dimers and trimers are shown in Figures S4A and S4B, respectively. In contrast to within-zone interactions, group-level comparisons for within-zone dimer and trimer interactions were not significant.

In general, when averaging over all brain regions (Figure 6A and B, left), general averages of dimer and trimer strengths progressively decreased from SNC, to NC, aMCI, and AD groups. This decrease, notably, was significant when comparing the two extreme SNC and AD groups. The effect was particularly strong for interzone trimer strengths (p = 0.005, Mann-Whitney U test, Bonferroni correction, for trimers), whose average value for the AD group not only decreased but changed its sign as it became negative. In contrast, within-zone trimer strengths remained strongly positive (Figure S4B). This means that, in the AD group, several regions are involved in a mixture of negative and positive trimer interactions. Positive interactions tend to synchronize the fluctuations of FC links, unlike negative interactions that tend to push them in an antiphase interaction. Furthermore, the mixture of positive and negative couplings results in a dynamic conflict scenario, known in the statistical mechanics as “frustration” (Vannimenus & Toulouse, 1977) and has been associated to disordered organization and a slowed-down relaxation to equilibrium (Mezard et al., 1988). The emergence of frustrated interzone trimer interactions is a strong qualitative discriminative marker of the AD group (see Discussion for possible interpretations of this finding).

The decrease of interzone trimer strengths and their switch to negativity in the AD group is confirmed also when focusing on individual brain regions, rather than the average (Figure 6B, right). Remarkably, strong decrease in trimer strengths were observed in regional subdivisions of the temporal pole and of the parahippocampal gyrus, along the hippocampus proper and amygdala. Some of these regions (entorhinal cortex in the parahimpocampal gyrus and the hippocampus), are among the first to be affected by neurofibrillary accumulation in AD pathology (Braak stages 1 and 2). In these same regions, we found a similar trend at the level of dimer strengths even when differences were not significant (Figure 6A, right). Of interest, the stronger effects at the level of dimer strengths were found in the cingulate gyrus, which are affected by early beta amyloid depositions and later on with neurofibrillary accumulation. Interestingly, the regions exhibiting the strongest effects at the level of trimers were not the ones with the strongest effects at the level of dimers (and vice versa; Figure 6A, right, vs. Figure 6B, right). The two analyses reveal thus complementary aspects of how pathology affects the spatiotemporal organization of functional interactions.

Lastly, we assessed differences on tetramer strengths across groups. In Figure 6 we show the average tetramer strengths for intrazone (Figure 6C) and interzone (Figure 6D) tetramers. In both cases, we observed a significant reduction of tetramer interactions from the SNC, to the NC, MCI, and AD groups. In the case of interzone tetramers, the drop in strength was large in the MCI group, with levels close to those in the AD group.

In summary, AD was associated with extensive reductions of not only dimer strengths, but more importantly, trimer and tetramer strengths. Furthermore, intergroup differences were salient when considering higher order trimer and tetramer compared to dimer interactions.

We have shown a large variety of changes associated with dFC across the cognitive spectrum from cognitively overperforming SNC subjects to AD. The rich set of complementary analysis approaches we deployed consistently converge toward a common message: AD is associated with a disordering of the rich spatiotemporal fluctuations that characterize healthy dFC.

It is worth noting that while BOLD activity misses many fast neuronal processes due to its slow sampling rate, what FC dynamics track are not neural level processes but variations of global brain state that can occur on much slower time scales. So dFC with a long TR accounts for variations of the way in which the repertoire of internal states is sampled, more than for variations of neural signals themselves. As a side note, these slow fluctuations are also what mean-field connectome-based whole-brain models are fit to reproduce via the stochastic sampling of their emergent repertoire of dynamic modes (Fousek et al., 2022; Hansen et al., 2015).

Our results showed that a pertinent description of dFC organization and its changes across groups can be formulated in terms of two anatomical zones segregating ventral from dorsal areas (Figure 6D). We found that the system spends less time in states with fluid Zone I dynamics and high global integration, visiting them more transiently, while it gets stuck on the contrary in less integrated states exhibiting Zone I hypersynchronization (Figure 2). At the dimer level, pairwise interactions between regions in different zones get more irregularly bursty (Figure 3). At the level of higher order trimers and tetramers, meta-connectivity analysis revealed a loss of coordination between the fluctuations of different sets of links, as quantified by dropping participation coefficients (Figure 4D). Trimer interactions between Zone I and Zone II, as well as tetramers, were weakened more distinctively than the interzone dimer interactions. Remarkably, regions in our limbic subnetwork for which conventional dimer analyses were not different between groups, showed a remarkably reduced involvement in trimer interactions between zones (Figure 6). Overall, these findings point together toward a “loss of structure” in dFC in parallel to the cognitive gradient across groups. This is in agreement with previous studies that showed a reduction of the complexity of spontaneous fluctuations of coordinated activity (Tait et al., 2020).

Nevertheless, even though being quite encouraging, a conclusive validation of our findings would require using larger cohorts, which preferably contains information on cortical thinning and PET scans of tau and Aβ depositions, to test whether their distributions correlate with the local network dynamics alterations we observe (thus establishing them as potential physiopathological causes of these changes) or not (advocating for alternative explanations; see later discussion). Similarly, our choice of regions and parcellations was arbitrary, generally based on the successful use of the same parcellation in previous modeling-based analyses of the same cohort (Zimmermann et al., 2018b). A better resolution fMRI from further cohorts would allow validating our results with finer and more extended parcellations, especially for the subcortical regions (Tian et al., 2020) that constitute the core of the limbic network on which we have focused.

Interestingly, our qualitative description emerges from radically different approaches to dFC parameterization: a state-based approach (the PBM method, by Thompson & Fransson, 2016); and a state-less approach (the random walk descriptions of dFC, by Battaglia et al., 2020, and Lombardo et al., 2020). The PBM method is firmly rooted in the developing field of temporal network theory (Holme & Saramäki, 2012). Temporal networks allow describing interregional communication as it unfolds in time, similarly to a call center, where operators can handle a multitude of brief first-contact calls at certain moments and dedicate extensive time to select customers at other times (Kovanen et al., 2013). Or to a primary school, where students interact in small groups during lectures and play in mixed larger groups in the playground during school breaks (Gemmetto et al., 2014). Eventually, even fluctuations between segregated or integrated states in brain systems at different scales (Pedreschi et al., 2020; Shine et al., 2016) give rise to network dynamics not dissimilar to these social systems. Note that our use of terms such as “burstiness” or “activation” (cf. Figure 3D and E) is also mediated from the jargon of temporal networks theory and should not be mistaken with the usual meaning of these terms in neuroscience, as they refer to FC link dynamics rather than to neuronal firing rates (exactly as we use the adjective “temporal” in the sense of “time dependent” and not in association with “temporal lobe”).

The dFC random walk approach (Arbabyazd et al., 2020; Battaglia et al., 2020; Lombardo et al., 2020; Petkoski et al., 2023) models rs dFC as a temporal network as well, but focuses on the variation from one network frame to the next, more than on the geometry of individual network frames. dFC is seen as a flow in network space and the nonrandomness of network reconfiguration was investigated via a time-to-time correlation approach known as meta-connectivity (Lombardo et al., 2020). In a dFC context in which the mode of coordination between regions is not frozen in time but changes smoothly, meta-connectivity reveals how the fluctuations of one or more regions modulate the degree of coordination between the fluctuations of other regions. In other words, meta-connectivity is an indicator of “many-body coordination.” Indeed, the terminology of dFC “dimers, trimers, tetramers” is reminiscent of perturbative diagrammatic expansions in statistical physics, such as the virial expansion (Landau & Lifshitz, 1980), in which clusters of increasingly large size account for progressively more elaborate and nonlinear patterns of many-body interactions. MC can thus be considered yet another form of high-order functional connectivity, adding up to a list of other approaches to track higher order coupling (Santoro et al., 2023; Torres et al., 2021) as hypergraph or homological methods (Battiston et al., 2020; Petri & Barrat, 2018; Sizemore et al., 2018), which have already identified synergistic aspects of human brain functioning (Luppi et al., 2022; Varley et al., 2023).

Unfortunately, both of the dFC methods implemented in this study provide results depending on specific parameter choices. For instance, concerning the state-less random walk approach, the selection of a window size remains ultimately arbitrary. The window size selected was short in contrast to other studies. However, our statistical analyses suggest that this window size results in similar discriminatory power as longer windows (Figure S2A). Furthermore, it is necessary to use short windows because the PBM method suggests that dwell times in consistent FC-state epochs are often short and thus dFC is intrinsically fast (Figure S2B). The need to track the covariance of fast FC fluctuations has inspired additional approaches analogous to MC, as edge-centric functional connectivity (eFC; Faskowitz et al., 2020). In this approach, covariance is estimated between individual events of instantaneous cofluctuation, without arbitrary windowing. However, we showed in Lombardo et al. (2020) that, despite the significant relation between MC and eFC, the use of a sliding window in the MC approach produces a smoothing effect that partially denoises the graph structure of interlink meta-connections, allowing a cleaner determination of modules and “meta-hub” nodes with large trimer strengths.

An additional aspect of the state-based PBM approach, is that it involves partially arbitrary steps as the choice of a number of states. The retrieved FC states depend on the extracting algorithm that depends on the brain parcellation and choice of regions of interest utilized. We found four states and increased dwell times in states with hyper-connectivity within Zone I. This finding of increased probability in AD of visiting hyper-connected states is in agreement with some state-based dFC studies (Gu et al., 2020), but in contrast with others (Fu et al., 2019; Schumacher et al., 2019), which instead find higher dwell times in disconnected states. Such discrepancies may arise because in the PBM method clustering of states is performed on activation patterns rather than on time-resolved functional networks. Our procedure has the advantage of showing that network dynamics is partially dissociated from node dynamics, with the possibility of hyper-connected FC modules arising both in presence of higher or lower activity of the nodes composing this module (Figure 2A). It may reduce the chance, however, of detecting extreme events along dFC or transient atypical network configurations that would be naturally assigned to separate clusters when directly clustering networks. Finally, the mentioned studies used reference parcellations with a larger number of regions or focused on more distributed network components, while here we particularly emphasize selected regions of interest, such as temporal and paralimbic cortices, known to develop epileptiform activity (Bakker et al., 2012; Cretin et al., 2016; Vossel et al., 2013). Thus, within the probed subsystem of interest, hypersynchrony may become particularly prominent and overexpressed (hence, the enhanced dwell time in hyper-connected FC states), a fact that has direct pathophysiological relevance.

Despite the arbitrary steps involved, both approaches independently provide sets of results with a high mutual consistency, making unlikely that our analyses reflect exclusively methods artifacts. Both methods confirm indeed that a dFC description in terms of two zones is pertinent, as the distinction between Zones I and II organizes the modular structure of both FC states in the state-based PBM approach (Figure 2A) and of the MC matrices in the state-free dFC random walk approach (Figure 5B and C). Furthermore, both methods confirm that the increased severity of cognitive decline across the four groups correlates with a reduced interzone coordination: more time spent in states with weaker integration (Figure 2B and C) and reduced interzone trimer strengths (Figure 6B). Such semantic agreement is remarkable, especially given the limitations of our approaches. Meta-connectivity analyses could be improved by seeking, beyond plain module detection, for a hierarchical community structure, that is often present in large networks (Jeub et al., 2018; Peixoto, 2014). State-based analyses could profit of better clustering approaches, as used by Rasero et al. (2018). However, while acknowledging these limitations, we found our four states and MC communities to be already highly interpretable, in term of the anatomical nature of the entangled links.

Particularly interesting is the fact that the weakening of interzone trimer interactions across the four groups decreases to such extent that some of these trimer switch from a positive to a negative value. As previously mentioned, the coexistence of negative and positive couplings in a graph or a hypergraph of interacting units is referred to in statistical physics as “frustration” (Toulouse, 1986), since it is associated with the emergence of conflicts preventing smooth relaxation to an equilibrium. To put these results in context, let us imagine that a dynamic FC link (a dimer FC_ij) is positively coupled to a second dimer FC_kl and negatively coupled to a third dimer FC_mn, and that the second and the third dimer simultaneously increase in strength (i.e., FC_kl and FC_mn get larger). Then the dynamics of FC_ij will “freeze” under the contrasting influence of the positive bias applied by FC_kl (pushing it to assumer stronger values), and the negative bias applied by FC_mn (pushing it to assume smaller values). Thus, the change of positive to a negative interzone influence—as the one signaled by the negative interzone trimer strengths of many limbic region within Zone I—gives rise to conflicts between the flows of Zone I and Zone II regions in AD patients, in contrast to control subjects where the fluctuations of the same regions are naturally synchronized.

In particular, in the context of cognition, Zone II regions such as the posterior cingulate cortex (pCC) have been postulated to play a regulatory role on the level of brain meta-stability, balancing “free-wheeling” internal cognition and focused outward attention (Leech et al., 2012; Leech & Sharp, 2014). In control groups, pCC has strong positive dimer coupling and moderately negative trimer coupling with regions in Zone I (Figure 6). This could allow the pCC to quickly coordinate with individual Zone I regions (and share information with them via direct positive FC dimers), while simultaneously “lowering the volume” of intrazone I communication (via pCC-rooted negative trimers with Zone I leaves). In AD subjects, this subtle equilibrium is lost, resulting potentially in perturbed integration of information within and between Zone I regions. Remarkably, pCC is also a key hub of the default mode network (Raichle et al., 2001), a system whose dFC had already been suggested as a biomarker in the conversion to AD (Jones et al., 2012; Puttaert et al., 2020).

Interestingly, our analyses on trimer strengths could detect intergroup differences within Zone I regions, for which the dimer analyses did not found significant differences. A possible explanation for the better sensitivity of trimer-based analyses could trivially be due to a larger sample size, as there were more possible trimers than dimers, resulting in similar average strengths but with a lower variance. However, another possibility could be that higher order interactions are readily affected by the pathology process earlier or at a higher degree than pairwise interactions. This fact is difficult to assess from our dataset, which is not longitudinal. Yet, this possibility is supported by our results showing that higher order trimer and tetramer terms convey in many cases genuinely new information, not redundant with dimer analyses. Indeed, even if we agree with other reports (Novelli & Razi, 2022) that dimer terms can sometimes explain trimer and tetramer term, we found in addition important trimer entanglement among otherwise individually weak dimers (Figure 5A) that lacked strong pairwise interactions between their dangling leaves (Figure 5B). Such genuine trimers cannot be explained by dimer motifs and describe thus a qualitatively different phenomenology, invisible to conventional FC analyses. Similar considerations apply to tetramers (Figure S3), which although generally weaker in strength than dimers and trimers, form an additional and pervasive background “medium” which also actively steer coordinated FC dimer fluctuations, with an overall influence degraded by the pathological process (Figure 5C and D). In the future, for an even better appreciation of pathology effects on higher order interactions, one may use methods that facilitates the generalization to arbitrarily high orders, even higher than the third or the fourth one, such as maximum entropy fitting (Ezaki et al., 2018; Savin & Tkačik, 2017) or other information theory approaches (Rosas et al., 2019).

Another question is what the mechanistic origin could be of the observed spatiotemporal complexity of dFC (and of its alterations). Previous studies have shown that structured dFC may emerge as an effect of global brain network dynamics to be tuned at a slightly subcritical working point (Arbabyazd et al., 2020; Glomb et al., 2017; Hansen et al., 2015), or as a consequence of cascades of neuronal activations (Rabuffo et al., 2021) that occur due to the flow on the manifold created by the symmetry breaking of the connectome (Fousek et al., 2022). However, these studies did not use very precise criteria when referring to their capacity to render dFC. In the future, the statistical descriptors of dFC alterations that we introduce here, such as regional spectra of trimer and tetramer strengths, may be used as more detailed fitting targets for the tuning of mean-field models aiming at explaining the circuit mechanisms for the emergence of higher order interactions. Such models, once fitted, may also allow reverse engineering the physiological changes that are responsible for the degradation of spatiotemporal dFC complexity along the SNC-to-AD spectrum.

It is likely that the dFC alterations we observe between groups are caused at least in part by underlying biological causes of AD, as the aggregation of misfolded proteins that cause cell death and atrophy (Soto & Pritzkow, 2018). However, not all the symptoms can be explained by these mechanisms. Among them, the existence of symptom severity fluctuating across hours in a way not accountable for sudden variations of amyloid load (Palop et al., 2006) or, yet, the phenomenon of cognitive reserve where subjects with virtually identical or even higher amount of amyloid load than others can maintain a very efficient cognition, (cf. Snowdon [2003] for the famous “Nun Study” or Rentz et al. [2010] for a review of other studies with similar conclusions). These findings suggest that neurodegeneration may coexist with compensations of unspecified nature that allow “cognitive software” to operate properly despite “hardware damage” (see, e.g., Petkoski et al. [2023] for examples of dynamic compensation in healthy aging, or Courtiol et al. [2020] for a similar phenomenon in epilepsy). Here, we propose the hypothesis that preserved dFC complexity may act as a possible form of cognitive reserve. We stress once again that, to check the soundness of this hypothesis, future analysis should rely on richer datasets that contain PET scans of tau and Aβ depositions, and possibly even a mechanistic model (Stefanovski et al., 2019, 2021) for their impact to the neuronal activity.

Ultimately, the degradation of dFC organizational complexity that we here described may not only correlate with cognitive decline but also, eventually, contribute to cause it. Indeed, a dFC with a complex organization could be the hallmark of brain dynamics implementing “healthy” cognitive processing. Computation can emerge from collective dynamics as long as this dynamics is sufficiently complex, that is, neither too ordered nor too random (Crutchfield, 2012; Crutchfield & Mitchell, 1995). More fundamentally, the existence of alternative information processing states—transient FC networks?—and of nonrandom transitions between these states—structured and complex dFC switching?—are two necessary conditions for whatever information processing system to perform computation (Turing, 1937). A speculative hypothesis is thus that the complexity of neural dynamics—and, more specifically the complexity of ongoing dFC, which is a measurable shadow of hidden neural processes—is an instrumental resource for cognitive information processing. Cognitive deficits in pathology could arise just in virtue of this resource becoming scarcer, because of less structured and more random dynamics. This phenomenon has been speculatively observed in hippocampal neuronal assembly dynamics in epilepsy (Clawson et al., 2023). In this line of thinking, preserved dFC complexity would act as a “dynamic reserve” allowing the implementation of elaborate neural computations (or “software patches”) to compensate for progressing neurodegeneration. Analogously, enhanced dynamic complexity could be the substrate for the superior cognitive performance achieved by subjects in the SNC group with respect to NC subjects. A more direct exploration of the link between dFC complexity and cognitive processing in the healthy and pathological brain will be needed to inquire into this suggestive hypothesis.

The study included 73 subjects between 70 and 90 years of age from the fourth wave of the Sydney Memory and Ageing Study (Sachdev et al., 2010; Tsang et al., 2013). The use of the database was approved by the Human Research Ethics Committee of the University Texas at Dallas. For detailed descriptive summaries on neuropsychological assessments for AD and amnesic aMCI, we refer the reader to Zimmermann et al. (2018b).

A specificity of our approach is the stratification of healthy controls with an additional “supernormal” category, putting our focus not only on mechanisms of disease but also on mechanisms of “health” based on cognitive performance. Results from 12 neuropsychological tests were combined in the following cognitive domains: attention/processing speed, memory, language, visuospatial ability, and executive function. In brief (Mapstone et al., 2017), we classified cognitive membership for each subject based on the composite Z scores as supernormal controls (SNC) or normal controls (NC). For this, the supernormal (SNC) group was defined as Z_mem > 1.35 SD (∼90th percentile) and Z_cog > 0.7 SD. The normal control participants are conservatively defined with Z_mem ± 0.7 SD (∼15th percentile to 85th percentile) of the cohort median. The classification of subjects as AD and aMCI described in Zimmermann et al. (2018b) was done by consensus and included the following: the amnesic MCI group was described by a cognitive decline at least in the memory domain (Z_mem and/or Z_cog < 1.5 SD below normative values), paired to subjective complaint of cognitive deficit and without deficits in activities of daily living (ADL). The AD group in presence of a diagnosis of Alzheimer’s disease according to DSM-IV criteria (American Psychiatric Association, 2000) assessed by a clinical expert panel that included significant cognitive decline in several cognitive domains in addition to significant decrease in ADLs (American Psychiatric Association, 2000; Zimmermann et al., 2018b).

Details about resting-state functional MRI acquisition and preprocessing can be found in Zimmermann et al. (2018b). We briefly mention, as relevant here that during the fMRI acquisition, participants were instructed to lie quietly in the scanner with their eyes closed. The TR used for the T2* weighted EPI sequence of time-resolved BOLD imaging was 2,000 ms. The acquisition time was of ∼7 min. Data from all MRI modalities was preprocessed using FSL and QA followed Smith et al. (2004). Subjects were removed if any of their scan acquisitions contained excessive artifacts including slice dropouts on the diffusion images (defined by zebra-like blurring or complete dropout; Pannek et al., 2012), the presence of orbitofrontal EPI signal dropout (Weiskopf et al., 2007), excessive motion on T1-images (i.e., ringing), or severe geometric warping. For details of additional fMRI preprocessing details (slice-timing correction, realignment and coregistration, linear detrending, head motion regression, probabilistic segmentation, spatial smoothing, etc.) please refer to Perry et al. (2017).

For structural and functional parcellation, the AAL atlas was used focused on 16 limbic regions (see Figure 6D) associated with early degeneration in AD according to Braak and Braak staging, as we did before (Zimmermann et al., 2018a). The regions of interest included: cingulate cortices (anterior, medial and posterior), parahippocampal gyrus (including entorhinal cortex), hippocampus proper, amygdala, and temporal pole (superior and middle). In this study, as pertinent given the spatial organization retrieved in many of the analysis results, we categorize regions as belonging: either to Zone I, including ventral regions (superior and medial portion of the temporal pole, parahippocampal gyrus, hippocampus proper, and amygdala in both hemispheres); or to Zone II, which included the six cingulate cortical regions (posterior, medial, and anterior) in both hemispheres (Figure 6D). This subdivision in two separate zones allowed us the categorization of network links from dimers to higher order arrangements (trimers, tetramers) determining “within zone” or “between-zone” interactions based on the relative zone membership of the different nodes involved. We remark that the delimitations of Zone I and Zone II are inspired from data-driven considerations (the spatial organization of FC-state centroids in Figure 2 and MC modules in Figure 4) rather than from a priori subdivisions.

In this study, we applied two complementary dynamic functional connectivity (dFC) approaches to investigate nonstationarity of BOLD signals and capture the recurring, time-varying, functional patterns. The first one was the so-called point-based method (PBM) introduced by Thompson and Fransson (2016), referred here as state-based dFC. This method assumes the existence of a small set of possible discrete FC configurations.

In this approach, BOLD signals of each subject were concatenated along the temporal dimension and transformed to z scores using Fisher’s z-transformation to stabilize variance prior to further analysis. Following Thompson and Fransson (2016), we applied a k-means clustering algorithm on the concatenated time series (Lloyd, 1982) to determine states based on global activity patterns (best partition out of 100 repetitions, max iterations 100). The optimal number of 4 clusters (k = 4) was validated based on detecting an elbow in the variation of the distortion score as a function of changing number of clusters k (Figure S1B). Based on the collections of activity patterns at times assigned to each of the states, we computed Pearson correlation matrices, yielding k state-specific FC matrix FC^(λ) (λ =1 … 4). A state was hence characterized by the centroid activation pattern of time frames within the state cluster and by its state-specific FC matrix (see Figure 2A and Figure S1A). To characterize the spatial properties of state-specific FC, we then used a graph-theoretical approach and measured global and local efficiencies (Achard & Bullmore, 2007; Latora & Marchiori, 2001) of the four FC^(λ) networks (Figure 2B) using the Brain Connectivity Toolbox (Rubinov & Sporns, 2010).

In a second approach, we assumed that FC networks are continually morphing in time, without priors on the existence of discrete state-switching events, following Battaglia et al. (2020), that conceptualized the evolution of FC as a stochastic walk in the high-dimensional space of possible network configurations. This stochastic walk, however, is not trivial, as different interregional links covary according to a specific higher order correlation structure called meta-connectivity (Lombardo et al., 2020). State-free and smoothly varying dFC temporal networks were extracted using a sliding-window approach, adopting the random walks and meta-connectivity approaches (Battaglia et al., 2020; Lombardo et al., 2020; Petkoski et al., 2023) released within the dFCwalk toolbox (Arbabyazd et al., 2020).

Following and based on the correlation matrix between “dimers” (dynamic FC links between two regions i and j), the entries MC_{ij, kl} of the MC matrix are either computed based on the dynamics of four regions involved in the links (ij) and (kl), or at least three regions, when the two considered dimers share a common vertex (e.g., i = k). MC can thus be seen as a collation of higher order interactions within the system, involving more than “two parts” (tetramers or trimers). In the case of a trimer, the region on which the two dimers converge to a “root” region, and the other two regions are the “leaves” of the trimer. In the case of a tetramer, each of the two nonincident dimers are called a “base.”

We used a graph-theory approach to quantify the communities of MC matrices. MC for all subjects were constructed and then averaged for each of the four subject’s groups (Figure 3B). To detect the modular structures of MC, we used the community Louvain algorithm (Rubinov & Sporns, 2011). We used a parameter Γ = 1.4, determined heuristically to yield a modular partition naturally interpretable in anatomical terms. To quantify the modularity changes across the groups, we computed the index of modularity (Q*) as measure of degree of intramodule connectivity. Since MC is a signed matrix, we applied disproportionate scaling to the positive and negative values of modularity indices to consider a lower contribution of negative meta-link weights to the index of modularity (Rubinov & Sporns, 2011). To quantify the degree of intermodular connectivity of group-averaged MCs, we computed the Participation coefficient of each dFC dimer node following (Guimerà & Amaral, 2005). This metric can be computed exactly as for an ordinary graph keeping in mind that FC links and meta-links among them are, respectively, nodes and links in the MC graph. The participation coefficient is close to one when meta-links of a link are distributed uniformly, therefore, integrated across MC modules and it is zero when all the meta-links of a link are segregated within its own MC module.

MC describes largely delocalized interactions but, for enhanced interpretability, it is important to describe the overall contribution of individual regions to the different higher order interactions. Hence, we defined various indices of meta-strength.

A trimer is defined between zones or within zones depending on the zones to which its leaves belong. If all leaves are in the same zone (independently from where the root is), then the trimer is considered within zone, otherwise it is considered between zones. For tetramers, we distinguished tetramers with base within a zone (if both bases are within zone dimers) or base between zones (if both bases are between zones). There are more combinatorial cases for tetramers that were ignored in this study for simplicity.

The authors wish to thank Leonardo Novelli and Giovanni Petri for constructive discussions about our results and methods.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00332.

Lucas Arbabyazd: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Software; Visualization; Writing – original draft; Writing – review & editing. Spase Petkoski: Conceptualization; Formal analysis; Methodology; Writing – original draft; Writing – review & editing. Michael Breakspear: Conceptualization; Data curation; Methodology; Writing – original draft; Writing – review & editing. Ana Solodkin: Conceptualization; Data curation; Investigation; Methodology; Writing – original draft; Writing – review & editing. Demian Battaglia: Conceptualization; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Resources; Software; Supervision; Validation; Writing – original draft; Writing – review & editing. Viktor Jirsa: Conceptualization; Data curation; Funding acquisition; Investigation; Resources; Supervision; Validation; Writing – original draft; Writing – review & editing.

Viktor Jirsa, Horizon 2020 Framework Programme (https://dx.doi.org/10.13039/100010661), Award ID: 945539. Viktor Jirsa, Horizon 2020 Framework Programme (https://dx.doi.org/10.13039/100010661), Award ID: 826421.

Aβ and Tau:

Symptoms of AD correlate with the accumulation of plaques and tangles formed by these two proteins (physiopathological causes).

Meta-connectivity:

Correlation between the fluctuations of functional links (which are themselves correlations between nodes’ fluctuating activity), i.e., a connectivity between connections.

S-graphlet:

State-specific functional connectivity estimated over time epochs belonging to a specific activity cluster, obtained via unsupervised clustering.

Limbic system and posterior cingulate cortex:

Brain regions affected by AD. The first include hippocampus (memory function), the second is a node of default mode network.

T-graphlet:

A temporal network in which every time frame is the S-graphlet of the state visited at each considered time instant.

Link burstiness, mean and total activation, memory:

How much link activation events last and how (ir)regular they are is captured by β, μ, τ, and λ coefficients.

Trimer strength:

How much a node is coordinating the fluctuations of incident links. Estimated from the (correlated) activity of triplets of regions.

Genuine trimers:

A trimer meta-connection that cannot be explained by the correlations that a root regions establishes independently with the two leaves.

Braak stages:

Neurodegeneration in AD (and other neurological diseases) affects progressively regions in precisely ordered sets, characterized by Heiko Braak and colleagues.