# Dynamic Functional Connectivity between order and randomness and its evolution across the human adult lifespan

Demian Battaglia, Thomas Boudou, Enrique C.A. Hansen, Diego Lombardo, Sabrina Chettouf, Andreas Daffertshofer, Anthony Mcintosh, Joelle Zimmermann, Petra Ritter, Viktor Jirsa

<figure><img src="image_1.png" /></figure>

## To cite this version:

Demian Battaglia, Thomas Boudou, Enrique C.A. Hansen, Diego Lombardo, Sabrina Chettouf, et al. Dynamic Functional Connectivity between order and randomness and its evolution across the human adult lifespan. NeuroImage, 2020, 222, pp.117156. 10.1016/j.neuroimage.2020.117156. hal-03092853

HAL Id: hal-03092853

https://hal.science/hal-03092853

Submitted on 5 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Dynamic Functional Connectivity between order and randomness and its evolution across the human adult lifespan

<figure><img src="image_4.png" /></figure>

Check for
updates

Demian Battaglia<sup>a, *</sup>, Thomas Boudou<sup>a, b</sup>, Enrique C.A. Hansen<sup>a, c</sup>, Diego Lombardo<sup>a</sup>, Sabrina Chettouf<sup>d, e, f</sup>, Andreas Daffertshofer<sup>f</sup>, Anthony R. McIntosh<sup>g</sup>, Joelle Zimmermann<sup>d, g</sup>, Petra Ritter<sup>d, e, 1</sup>, Viktor Jirsa<sup>a, 1</sup>

## ARTICLE INFO

### Keywords:

Dynamic functional connectivity

Resting-state

fMRI

Aging | anomalous diffusion

## ABSTRACT

Functional Connectivity (FC) during resting-state or task conditions is not static but inherently dynamic. Yet, there is no consensus on whether fluctuations in FC may resemble isolated transitions between discrete FC states rather than continuous changes. This quarrel hampers advancing the study of dynamic FC. This is unfortunate as the structure of fluctuations in FC can certainly provide more information about developmental changes, aging, and progression of pathologies. We merge the two perspectives and consider dynamic FC as an ongoing network reconfiguration, including a stochastic exploration of the space of possible steady FC states. The statistical properties of this random walk deviate both from a purely “order-driven” dynamics, in which the mean FC is preserved, and from a purely “randomness-driven” scenario, in which fluctuations of FC remain uncorrelated over time. Instead, dynamic FC has a complex structure endowed with long-range sequential correlations that give rise to transient slowing and acceleration epochs in the continuous flow of reconfiguration. Our analysis for fMRI data in healthy elderly revealed that dynamic FC tends to slow down and becomes less complex as well as more random with increasing age. These effects appear to be strongly associated with age-related changes in behavioural and cognitive performance.

## 1. Introduction

Recent studies emphasized the structured temporal variability of resting-state (rs) and task Functional Connectivity (FC; Tagliazucchi et al., 2012; Allen et al., 2012; Liu and Duyn, 2013; Hutchison et al., 2013; Chen et al., 2015; Preti et al., 2017; Gonzalez-Castillo and Bandettini, 2018), whose study is the defining focus of a research direction recently designated as “chronnectomics” (Calhoun et al., 2014). If rs-FC

is dynamic, a wealth of information may be lost by averaging over long imaging sessions, and averaged temporal variability might be mistaken as inter-subject variability. Temporal FC variability—which we will refer to as *dynamic FC* (dFC)— may carry an inherent meaning. It has been suggested to manifest ongoing cognition at rest (Gonzalez-Castillo et al., 2015) with an immediate impact on cognitive performance (Bassett et al., 2011; Braun et al., 2015; Shine et al., 2016; Cohen, 2018) and attentional or awareness levels (Kucyi et al., 2017; Cavanna et al., 2018; Lim et al.,

Abbreviations: rs, resting-state; RSN, resting-state network; fMRI, functional magnetic resonance imaging; BOLD, blood oxygen level dependent; SC, structural connectivity; FC, functional connectivity; dFC, dynamic functional connectivity; DFA, detrended fluctuation analysis; MoCA, Montreal Cognitive Assessment; SO, spectral overlap.

* Corresponding author.

E-mail addresses: demian.battaglia@univ-amu.fr (D. Battaglia), thomas.boudou@ensta-paristech.fr (T. Boudou), ehansen@biologie.ens.fr (E.C.A. Hansen), diego.lombardo@univ-amu.fr (D. Lombardo), sabrina.chettouf@charite.de (S. Chettouf), a.daffertshofer@vu.nl (A. Daffertshofer), rmcintosh@research.baycrest.org (A.R. McIntosh), jzimmermann@research.baycrest.org (J. Zimmermann), petra.ritter@charite.de (P. Ritter), viktor.jirsa@univ-amu.fr (V. Jirsa).

¹ Shared last authorship.

<footnote><sup>a</sup> Université Aix-Marseille, INSERM UMR 1106, Institut de Neurosciences des Systèmes, F-13005, Marseille, France</footnote>

<footnote><sup>b</sup> ENSTA ParisTech, F-91762, Palaiseau, France</footnote>

<footnote><sup>c</sup> Institut de biologie de l'Ecole normale supérieure (IBENS), École Normale Supérieure, CNRS, INSERM, PSL Université Paris, F-75005, Paris, France</footnote>

<footnote><sup>d</sup> Brain Simulation Section, Department of Neurology, Charité Universitätsmedizin and Berlin Institute of Health, D-10117, Berlin, Germany</footnote>

<footnote><sup>e</sup> Bernstein Center for Computational Neuroscience, D-10117, Berlin, Germany</footnote>

<footnote><sup>f</sup> Faculty of Behavioural and Movement Sciences, Vrije Universiteit Amsterdam, 1081 BT, Amsterdam, the Netherlands</footnote>

<footnote><sup>g</sup> Rotman Research Institute, Baycrest Centre, Toronto, Ontario, M6A 2E1, Canada</footnote>
2018). It may reflect the sampling of an internal repertoire of alternative dynamical states (Hansen et al., 2015; Golos et al., 2015; Cabral et al., 2017a; Glomb et al., 2017). From a biomarking perspective, pathological conditions, such as Alzheimer's dementia, schizophrenia and other mental disorders (Jones et al., 2012; Damaraju et al., 2014; Braun et al., 2018), or differences in general attributes like gender (Yaesoubi et al., 2015) or development and aging (Hutchison and Morton, 2015; Qin et al., 2015; Davison et al., 2016; Schlesinger et al., 2016; Chen et al., 2017; Viviano et al., 2017), may alter dFC more than they affect time-averaged FC.

The growth in the number of dFC studies, based on both fMRI and electrophysiological signals, has been paralleled by an increasing number of possible technical approaches to estimate dFC (Preti et al., 2017). A non-exhaustive list ranges from sliding window approaches (Allen et al., 2012), to statistical modelling of signals (Lindquist et al., 2014) and state transitions (Baker et al., 2014; Vidaurre et al., 2016; Cabral et al., 2017b), temporal network approaches (Thompson and Fransson, 2016), or the study of coactivation patterns (CAPs: Chen et al., 2015; Matsui et al., 2016). However, several concerns have been raised on whether dFC reflects genuine neural network dynamics or rather artefactual fluctuations, linked, e.g., to head motion (Laumann et al., 2016) or signal processing aspects (Leonardi and Van de Ville, 2015). There are also several statistical concerns about whether resting-state FC is really non-stationary (Zaleski et al., 2014; Hindriks et al., 2016) or whether discrete connectivity states exist that might be reliably extracted (Shakil et al., 2016; Liégeois et al., 2017). In fact, while “FC clusters” can always be extracted using ad hoc algorithmic methods, as of yet that is not evident that such clusters correspond to well-defined, distinct attractor states (Zaleski and Breakspear, 2015).

Here, we introduce yet another way to look at dFC, which, we believe, circumvents some of the concerns mentioned above on the difficulty of assessing the actual non-stationarity of FC fluctuations in time. We do not attempt segmenting dFC in a sequence of sharp switching transitions between FC states but instead describe it as a smooth flow across continually morphing connectivity configurations. Conventional analyses of static FC emphasize the spatial structure of FC networks discarding most information about time. We adopt the opposite approach: de-emphasizing space and collapsing FC networks to a “point” in the space of possible FC network realizations. And, we interpret its erratic evolution as a random walk, a stochastic exploration of a high-dimensional space. An iconic example of random walk is the one of a grain of dust floating in the air, made visible through a beam of shining light. By sampling the position of this grain at different moments in time, one may observe that the displacements of the grain over a fixed observation interval are not constant in time. On the contrary, the distances travelled over different time intervals (or, equivalently, the speeds of motion) give rise to a distribution, which will be (close to) a Gaussian with a well-defined peak. Thus, the mean speed of the moving particle constitutes a first metric for the quantitative characterization of its random walk. As a second step, one can consider the statistical properties of the fluctuations of random walk speed around this mean. Not all the random walks are identical and exploration paths with different shapes (i.e., fractal geometries) might be generated, depending on the degree of correlation between the step lengths travelled at consecutive times (Mandelbrot and Van Ness, 1968; Mandelbrot, 1983). In some random walks –such as the one of a grain of dust, or of a particle suspended in water, originally studied by Robert Brown (1828)– the speed of consecutive exploration steps are uncorrelated (“memoryless”), giving rise to unstructured uniformly space-filling paths. In some other cases –such as the spread of infectious diseases in contemporary times (Brockmann et al., 2006) or foraging for resources in an ecosystem (Viswanathan et al., 1999)–, consecutive random walk steps have correlated speeds (“memory”), giving rise to paths with a characteristic alternation between exploration of localized clusters and longer jumps. Methods like detrended fluctuation analysis (DFA, Peng et al., 1993) serve to characterize the geometry of empirically observed random walk paths, so to

estimate their degree of memory and their deviation from Gaussian uncorrelated walks. Such approaches have been successfully applied to detect anomalous scaling properties of fluctuations in neural datasets at different spatial and temporal scales (He, 2014), associated to critical behavior in both rest and task conditions (Linkenkaer-Hansen et al. 2001; Chialvo, 2010; Palva et al., 2013). Here, with a picture in mind of resting state dFC as a random walk, we will study its typical speed and the shape of exploration paths in FC space.

We extend previous results that indicated how ongoing FC fluctuations implement a complex random walk, endowed with non-trivial statistical properties. The features of this dFC random walk appear to be intermediate between two possible trivial null hypotheses: a first scenario of “order” in which stationarity is strictly imposed, and a second scenario of “randomness” in which the time ordering of the observed time-resolved FC matrices is shuffled to destroy any long-range sequential correlations. Using DFA analysis and statistical comparisons with suitable surrogate dFC streams, we find that dFC deviates from both “order” and “randomness”, being thus “complex” (Crutchfield, 2011).

We also investigate how random walk properties of rest and task dFC may be modified over the healthy human adult lifespan. As we age, our brain undergoes characteristic structural and functional changes, with a tendency toward increased structural disconnection (Salat, 2011), disruptions in rs-FC (Andrews-Hanna et al., 2007; Betzel et al., 2014) and modified structural-to-functional connectivity inter-relations (Zimmermann et al., 2016). Analogously, changes in dFC have been reported at the level of the temporal stability of FC network modules (Davison et al., 2016; Schlesinger et al., 2016), general or specific network variability (Qin et al., 2015; Chen et al., 2017), “FC state” occupancy (Hutchison and Morton, 2015; Viviano et al., 2017), and complexity of phase synchrony (Nobukawa et al., 2019). We complement these previous findings and show that dFC random walks may occur at an increasingly reduced speed and complexity with age, slowing down and becoming increasingly more “random”. These reductions in dFC speed and complexity correlate with the level of general cognitive and behavioural performance, as probed by both standard clinical assessments of cognitive impairments (Nasreddine et al., 2005) and a simple visuomotor coordination task (Houweling et al., 2008). For a refined analysis of dFC random walk alterations along with more subtle and task-specific cognitive decline, see a related paper by Lombardo et al. (2020). The fact that slowing down and complexity loss in dFC are associated with degraded general performance can be linked to prominent theories of cognitive aging, speculatively establishing alterations of dFC random walk properties as novel imaging correlates of processing speed reduction (Salthouse, 1996; Finkel et al., 2007) and de-differentiation (Baltes, 1980; Sleimen-Malkoun et al., 2014).

## 2. Materials and methods

### 2.1. Experimental subjects

Overall $N = 85$ healthy adult subjects ($N = 53$ females, $N = 32$ males) were recruited at Charité - Universitätsmedizin Berlin to voluntarily participate in rs-fMRI and DSI scans and, for a subset of them, also in a visuomotor 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. These two subsets of subjects were scanned in the framework of initially independent studies (using the same MR-scanner) and not initially intended to study dFC, but could be combined at least for subsets of the analyses, to increase sample size whenever possible.

All subjects had no self-reported neurological, psychiatric, or somatic conditions. Subjects in the ‘rs + tasks’ subset, in addition to resting state analyses, were tested behaviourally with a visuomotor task (see later). Furthermore, for 20 of them we also assessed their general cognitive function with the Montreal Cognitive Assessment (MoCA) (Nasreddine et al., 2005). For all the rs analyses of Figs. 1–5 and S1–4, 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 with lower than median age (overall N = 42, 18-42 yrs, median age = 24 yrs), and an 'Older group' composed of subjects with larger than median age (overall N = 42, 47-80 yrs, median age = 63 yrs).

In addition to general exclusion criteria for participation in an MRI experiment, we excluded subjects with a self-reported musical background, as musical training may affect the performance of rhythmic visuomotor tasks. 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 before data collection. The studies were performed in accordance with the local medical ethics committee protocol at the Charité Hospital (Berlin, Germany).

## 2.2. MRI acquisition

Magnetic resonance imaging (MRI) acquisition was performed on a 3T Siemens Tim Trio scanner. Every subject was scanned in a session that included a localizer sequence (3, 8 mm slices, repetition time [TR] = 20 ms, echo time [TE] = 5 ms, voxel size = 1.9 × 1.5 × 8.0 mm, flip angle [FA] = 40°, field of view [FoV] = 280 mm, 192 mm matrix), a T1-weighted high-resolution image (MPRAGE sequence, 192, 1 mm sagittal slices, voxel size 1 × 1 × 1 mm, TR = 1940 ms, TE = 2.52 ms, FA = 9°, FoV = 256 mm, 256 mm matrix), a T2 weighted image (2:16 min, 48, 3 mm slices, voxel size 0.9 × 0.9 × 3 mm, TR = 2640 ms, TE1 = 11 ms, TE2 = 89 ms, FoV 220 mm, 256 mm matrix), followed by diffusion weighted imaging (61, 2 mm transversal slices, voxel size = 2.3 × 2.3 × 2.3 mm, TR = 7500, TE = 86 ms, FoV 220 mm, 96 mm matrix). Subjects were then removed from the scanner to have their EEG cap put on, and then simultaneous fMRI-EEG images were acquired in a single run (BOLD T2* weighted, 32, 3 mm transversal slices, voxel size = 3 × 3 × 3 mm, TR = 1940 ms, TE = 30 ms, FA = 78°, FoV = 192 mm, 64 mm matrix). Five dummy scans were automatically discarded by the Siemens scanner.

During resting-state scans, subjects were to remain awake and reduce head movement. Head cushions served 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, 20 min of uninterrupted rs scan were performed. For the 'rs + task' subset, 5 min of rs were collected before 20 min of task acquisition (see later), and then further 5 min after the task. We ignored differences between the two rs blocks by concatenating them prior to subsequent analysis.

## 2.3. fMRI processing

fMRI data were pre-processed following Schirner et al. (2015). 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. As an additional correction measure, we further regressed out six FSL head motion parameters from the measured 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 (Desikan et al., 2006). See Table S1 for the regions of interest. We neither performed a slice-timing correction, smoothing, normalization of BOLD intensities to a mean, nor global regression.

## 2.4. Visuomotor coordination task

The N = 36 subjects in the 'rs + tasks' subset performed a visuomotor coordination task while in the scanner (Chettouf et al., 2020). The task

followed an unimanual paradigm, which was adapted from a bimanual paradigm introduced in (Houweling et al., 2008). 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 a 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 to achieve synchrony.

Behavioural measures were collected (one performance score per trial) based on the frequency locking of the two rotating circles. If the two disks rotated perfectly in-synchrony (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 Spectral overlap (SO) between the power spectra of the two forces, $P_x$ and $P_y$, as described in detail in (Daffertshofer et al., 2000):

$$SO = \frac{2 \int P_x(f) P_y(\rho f) df}{\int [P_x^2(f) + P_y^2(f)] df}$$

with $\rho = 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; here we focused 'only' on the average performance over trials, ignoring all learning effects.

## 2.5. 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; visuospatial 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 (Nasreddine et al., 2005).

## 2.6. Control resting state fMRI dataset

To verify that the main results of our analyses did not apply just to our specific dataset, we also repeated some of the analyses on an independent control dataset, selected from rs-fMRI data released as part of the Human Connectome Project (HCP), WU-Minn Consortium. We used the same selected subjects used for statistical benchmarking of functional network analyses in Termenon et al. (2016). This sample includes 99 young healthy adults from 20 to 35 years old (54 females). Each subject underwent two rs-fMRI acquisitions on different days and here, not being interested in test-retesting issue, we used only data from the first day sessions. For these sessions, TR = 720 ms and resting state scan duration was of 14 min and 24 s. More details on data acquisition and pre-processing for this control dataset can be found on Termenon et al.
(2016).

## 2.7. Extraction of time-dependent functional connectivity and dFC matrices

In brief, we estimated the sequence of time-dependent Functional Connectivity matrices $FC(t)$ –or dFC stream– by sliding a temporal window of fixed duration $\tau$ (cf. Allen et al., 2012) and by evaluating zero-lag Pearson correlations between resting-state BOLD time series from different brain regions i and j:

$$FC_{ij}corr[BOLD_{i}(t'), BOLD_{j}(t')] \text{ over the interval } t - \frac{\tau}{2} \le i \le t + \frac{\tau}{2}$$

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).

To evaluate the dFC matrices of Figs. 2 and 4 and S2 we introduced a notion of similarity between $FC(t)$ matrices following (Hansen et al., 2015), based on the Pearson correlation between the entries of their upper-triangular parts:

$$dFC(t_1, t_2) = corr[UpperTri(FC(t_1)), UpperTri(FC(t_2))]$$

The dFC matrices thus depend on the window-size $\tau$ adopted when

evaluating the dFC stream. To perform the two-dimensional projections of the sequence of $FC(t)$ matrices in Fig. 2, the vectors $UpperTri(FC(t))$ served as input features into a t-Stochastic Neighborhood Embedding algorithm (exact method) as described by Hinton and Van der Maaten (2008), with default perplexity = 30 and exaggeration = 4 parameters in the employed MATLAB® (MathWorks R2017b) implementation. To compare projections of actual empirical dFC and dFC evaluated for surrogate data (see below) we determined a common projection based on a unified training set combining empirical and surrogate $FC(t)$ matrices.

## 2.8. Surrogate dFC models

We compared dFC streams between actual empirical data and two different types of surrogate, to probe deviations from two alternative null hypotheses.

We constructed phase-randomized surrogates following Hindriks et al. (2016). In brief, the multi-variate time-series of BOLD signals were first Fourier transformed, extracting time-dependent amplitude and phases at different frequencies for every region. Subsequently, the Fourier phases were drawn from a uniform distribution $\sim(0, 2\pi)$ before applying an inverse Fourier transform yielding surrogate time-series (Kantz and Schreiber, 2004; chap. 7.1.2) with unaltered power spectra. Note, however, that additional care needs to be taken to preserve the overall

**Fig. 1.** Streams of dynamic Functional Connectivity (dFC). (A) We adopt a sliding-window approach to estimate temporal changes of dynamic Functional Connectivity from fMRI data (resting-state or task). We call the resulting smooth sequence of time-resolved $FC(t)$ matrices a dFC stream. We then measure at any time $t$ the dFC speeds $V_{dFC}, \tau(t)$ as the variation of $FC(t)$ observed between a time $t$ and a time $t + \tau$, where $\tau$ is the size of the window used to estimate $FC(t)$ matrices. (B) The degree of similarity (inter-matrix correlation) between $FC(t)$ networks observed at different times can then represented into a recurrence matrix, or dFC matrix, whose entry $dFC(t_1, t_2)$ reports the correlation distance between the functional networks $FC(t_1)$ and $FC(t_2)$ estimated respectively at times $t_1$ or $t_2$. The block structure of dFC matrices reflects the inhomogeneous speed of variability of $FC(t)$ along the flow of the dFC stream. The quantification of dFC speeds allows answering different questions about the statistical properties of dFC streams, summarized by two graphical cartoons: first, (C) whether $FC(t)$ matrices tend to cluster into FC states dissimilar between them, but internally similar (stationarity or non-stationarity of dFC); second, (D) whether there are sequential correlations in the dFC stream (e.g. short “jumps” followed by other short “jumps” with high probability), which could be possible even in the case of stationarity. In panel (D), exclamation marks highlight that “dFC knots” –epochs of sequentially correlated short jumps– are not clusters since they can have an internal variance larger than the variance between time-resolved $FC(t)$ networks separated by “dFC leaps”.

<figure><img src="image_5.png" /></figure>

Time-series of resting-state BOLD

Sequence of time-dependent FC(t)
("dFC stream")

<figure><img src="image_6.png" /></figure>

Recurrence matrix of FC(t)
("dFC matrix")

Rate of variation of FC(t)
(time-series of "dFC speed")

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 2. Empirical vs surrogate dFC streams. We compare dFC streams evaluated from actual empirical fMRI data (A) with surrogate dFC streams evaluated: from phase-randomized BOLD time-series surrogates (B), compatible with a null hypothesis of stationarity of the dFC stream; and from time-shuffled surrogates (C), compatible with an alternative null hypothesis of lack of sequential correlations in the dFC stream. We show in panels D–F representative dFC matrices (for a given window size $τ = 52$ s) for: (D) empirical data; (E) phase-randomized data; and, (F) time-shuffled data. To the right of the dFC matrices we also show a distance-preserving non-linear projection in two dimensions of the associated dFC stream (using the t-SNE algorithm; every dot corresponds to a specific observations of FC(t) and the path connecting the dots indicates the temporal order in which the different networks are sequentially visited). Note that the projections of individual FC(t) networks in panels (D) and (F) are identical by construction but visited in different orders. These projections make visually evident the stochastic walk nature of dFC streams. (G) Shown here are distributions (smoothed kernel-density estimator) of resting-state dFC speeds for empirical and surrogate ensembles (averaged over subjects), pooled over three distinct window-size ranges: long windows (top, 60–210 s); intermediate windows (middle, 15–60 s); and short windows (bottom, 6–15 s). Distributions for surrogate data significantly differ at the global level from distributions for empirical data in most cases (differences between empirical and time-shuffled distributions in red, between empirical and phase-randomized in green color; stars denote significant differences under two-sided Kolmogorov-Smirnov statistics: *, $p < 0.05$; **, $p < 0.01$).</figcaption></figure>

covariance matrix, as required for a null hypothesis of stationarity. Hence, we used the same random sequence of phases for all the time series, i.e. all the regions, per epoch. We adopted the MATLAB code kindly provided by Rikkert Hindriks and incorporated it in a dFC toolbox soon to be released (to be described in a MethodsX article in preparation). This approach guarantees that the covariance matrix of the original empirical data is maintained but destroys systematic deviations from stationarity. The phase-randomized surrogates are thus compatible with a null hypothesis of stationary inter-regional FC. Once phase-randomized BOLD time-series have been generated, the dFC stream can be constructed based on them as for the original data.

We then generated time-shuffled surrogates. After computing a dFC stream on the actual empirical data, we generated a time-shuffled version by randomly permuting the order of the FC(t) timeframes but maintaining them individually unchanged. By this, the mean and variance of each of the FC connections independently are preserved but any sequential correlations are disrupted. Time-shuffled surrogates are thus compatible with a null hypothesis of absence of sequential correlations in the dFC stream.

When presenting results for surrogate ensemble we usually generated 1000 different random realizations of each of the surrogate types for each subject and present average results or over these different realizations, unless differently specified.

## 2.9. Analysis of dFC speeds

Using correlation as a measure of similarity between matrices implied to use the correlation distance between two observations of functional

connectivity as a measure of the amount of change between two FC(t) observations. By measuring the distance between two FC observations separated by a fixed amount of time set to be equal to the window-size W we thus defined the instantaneous global dFC speed as:

$$V_{dFC,t}(t) = 1 - dFC(t, t+\tau)$$

We refer to this dFC speed as “global” because it is evaluated by comparing whole-brain FC matrices. Note that by this definition hte dFC speed depends on the chosen window size $τ$. We improved the dFC speed histogram estimates by increasing the number of sampled speed observations and avoided potential aliasing artifacts due to the use of a single window (Leonardi and Van de Ville, 2015) by pooling window sizes, given that speed distributions for close enough $τ$ were similar. We could realise that for a vast majority of subjects and for most dFC speed bins, as binned counts in the histograms extracted at contiguous window sizes were statistically indistinguishable. For 76 out of the 85 included subjects, the binned dFC speed counts in the histograms extracted at window sizes $τ_k$ and $τ_k + 1$ were statistically indistinguishable (overlapping confidence intervals, Agresti-Coull estimation) for at least 12 out of 20 bins, where $τ_k$ denote the window sizes included (this holds for any of the considered $τ_k$-s, ranging from 3 to 150 TRs for rs scans in the ‘rs-only’ subset or HCP control dataset and task scans in the ‘rs + task’ subset and from 3 to 105 TR in the –shorter– resting-state scans for the ‘rs + task’ subset). Given this substantial degree of redundancy between speed distributions for contiguous window sizes, we chose to pool speed samples to form just three histograms, over three (arbitrary) window ranges: a “short” window-size range, ~6 < $τ$ < ~15 s (3–8 TRs); a “mid” window-size range, ~15 s < $τ$ < ~60 s (9–32 TRs); and a “long” window-size range, ~60 s < $τ$ < ~210 s (33–105
TRs). Window pooling leads to smoother dFC speed histograms. The three different ranges were chosen to construct in each range histograms with close numbers of speed observations after pooling (note that for longer windows less speeds are computed because there are less non-overlapping windows pairs, so that more window sizes must be pooled to reach a number of observations close to the short windows range). A similar window-pooling strategy is used as well in our related study by Lombardo et al. (2020), where we also introduce metrics of dFC speed restricted to specific networks of interest.

Correlations and scatter plots between age or cognitive scores and dFC speeds were constructed based on the median of dFC speed distributions, either computed at single window sizes or pooled, depending on the different analyses. The same procedures were followed for all dFC speed analyses for the actual empirical and for the two types of surrogate data. For comparing dFC speed histograms between empirical and surrogate data (as in Fig. 3 and S3, S4 or S5), we used a same binning in 20 bins for the three types of data for every subject, but adapted the binning to each specific subject’s range, always keeping the same number of bins. Next, we compared subject-by-subject the normalized counts between empirical and surrogate histograms, proceeding from the leftmost to the rightmost bin. Finally, we tested in how many subjects the dFC counts – bin-by-bin – for empirical data were under- or over-estimated (lack of overlap between 95% Coull-Agesti confidence intervals for the count) with respect to histograms for a given surrogate type.

## 2.10. Detrended fluctuation analysis

Detrended Fluctuation Analysis (DFA) allows for detecting intrinsic statistical self-similarity embedded in a 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, the latter requiring however the stationarity of the analysed signal. See Witt and Malamud (2013) or Metzler et al. (2014) for a review. 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.

To perform DFA we first evaluated dFC streams (using a given window size $\tau$) and then evaluated its instantaneous dFC increments:

$$v_{dFC,τ}(t) = 1-dFC(t, t + δt)$$

where $δt$ corresponds here to the minimum possible sliding time-shift in our discretely sampled time-series, i.e. 1 TR (one data point of shift). Note that these instantaneous increments $v_{dFC,τ}(t)$ continue to depend on the window-size $\tau$ because the dFC streams are computed adopting a specific window-size. Therefore, we perform a different fractal scaling analysis over each of the dFC streams evaluated for different window sizes $\tau$.

To perform DFA, we first converted the time-series of instantaneous dFC increments $v_{dFC,τ}(t)$ into an unbounded process:

$$D_{\tau}(t_i) = \sum_{l=1}^{i} v_{dFC, \tau}(t_i)$$

Let $K$ denote the number of samples in the time series, that are split into $M$ non-overlapping segments $q = 1...M$ of length $k$ each, with $M = K/k$. For each segment $q$ the fluctuation strength was computed as the squared difference between $D_{\tau}(t)$ and its trend $D_{\tau}^{(trend)}(t)$ (in the linear case this is the regression line of $D_{\tau}(t)$ over the interval $t = 1...k$):

$$F_{q}^{2}(k) = \frac{1}{k} \sum_{l=0}^{k-1} [D_{\tau}(t_{q+l}) - D_{\tau}^{(trend)}(t_{q+l})]^{2}$$

In the case of scale-free correlation this fluctuation strength scales with segment size $k$. That is, (on average) one finds a linear power law of the form:

$$\log F_{q}(k) = \alpha_{DFA} \log k + C$$

In Fig. S6 we denote $F_{DFA}(k) = \langle F_{q}(k) \rangle$. The scaling parameter $\alpha_{DFA}$ is the primary outcome measure of DFA. In the case of the scale-free processes with the aforementioned power law, $\alpha_{DFA}$ resembles the Hurst exponent (Metzler et al., 2014), leading to the interpretation:

* $0 < \alpha_{DFA} < 0.5$: $v_{dFC,τ}(t)$ displays anti-persistent fluctuations

* $\alpha_{DFA} = 0.5$: $v_{dFC,τ}(t)$ displays uncorrelated Gaussian fluctuations (“white noise” or, equivalently, D$\tau$ resembles Brownian motion)

* $0.5 < \alpha_{DFA} < 1$: $v_{dFC,τ}(t)$ displays persistent fluctuations (approaching “pink noise” when $\alpha_{DFA}$ is close to 1)

* $1 \le \alpha_{DFA}$: $v_{dFC,τ}(t)$ is non-stationary (strictly speaking, DFA is undefined in this case)

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 the output value $\alpha_{DFA}$ could not be interpreted as a scaling exponent. Following (Ton and Daffertshofer, 2016; see https://github.com/marlo w17/FluctuationAnalysis for a dedicated toolbox), 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 fitted 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 $F_{q}(k)$ – beyond its mean value $\langle F_{q}(k) \rangle$ – 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 $k$. To perform model selection, the toolbox then computes the corrected Akaike Information criterion, AIC$_c$, for each one of the tested models:

$$\text{AIC}_c = 2\log \mathcal{L}_{\text{max}} - 2p + \frac{2q(p+1)}{M-p-1}$$

where $p$ is the number of free parameters to fit in the model and the number of points used to estimate the density of $F_{q}(k)$. 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 AIC$_c$ was selected as the relatively best one, and if this was the linear one, the corresponding $\mathcal{L}_{\text{max}}$-fitting parameter was considered as $\alpha_{DFA}$.

We used a range of $10 < k < 80$, 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. A genuine power-law scaling in the DFA of subjects could be found for all subjects in the ‘rs-only’ and ‘rs + task’ subsets in at least 80 out of the 105 different window-sizes (3–105 TRs) used to estimate the dFC streams. Given this general evidence for widespread power-law scaling of the $v_{dFC,τ}(t)$ increments in all subjects, during both resting-state and task scans (apart from sporadic exceptions), as well as in the HCP control data, we computed $\alpha_{DFA}$ exponents in all cases.

## 2.11. dFCwalk toolbox

We have compiled a toolbox of MATLAB® functions to perform most dFC analysis operations described in the present study. It will be described in detail in a MethodsX paper in preparation and, waiting for release, it can be obtained upon reasonable request to the authors.
# 3. Results

## 3.1. Dynamic Functional Connectivity as a stochastic walk

A widespread manner to extract FC from rs- or task-fMRI is to parcellate the brain into *N* macroscopic regions – we follow here a parcellation described first by Desikan et al. (2006) (see Table S1) – and to compute pairwise linear correlations $FC_{ij} = \text{Corr}[a_i(t), a_j(t)]$, between the region-averaged time series of neural activity $a_i(t)$ and $a_j(t)$ of regions *i* and *j*, based on the entire fMRI imaging session (tens of minutes). The result of this procedure is a square *N*-times-*N* matrix of static FC. To estimate how FC fluctuates in time around the session averaged static average, we adopted a common sliding window approach – followed e.g., by Allen et al. (2012) –, repeating the FC construction separately for each time-window (of a fixed duration $\tau$) and generated a time-ordered sequence of $FC(t)$ matrices, or dFC stream (Fig. 1A).

Next, we studied the similarity between time-resolved networks observed at different times. To quantify the amount of FC network variation, we introduced a metric of similarity between FC matrices and evaluated the so-called *dynamic Functional Connectivity (dFC) matrices* (cf. Hansen et al., 2015). The dFC matrix entries $dFC_{ab}$ provide the normalized correlation $corr[FC(t_a), FC(t_b)]$ between any two $FC(t)$ networks observed at times $t_a$ and $t_b$, as depicted in Fig. 1B. A rs-dFC stream for a representative resting-state recording is shown in Fig. 2A, with its associated dFC matrix in Fig. 2B (to the left). From the inspection of this matrix, we can recognize that the rate of variation of $FC(t)$ matrices was not constant along time, but rather heterogeneous. The associated dFC matrix of Fig. 2B (left) displayed characteristic patterns composed out of square-shaped red-hued blocks, corresponding to epochs of transiently increased similarity between consecutive $FC(t)$ network frames. Such epochs of relative FC stability increase –or ‘dFC knots’– were intertwined with transients of relative instability, shown in Fig. 2B (left) by light green or blue stripes in the dFC matrix, denoting strong dissimilarity from previously visited $FC(t)$ networks. During such transients –or ‘dFC leaps’– $FC(t)$ quickly morphed before stabilizing again into the next dFC knot. Knots and leaps could be observed when computing dFC streams and matrices over the whole broad range of window-sizes we tried (between 6 s and 5 min). Additional examples for representative subjects and window sizes are shown in Fig. S2A.

To provide a more quantitative description of the heterogeneous $FC(t)$ change, we determined the rate of change of FC networks along the dFC stream. As said, we considered the dynamics of FC as a stochastic exploration of the space of possible FC configurations and assumed pairs of $FC(t)$ matrices to be separated by an observation window equal to the window-size itself $\tau$ used for $FC(t)$ estimation. We defined the global dFC speed at time $t$ as the quantity $V_{dFC,t}(t) = 1 - \text{corr}[FC(t), FC(t + \tau)]$. As illustrated on the bottom of Fig. 1A, this global dFC speed can be interpreted as the distance travelled in FC space between two ‘stroboscopic’ observations at times $t$ and $t + \tau$ (i.e. over a fixed time interval corresponding to the closest possible interval separating two windows without overlap). Hence, the time-resolved $FC(t)$ matrix may be seen as performing a stochastic walk in the space of possible FC network configurations. The global dFC speed $V_{dFC,t}(t)$ thus informs us about how fast and far away is the time-resolved FC network moving along the stochastic path given by the observed dFC stream, at a time $t$.

We sampled the statistical distributions of $V_{dFC,t}(t)$ for different subjects and values of $\tau$ (see Materials and Methods). We also computed window-pooled dFC speed, simultaneously mixing estimations from different window sizes within three different ranges–long (one to 3 min), intermediate (tens of seconds) and short (6–15 s) window sizes–, also to avoid detection of artefactual fluctuations due to the use of a unique fixed window (Leonardi and Van de Ville, 2015). A distribution of resting-state pooled dFC speeds over all subjects is shown by the blue curve in Fig. 2C, for the three long, intermediate and short window-size ranges (from top to bottom). Analogously, examples of pooled dFC distributions sampled

for single subjects are given in Fig. S2B. For all considered $\tau$-s (both at the group and at the single subject level), these global dFC speed distributions displayed a clear peak near a median value $V_{dFC}$, i.e. the typical dFC speed. The dFC speed distributions were generally skewed and/or kurtotic, deviating from Gaussianity (over three quarters of computed distributions, for different subjects and window sizes, Lilliefors test, $p < 0.05$) and had generally fat tails, in particular evident for long or short window sizes. Before we discuss the implications of these deviations from Gaussianity, we here want to highlight these fat tails indicate that FC reconfiguration events of an anomalously small (low dFC speed) or large (high dFC speed) size are observed with anomalously large probability.

Overall, the FC reconfiguration speed along the dFC stream appeared not to be constant. To show this, we performed a non-linear distance-preserving projection of the sequences of $FC(t)$ matrices observed along the dFC stream into a space of lower dimension. We used two-dimensional projections of dFC streams via t-Stochastic Neighborhood Embedding (t-SNE, Hinton and van der Maaten, 2008). In these plots, each dot corresponds to the projection in two dimensions of a different time-resolved FC network and temporally consecutive dots are linked by a line (see Materials and Methods). The dFC stochastic walk can thus be explicitly visualized, as in the example of Fig. 2B (right), showing the projection of the dFC stream associated to the dFC matrix plotted on the left. In this t-SNE projection, $FC(t)$ within dFC knots form smooth and continuous segments, interrupted by a few cusp points, associated instead to dFC leap events. Here we would like to note that dFC speed fluctuations did not reflect mere head motion artifacts because the size of instantaneous dFC variations (estimated using two alternative forms of motion correction, cf. Materials and Methods) did not correlate with the size of instantaneous head displacements; see Fig. S1A. After motion correction, large dFC speeds could still be detected even in absence of head movement. Conversely, large head movements could occur without big changes of FC. Therefore, dFC fluctuations are most likely not artefactual. This is in line with observations of extensive dFC even in the spontaneous fMRI activity of deeply anesthetized mice in which head was completely fixed (Melozzi et al., 2019).

We finally note, that the fact alone of observing changes in dFC speed along time does not imply that these changes are non-trivial. Indeed, to assess “non-triviality”, one must be able to prove deviations of the statistical properties of empirical dFC streams from trivial null hypotheses. Below we discuss possible deviations from two alternative null hypotheses of “order” and “randomness”. To anticipate, empirical dFC streams have dFC speed distributions lying between these two null hypotheses (Crutchfield, 2011; see Discussion).

## 3.2. Dynamic Functional Connectivity deviates from “order”

$FC(t)$ variability might not be indicative of genuine dynamics but represent fluctuations around an underlying “order” described by an unchanging, static FC. This “order” scenario corresponds to a null hypothesis of FC stationarity. In contrast to this is the possibility that a multiplicity of separate FC states exists. As represented in the graphical cartoon of Fig. 1C, deviations from stationarity would imply that well-defined clusters of $FC(t)$ matrices exist along the stochastic walk described by the dFC stream. The step lengths travelled in FC space –measured by dFC speed– would thus be significantly shorter between $FC(t)$ matrices belonging to a same cluster than the distance travelled between $FC(t)$ matrices belonging to different clusters. A possibility is thus that at least some of the red-hued “knots” observed in dFC matrices (cf. Fig. 2B) correspond to well-separated clusters, associated with distinct FC states. To test for this possibility, we generated phase-randomized surrogate dFC streams, following Hindriks et al. (2016). In this type of surrogates (Fig. 2D), a phase-shift–randomized at every time-point but identical for all regions– is applied to fMRI time-series, such to preserve by construction the static time-averaged FC of the measured session but to destroy any coherent fluctuations around it that may result
in deviations from stationarity (see Materials and Methods and Hindriks et al. (2016)). When computing dFC streams, dFC matrices and dFC speed distributions for phase-randomized surrogates, we found that $FC(t)$ fluctuations were not suppressed. Phase-randomized dFC matrices still appeared to have blocks and stripes reminiscent of “knots” and “leaps” in empirical dFC streams (Fig. 2E, left), generally not occurring, however, at the same times as knots and leaps in the original data. Analogously, when performing a low-dimensional projection via the t-SNE algorithm, phase-randomized dFC streams still gave rise to stochastic exploration paths alternating continuous sections with discrete jumps (Fig. 2E, right). The observed larger-than-average fluctuations of $FC(t)$ were well compatible with the null hypothesis of stationarity. However, when computing the distributions of dFC speed for phase-randomized surrogates, we found that they generally differed from empirical distributions, both at the group level and at the level of single subjects. The green dashed curves in Fig. 2C show distributions of resting-state pooled dFC speeds over phase-randomized surrogates for all subjects, for the three long, intermediate and short window-size ranges (from top to bottom). Once again, these distributions were skewed and/or kurtotic, tending to deviate from Gaussianity (Lilliefors test, $p < 0.01$ for long and short windows, not significant for intermediate windows). Most importantly, these pooled dFC speed distributions for phase-randomized surrogates were statistically different from equivalent empirical distributions for two out of three window-size ranges (two-sided Kolmogorov Smirnov test, Bonferroni corrected, $p < 0.01$ for long windows; $p < 0.05$ for intermediate windows; not significant for short windows). The distribution mode was smaller than for the empirical dFC streams.

We performed comparisons at the level of single-subject dFC speed histograms. As shown in Fig. 3A, the probability of observing speeds in the slow range was significantly smaller in empirical distributions than in phase-randomized distributions for a majority of subjects. At the same time, the probability of observing speeds in the fast range was significantly larger (95% binomial confidence interval speed bin-by-speed bin comparison). We obtained similar results when analyzing single subject distributions of dFC speeds during task blocks, rather than resting state (Fig. S3A). Hence, both at the group and single-subject levels, phase-randomized dFC speed distributions appeared shifted toward slower speeds, with respect to empirical dFC streams. Note that this deviation of empirical dFC speed distributions from the stationary “order” scenario still does not imply that red-hued “knots” in the dFC matrix are actual clusters and that “leaps” are FC state transitions. Indeed these “knots” and “leaps” are visible even in dFC matrices for stationary phase-randomized surrogates (Fig. 2D), as previously mentioned. Furthermore, by construction, phase-randomized surrogates must yield the same covariance matrix of empirical data and, therefore, the level of $FC(t)$ variability for empirical and phase-randomized data should be the same. How can these observations be reconciled with a statistically significant difference in dFC speed distributions? We will comment on this apparent contradiction in the Discussion. However, for now, we would like to remark that deviations from the stationary order expressed by the phase-randomized surrogates do not mean automatically that dFC streams are non-stationary. Still, it means that they are “differently stationary”.

We also note that the results obtained for our original dataset hold even for an independent dataset mediated from the Human Connectome Project (HCP, see Materials and Methods). Indeed, as visible in Fig. S4A, distributions of dFC speeds for phase-randomized surrogates are shifted toward slower speeds than for empirical control data, both at the group level (Fig. S4A) and at the single-subject level (Fig. S4B) and at all window-size ranges.

### 3.3. *Dynamic Functional Connectivity deviates from “randomness”*

The notion of dFC speed does not capture generic variance across $FC(t)$ matrices observed at any point in time but specifically focus on *sequential variability*, i.e. on variability occurring along a time-ordered dFC stream between a matrix $FC(t)$ and a second one at a fixed time-

<figure><img src="image_10.png" /></figure>

**Fig. 3.** Empirical dFC streams lie “between order and randomness”. We compare empirical and surrogate histograms of resting-state dFC speeds, pooled over three different window-size ranges (long, intermediate, and short) at the single-subject level. We perform comparisons speed bin by speed bin (overall 20 speed bins, ranked from lower to higher speeds), checking whether dFC speeds within a given bin are observed with a probability significantly above chance-level (red), significantly below chance level (blue) or compatibly with chance level (green), according to a selected null hypothesis: stationarity for comparison with phase-randomized surrogates; and lack of sequential correlations for comparison with time-shuffled surrogates. (A) Comparison with phase-randomized surrogates indicates that, in empirical resting-state data, lower (higher) than median dFC speeds are often under- (over-) represented. These effects are particularly evident in the long time-windows range (leftmost plot). (B) Comparison with time-shuffled surrogates indicates that, in empirical resting-state data, lower (higher) than median dFC speeds are often over- (under-) represented, i.e. a reverse pattern relative to phase-randomized surrogates. (C–D) When separating subjects into two age groups (younger or older than the median), the comparison patterns revealed by panels A and B are confirmed, but crisper for young subjects and more blurred for older subjects. Overall, if we dub as “order” the null hypothesis of static average FC (i.e. phase-randomization) and as “randomness” the null hypothesis of temporally uncorrelated dFC fluctuations (i.e. time shuffling), the statistics of resting-state dFC fluctuations appear to lie “between order and randomness” (i.e., they are “complex”).

distance $FC(t+t)$. The null hypothesis of stationarity imposes a static covariance matrix but does not make differences about the “when” relatively larger and smaller fluctuations occur. “Knots” and “leaps” are localized in specific time-epochs. Waxing and waning knots have a beginning, a duration and an ending, depicted as blocks in the dFC matrix. To test for sequential aspects in the empirical dFC streams, we designed a second trivial scenario of maximal sequential “randomness”. In this alternative null hypothesis –associated with shuffled surrogate dFC streams–, $FC(t)$ means and variances are preserved, while sequential correlations along the stream are destroyed by randomly shuffling the order of timeframes in the empirical dFC stream (Fig. 2F). Contrary this scenario is the possibility that “flight lengths” in the space of FC matrices are sequentially correlated. For instance, as sketched in the cartoon of Fig. 1D, short steps may be followed by short steps with a larger than chance probability –a property of the stochastic walk known as persistence (Witt and Malamud, 2013; Metzler et al., 2014; see later). If “knot”
epochs are sufficiently long-lasting, FC(t) may evolve considerably –and do so smoothly and very gradually– through the composition of a multiplicity of small steps. In this way, transient slowing downs would not be automatically associated to the emergence of FC clusters (as in Fig. 1C). Distances between FC(t) matrices visited during a slowing-down epoch could indeed be even larger than distances between matrices visited during different slowed-down transients (compare e.g., the distances marked by an exclamation mark in Fig. 1D). However, slowing down transients and transient accelerations would still appear in the dFC matrix as red-hued “knot” blocks and blueish “leap” stripes. Incidentally, this discussion confirms once again that the appearance of “knots” and “leaps” in the dFC matrix is not proof per se of the existence of FC clusters and states. But are the “knots” indicative of significant transient slowing-downs with respect to chance level set by the “randomness” null hypothesis? The full proof will require a direct, explicit study of the persistence of the dFC stochastic walk (see below). Here, we first study differences at the level of the distribution of dFC speeds, as we already did with the “order” null hypothesis.

In Fig. 2E we show the dFC matrix and the two-dimensional projection associated to a typical shuffled dFC stream. The dFC matrix (Fig. 2E, left) appears powdery and scattered without visible knot and leap patterns. This fact is not surprising, because this matrix contains precisely the same entry values of the original empirical dFC matrix in Fig. 2B but in a randomly permuted order. Analogously, the individual FC(t) matrices of the shuffled dFC stream are identical to the matrices composing the empirical dFC stream but appear in a permuted order. As a result, in the t-SNE projection of Fig. 2E (right), the points associated with each of the individual FC(t) timeframes are precisely identical to the ones in Fig. 2B (right). The path linking them is erratic and unstructured, much related to Brownian motion.

When computing the distributions of dFC speed for shuffled surrogates we found that they generally differed from empirical and phase-randomized distributions, both at the group level and at the level of single subjects. The red dashed curves in Fig. 2C show distributions of resting-state pooled dFC speeds over shuffled surrogates for all subjects, for the three long, intermediate and short window-size ranges (from top to bottom). These distributions had a residual skewness but differed significantly from a Gaussian only for short window sizes (Lilliefors test, *p* < 0.05 for short windows, not significant for other windows). Once again, they significantly differed from matching empirical dFC speed distributions in all the window ranges (two-sided Kolmogorov Smirnov, Bonferroni corrected, *p* < 0.01 for long and intermediate windows; *p* < 0.05 for short windows). The distribution modes were larger than for empirical dFC streams. When performing speed bin-by-speed-bin comparisons at the level of single-subject dFC speed histograms, the pattern was reversed with respect to the comparison with phase-randomized surrogates. Notably, as shown in Fig. 3B, the probability of observing speeds in the slow range was significantly larger in empirical distributions than in shuffled distributions for a majority of subjects. At the same time, the probability of observing speeds in the fast range was significantly smaller (95% binomial confidence interval speed bin-by-speed bin comparison). We obtained similar results when analysing single subject distributions of dFC speeds during task blocks, rather than resting state (Fig. S3B).

Both at the group and single-subject levels, shuffled dFC speed distributions appeared shifted toward faster speeds relative to empirical dFC streams. Similar results hold as well for the control HCP resting state fMRI dataset, as shown by Fig. S4A (for group-level differences) and S4C (for single subject-level differences). Summarizing, the results of comparisons of the empirical dFC speed distributions with the two alternative types of surrogates, we conclude that, in both rest and task conditions, the empirical dFC speed distribution lies between the trivial “order” and “randomness” scenarios. It is thus reflecting a non-trivial disorder, a.k.a. Complexity (see Discussion).

### 3.4. Slowing of dFC through the human adult lifespan

Our fMRI dataset included subjects over a wide age range of 18–80 years. This allowed us to demonstrate that dFC stream properties are

affected by “healthy aging”, as none of the subjects were diagnosed with a pathological decline of cognitive abilities. Fig. 4A shows dFC matrices for representative subjects of different ages (see Fig. S2A for different τ-s). It is visually evident that the typical duration of dFC knots varied with subject age, seeming to become longer for older subjects.

To quantify this visual impression, we computed in Fig. 4B resting-state pooled dFC distributions separately for the group of subjects younger (blue curves) and older than median age (magenta curves), for the three long, intermediate and short window ranges (from left to right). For all the three window ranges, the distribution of dFC speeds for older subjects was significantly shifted toward slower values (one-sided Kolmogorov-Smirnov, *p* < 0.01 for long and intermediate windows, *p* < 0.05 for short windows), reflecting the longer duration of the slowing-down epochs associated to knots in the dFC matrix. Similar results held for task dFC speed as well.

At the level of single-subject pooled dFC speed distributions, we tracked for every subject the position of the distribution medians, giving the typical dFC speed $v_{dFC}$ and computed its correlation with subject age over the different window ranges. As shown by the scatter plot in Fig. 4C, $v_{dFC}$ significantly decreased with age (bootstrap with replacement c.i., see caption for *p* values) for all three pooled window ranges, for resting-state –with even stronger correlations– for task dFC speeds. The correlations of dFC speed with age where robust even when considering single window estimations, without pooling. As visible on the top of Fig. 4D (blue ranges), single window $v_{dFC}$-s correlated negatively with subjects age over a broad range of sizes ranging from very-short window sizes (6–15 s) up to window sizes of several minutes.

We also found that these negative correlation between $v_{dFC}$ and age held as well for both phase-randomized and shuffled surrogates, and over broad window ranges well matched the ones found for empirical data (see green and red ranges in Fig. 4D). This suggests that dFC slowing down may be linked to a general reduction of FC(t) variance, which is retained as well when constructing the surrogate dFC streams (see Discussion). More specifically, we found that the differences between surrogate and empirical dFC speed distribution tend to smear out with increasing age. In Fig. 3C and D (for the intermediate window range) and in Fig. S5 (for short and long window ranges), we repeated the speed bin-by-speed bin comparison of empirical with surrogate single subject dFC distributions, separating now, however the subjects into two age groups, younger or older than the median age. The patterns of comparison with surrogate distributions have the same overall directions for younger or older subjects: slower speeds tend to be under-represented (over-represented) and faster speeds over-represented (under-represented) when comparing empirical with phase-randomized (shuffled) dFC distributions. However, for the older subjects, the fraction of subjects for which the bin-by-bin comparisons with surrogates were not substantially increased (cf. the broader green central band in Fig. 3C and D or S5). In a sense, therefore, the dFC speed distributions get more “trivial” with aging (i.e., less complex).

### 3.5. Dynamic Functional Connectivity is an anomalous stochastic walk

As previously said, stochastic processes can be “memory-less” –the next step is uncorrelated from the preceding– or display long-range correlations (Witt and Malamud, 2013; Metzler et al., 2014), positive (persistence) or negative (anti-persistence). The classic “Drunkard’s walk”, associated with Brownian motion (or white noise), corresponds to a Gaussian, memoryless process. However, other stationary stochastic processes can display long-range sequential correlations resulting in different statistics of fluctuation, toward “pink noise” (for positive correlations and persistence) or “blue noise” (for negative correlations and anti-persistence). Such anomalous –i.e., deviating from Gaussianity– are common in a variety of contexts (see Discussion) and has been already identified in resting-state and task fluctuations of both electrophysiological and fMRI signals (Linkenkaer-Hansen et al., 2001; Van de Ville et al., 2010; He, 2014).
<figure><img src="image_11.png" /><figcaption>**Fig. 4.** The speed of dFC streams slows down with aging. (A) We show here representative resting-state dFC matrices for three subjects of different ages. With growing age (top to bottom), red block associated to dFC “knots” become longer-lasting. (B) Correspondingly, the distributions of dFC speeds, over all three studied ranges of window sizes, are shifted toward lower values for older relative to younger subjects (one-sided Kolmogorov-Smirnov statistics: *, *p* < 0.05; **, *p* < 0.01). At the single-subject level, median dFC speeds significantly correlate with subject's age. (C) Scatter plots of single subject age vs median dFC speeds (empirical data), for three pooled window size ranges (long to short, from left to right) and for both resting-state (blue dots) and task (red dots) fMRI scan blocks. Significant correlations with age are found for all three window ranges and for both rest and task dFC speed analyses (bootstrap with replacement confidence intervals for Pearson correlation: *, *p* < 0.05; **, *p* < 0.01; ***, *p* < 0.001). (D) Significant age correlations occur robustly as well for single-window dFC speed estimations over wide continuous ranges of window sizes. We report ranges of window sizes in which statistically significant correlations between median dFC speed and age are detected (a “-” sign indicates negative correlation). To the left, resting-state; to the right, task dFC speeds. Significant correlations (bootstrap, *p* < 0.05) are found not only for empirical data but also for both types of surrogate data in widely overlapping window ranges.</figcaption></figure>

To characterize deviations from Gaussianity in the FC(t) fluctuations along the dFC stream, we used a quantitative approach, Detrended Fluctuation Analysis (DFA, Kantelhardt et al., 2001) to quantify the degree and type of long-range correlations along the dFC stream of FC(t) networks. The DFA procedure quantifies the strength of auto-correlations in a sequence by detecting a power-law scaling – described by a scaling exponent $\alpha_{DFA}$ – in the divergence of a quantity $F_{DFA}(k)$, probing the strength of the fluctuations of the sequence at different scales of observation $k$ (see Material and Methods). A value of $\alpha_{DFA} = 0.5$ corresponds to a Gaussian white noise process, in which the standard deviation of the fluctuations grows as $\sqrt{N}$ after $N$ uncorrelated steps. In contrast, larger values $0.5 < \alpha_{DFA} \le 1$ correspond to anomalously persistent and $0 < \alpha_{DFA} \le 0.5$ anomalously anti-persistent fluctuations. Here, specifically, we measured along the dFC stream sequences of instantaneous increments $v_{dFC,t}(t) = 1 - \text{dFC}(t, t+\delta t)$, where $\delta t$ correspond to the time-step between one FC(t) network frame and the following (i.e., 1 TR). This choice yielded a description of the dFC stream as close to continuum in time as possible. We then applied the DFA procedure on sequences of $v_{dFC,t}$ supplemented by a (Bayesian) model-selection step (see Material and Methods), which allows for discarding subjects, for which a genuine power-law scaling is not present (Ton and Daffertshofer, 2016). Note that the DFA exponents $\alpha_{DFA}$ will continue depending on the window $\tau$ chosen for FC(t) estimation. Indeed, instantaneous increments $v_{dFC,t}(t)$ –analogously to “stroboscopic” dFC speeds $V_{dFC,t}(t)$– are necessarily evaluated on a given input dFC stream. Therefore, different dFC streams are obtained for different choices of $\tau$. It will thus be necessary to study how $\alpha_{DFA}$ depends on $\tau$, since different fluctuation statistics could be found for the different ways of measuring FC(t) (see Discussion for possible alternatives and their drawbacks).

We show in Fig. 5 the results of DFA analyses on our data. First of all, we found robust power-law scaling relations in $v_{dFC,t}$ fluctuations over

most subjects and window sizes. Fig. S6 shows examples of robust power-law scaling under DFA for representative subjects and window sizes $\tau$. We also found that for window sizes $\tau \ge 20s$ (i.e., roughly over the intermediate and large windows ranges) the exponent $\alpha_{DFA}$ was systematically and significantly (lack of overlap between 95% mean c.i.) larger than 0.5, for both empirical resting-state and task data (sample mean of $\alpha_{DFA}$ is given by the blue curves in Fig. 5A), indicating persistence of dFC fluctuations. On the contrary, for window sizes $\tau \le 15s$ (i.e., over the short windows range), the exponent $\alpha_{DFA}$ was significantly smaller than 0.5, denoting anti-persistence. For nearly all window sizes, dFC stream was thus significantly deviating from a Gaussian random walk displaying important long-range correlations. Remarkably, we found the same pattern of deviation from Gaussianity as a function of the chosen window size for the control HCP resting state dataset (see Fig. S4D). In this case, for window sizes $\tau \ge 100s$, the average $\alpha_{DFA}$ exponent indicated a persistence even larger than for our aging dataset.

When performing DFA analysis on surrogate dFC streams, we found by construction that the “randomness” shuffled surrogates corresponded for all probed windows to a Gaussian uncorrelated walk with $\alpha_{DFA} \sim 0.5$ (sample mean given by the red curve in Fig. 5A). Within computational limits, the “order” phase-randomized surrogates (sample mean given by the green curve in Fig. 5A) had an $\alpha_{DFA}$ spectrum statistically indistinguishable from the one of empirical data, in agreement with the literature (Dingwell and Cusumano, 2010), and providing therefore a robust benchmark to probe for eventual correlations between $\alpha_{DFA}$ and behavior or cognition.

Note that the $\alpha_{DFA}$ exponents measured for dFC streams strongly differ for analogous exponents estimated from head motion time-series (cf. Fig. S7A), hinting at their origin in genuine fluctuations of neuronal activity. Also note that deviations from Gaussianity in fluctuations of dFC streams may be linked to the underlying spectral properties
<figure><img src="image_12.png" /><figcaption>**Fig. 5.** The dFC stream is an anomalous stochastic walk. (A) We performed Detrended Fluctuation Analysis (DFA) of the time-series of instantaneous variations of FC(t) along the resting-state and task dFC streams. For every window size, such analysis quantifies whether the stochastic fluctuations of FC(t) are uncorrelated, as in a Gaussian random walk (corresponding to a DFA exponent $α_{DFA} ∼ 0.5$, akin to “white noise”); positively correlated as in a persistent anomalous stochastic walk (corresponding to $0.5 < α_{DFA} < 1$, akin to “pink noise”, observed for most intermediate and long window sizes); or negatively correlated as in an anti-persistent anomalous stochastic walk (corresponding to $α_{DFA} < 0.5$, akin to “blue noise”, observed for short window sizes). Cartoons of these three types of stochastic walk are represented on the top of the plot. Lines report the median $α_{DFA}$ and shaded contours the 95% confidence interval over subjects (median ± 1.96*standard deviation of sample mean). The max $τ$ used for single-window $α_{DFA}$ calculation was shorter than for dFC speed analyses and trimmed to 210 s (excluded $τ$ range shaded in gray). (B) Distributions of resting-state $α_{DFA}$, pooled over all subjects and window sizes within distinct ranges (long windows on top and short windows on bottom) are shifted from persistence or anti-persistence toward the uncorrelated randomness value $α_{DFA} ∼ 0.5$ for older relative to younger subjects (one-sided Kolmogorov-Smirnov statistics: *, $p < 0.05$; **, $p < 0.01$). (C) Significant correlations (bootstrap, $p < 0.05$) between $α_{DFA}$ and age occur over selected ranges of window sizes (a “-” or “+” sign indicate negative or positive correlations). On top, resting state; bottom, task dFC speeds. Robust negative correlations are found only for empirical data. (D) Scatter plots of age vs $α_{DFA}$ (empirical data), for a representative $τ$ in the long window range and for both resting state (blue dots) and task (red dots) fMRI scan blocks.</figcaption></figure>

of fMRI multivariate time-series (cf. Ciuciu et al., 2012, 2014; see Discussion).

### 3.6. Loss of complexity of dFC through the human adult lifespan

We then studied how the scaling properties of dFC stream fluctuations varied with age. We show in Fig. 5B sample resting-state distributions of $α_{DFA}$, pooled over the long (top) and short (bottom) window size ranges, separately for the two groups of subjects with smaller (blue curves) or larger (magenta curves) than median ages. For younger subjects, the distributions for the long (short) window ranges are peaking at $α_{DFA}$ values well above (below) the Gaussian expected value at $~0.5$. For older

subjects, however, the distributions were shifting significantly toward 0.5, indicating Gaussianity (one-sided Kolmogorov Smirnov, Bonferroni corrected, $p < 0.05$ for long windows; $p < 0.01$ for short windows; not significant for intermediate windows). Interestingly, the median $α_{DFA}$ increased with aging in the short windows range, correspondingly to reduced anti-persistence and decreased in the long windows range, reducing persistence. In other words, aging lead invariably to replacing complexity with randomness.

Moving from group to subject-by-subject analyses, Fig. 5C shows the window size ranges for which correlations between subject-specific $α_{DFA}$ and age were significant. For both empirical resting-state and task fMRI, correlations were significantly positive for short window sizes. There were also significantly negative correlations for several windows in the intermediate and long window size ranges. The range of significance was broader for resting-state than for task dFC streams.

We also performed a similar analysis for surrogate dFC streams. There were no correlations with age for shuffled dFC streams for which $α_{DFA}$ was always fluctuating around the Gaussianity value of $~0.5$. For phase-randomized surrogates whose $α_{DFA}$-s were very close to the empirical data, although not precisely identical, we also found significant positive correlations in the short window range. However, the negative correlations with age in the intermediate and long window sizes ranges were tendential but not significant (even before applying multiple comparison correction). The replacement of complexity by randomness is thus more marked in empirical than in phase-randomized data.

Finally, in Fig. S7B we show that the $α_{DFA}$ computed from head motion time-series (cf. Fig. S7A) were not correlating significantly with age. Therefore, variations of $α_{DFA}$ with age for empirical dFC streams do not reflect –or at least, not uniquely– variations of head motion fluctuations (see Discussion).

### 3.7. Dynamic Functional Connectivity correlates with task and cognitive performance

We tested whether dFC stream properties –as its speed $V_{dFC}$ and the scaling properties of its instantaneous increments summarized by $α_{DFA}$– were indicative of behavioural or cognitive performance. We analysed in particular correlations with fluency in a simple visuomotor task or with the score obtained general in a common clinical assessment –the Monreal Cognitive Assessment (MOCA), probing several cognitive domains affected in age-related dementias (Nasreddine et al., 2005). In short, we found that speed of dFC measured during the task –but not at rest– correlated positively with task performance, measured by the spectral overlap (SO, Fig. 6A and Figs. S8A–B). We also found that $α_{DFA}$ measured at rest correlated positively with the general composite MOCA score (Fig. 6B) and that variations of $α_{DFA}$ between rest and task could correlate with visuomotor task PLI (Figs. S8C–D). More specifically, concerning dFC speed analyses, we found that window-range pooled task $V_{dFC}$ correlated with task’s SO for all three long (Fig. S8A), intermediate (Fig. 6A) and short (Fig. S8B) window ranges (bootstrap with replacement c.i., see captions for p values). Correlations were particularly strong in the intermediate window ranges. Considering single-window $V_{dFC,t}$, we also found significant positive correlations over wide ranges of short, intermediate and long window sizes, as shown in Fig. 6C (blue ranges). Therefore, the faster were dFC streams during the task and the best manual movement could remain phase-locked with the displayed visual motion.

The correlations between dFC speed and SO also held for surrogate task dFC streams of both phase-randomized and shuffled types, over window size ranges close to the ones for empirical data (Fig. 6C, red and green ranges). Therefore, as in the case of correlations with age (Fig. 4D), it seems that SO correlations with dFC speed can be accounted for by a general increase of dFC variance, which is preserved when constructing the surrogates dFC streams (see Discussion). Correlations between resting-state dFC speeds and SO were not significant.

Concerning DFA analyses, we found (Fig. 6D) that resting-state $α_{DFA}$ correlated positively with the general MOCA score over a substantial
dFC knots are not necessarily indicating the existence of compact clusters separated by long gaps (cf. the comparison between Fig. 1C–D). Still, it may be that at least some dFC leaps are large enough to be qualified as proper FC state transitions (cf., the over-expression of large speeds in Fig. 3). We insist here once again on the fact that the debate “stationary vs. non-stationary” (see e.g., Hindriks et al., 2016) may not be the more pertinent when attempting to capture the nature of dFC specificities. It instead appears that what makes empirical dFC unique relative to surrogate ensembles are its specific sequential correlation properties, deviating from trivial Gaussian stationarity.

We used fractal scaling analysis to assess “how” empirical dFC streams manifest long-range sequential correlations. We identified for most window sizes deviations of empirical dFC streams from Gaussianity, finding strong evidence for long-range correlations, positive –denoting persistence– or negative –denoting anti-persistence– depending on the window size. We detected a nearly identical window-size dependent spectrum of long-range correlations for a fully independent dataset mediated from the HCP initiative, hinting at the fact that our findings are general for resting state fluctuations of the BOLD signal and not just restricted to our specific study.

A possibility is that higher than Gaussian $\alpha_{DFA}$ exponents are the ultra-
slow range extension or at least the surviving observable shadow of faster
scale-free microstate fluctuation dynamics (Van de Ville et al., 2010) that
cannot be tracked with the too-low temporal resolution of fMRI. We
observed this persistence of dFC streams very robustly over the whole
intermediate and long window size ranges. Over time windows from
~20s up to several minutes, exponents are safely lying above the
Gaussian value of $\alpha_{DFA} \sim 0.5$. The extension of the range for which
persistence holds also means that, although the exact dFC speed and $\alpha_{DFA}$
values may quantitatively change for different window sizes, our random
walk analyses are describing essentially the same qualitative phenomenon,
as long as the adopted window size is larger than $\tau > 20$s.

Long-range persistence has been observed since long in contexts as
diverse as water flooding (Hurst, 1951; Mandelbrot and van Ness, 1968),
fluctuations of musical rhythms (Hennig et al., 2011), foraging in eco-
systems (Viswanathan et al., 1999), prices on the stock market (Lo,
1991), internet traffic (Cleveland and Sun, 2000), etc. In neuroscience;
long-term persistence has been identified in resting-state and task fluc-
tuations of electrophysiological and imaging signals (Linkenkaer-Hansen
et al., 2001; Van de Ville et al., 2010; He, 2014) and found to correlate
with behavior (Palva et al., 2013; Ciuciu et al., 2014). Anti-persistence is
less frequently discussed, but it has also been routinely found, for
instance, once again in nonlinear dynamics (Penna et al., 1995), sportive
scores series (Gabel and Redner, 2011) or in heartbeat (Peng et al.,
1993), more relevant here since it may be a potential contaminating
artefact of physiological but not neural origin in the BOLD signal.

As a matter of fact, our observation of multifractal scaling in dFC is
not completely unexpected. Quite on the contrary, it is well known that
fMRI time-series at rest and during tasks display very characteristic
multifractal spectral properties, relating to both behavioural perfor-
mance and pathological alterations and modulated by task difficulty or
aging (Maxim et al., 2005; Ciuciu et al., 2012; He, 2014; Churchill et al.,
2016; Dong et al., 2018). Tools way more sophisticated than the ones we
adopt here have been used to characterize and confirm multifractality in
fMRI signals (Ciuciu et al., 2017; La Rocca et al., 2018). Furthermore, the
multifractal properties of fMRI signals hold not only at the level of uni-
variate spectra but also of cross-spectra therefore translating into specific
signatures even at the level of static networks, once again with behav-
ioural correlates (Ciuciu et al., 2014). Therefore, our findings reconnect
with solid, cumulating evidence about the fact that fluctuations of co-
ordinated neural activity have non-Gaussian components. Most likely,
the multifractality of the original underlying signals contributes to the
multifractality of the derived dFC streams. It is however only when dis-
cussing multifractality directly at the higher level of time-resolved net-
works, that the interpretation of dFC as an anomalous random walk in FC
space naturally arises. Beyond merely proposing variants of previously

used quantitative biomarkers, the fact itself of lifting up multifractal
analyses from the level of activations to the level of the “chronnectome”
enables a qualitatively new vision of what dFC really is: a non-random,
complex search for a coordinated brain state –captured by the transient
FC(t)– within a large dimensional space of possible dynamic
configurations.

Persistent random walks have been associated to optimality in local
search for resources during foraging in ecosystems (Viswanathan et al.,
1999) but also for “mental space search” of optimal strategies during
bidding (Radicchi et al., 2012). Intriguingly, we observed that higher
$\alpha_{DFA}$ exponents –i.e. deviating more from unstructured Gaussian random
walk and more approaching persistent walks of the Levy type (Metzler
et al., 2014)– were associated with better general cognitive performance
(Fig. 6B–D). Following up on our interpretation of dFC as an anomalous
stochastic walk, we may thus speculate that dFC serves as a neural pro-
cess to efficiently ‘forage’ for cognitive processing resources. Specifically,
the non-trivial fluctuations of $FC(t)$ may play the functional role of
searching for FC patterns adapted to the information routing demands
(Battaglia et al., 2012; Kirst et al., 2016) of ongoing mental
computations.

The anti-persistence of dFC at very short window-sizes also correlated
with cognitive performance. However, this correlation was not anymore
holding at the partial correlation level when we regressed age out as a
common covariate (unlike the stronger positive correlations found for
intermediate and long window ranges, which robustly survived even at
the partial correlation level). This anti-persistence was not specific to
empirical dFC, since it was found as well in phase-randomized surrogates
and varied with age. Since “blue” anti-persistent noise is associated to
heart-beat variability (Peng et al., 1993) and that the fractal scaling of
heart-beat dynamics is also affected by aging (Iyengar et al., 1996) –as
well as in the case of many other physiological signals (Goldberger et al.,
2002)– is likely that dFC estimation over these very short scales is
affected by physiological artifacts of not neural origin. However,
persistence observed for window size $\tau > 20$ s is more robust, and the
measured DFA exponents are very different from the ones concomitantly
measured from head-motion (cf. Fig. S7).

Ultimately, a strong indication of the non-artefactual nature of
empirical dFC is the fact that its random walk properties do correlate to a
certain extent with behavior and cognitive performance. 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 information processing (Salthouse, 1996; Finkel et al., 2007).
Cognitive aging has been associated with deficits in disengaging from
active brain functional states, more than to alterations of the states
themselves (Clapp et al., 2011; Cashdollar et al., 2013). Aging affects dFC
streams by reducing their speed. Thus, it may seem that slowing down of
cognition is paralleled by slowing of dFC. To move beyond mere con-
jectures, we do believe that more and better-adapted experiments should
be designed to probe speed of processing or task switching. Indeed, the
fact that dFC speed correlations with task performance are significant
during task fMRI but only tendential during resting-state may also hint at
the fact that dFC “moves slower” because motion during the task is
different, if not slower. Theoretical studies on task dynamics in behavior
and their neural correlates emphasize the emergence of low-dimensional
subspaces holding task-specific flows (see Huys et al., 2014 for a review),
which act as generative models for cognitive and behavioural processes.
Structured Flows on Manifolds (SFMs) are the mathematical represen-
tation of the deterministic features underlying behavior, including mul-
tistability, convergence/divergence of trajectories, as well as
task-specific stability and robustness (Pillai et al., 2017). They arise
from interactions of coordinated brain activations within the
neuro-skeletonamuscular system including visuomotor tasks (Fink et al.,
2008; Sleimen-Malkoun et al., 2014), multi-limb coordination (Fink
et al., 2000a, 2000b), multisensory integration (Lagarde et al., 2006) and
learning (Schöner et al., 1992; Zanone et al., 1997), but causal contri-
butions are difficult to disentangle. It is tempting to interpret the
“randomness replacing complexity” pattern observed here at the level of DFA analyses as degradations of SFMS and a novel form of de-differentiation, a collective name encompassing several kinds of destructuring and loss of complexity occurring in healthy aging at the level of behavior, cognitive strategies and brain activity (Baltes et al., 1980; Lemaire and Arnaud, 2008; Sleimen-Malkoun, 2014).

Interestingly, we also found that performance in our visuomotor task correlated with the capacity to actively modify the DFA exponent found at rest toward a different functioning value during task (Figs. S8C–D). Our finding suggests that brain systems actively tailor the scaling properties of emergent dFC to the specific task demands and that the capacity to do so is important to explain the achieved performance. Similar results have been found at the level of changes in the scaling properties of static FC correlations (Ciuciu et al., 2014), and we now further extend them at the level of dFC analyses. With our dataset and analyses did not allow for probing these and other hypotheses beyond speculations. In particular, the MOCA global cognitive score provides only a very rough summary characterization of the many facets of cognitive performance and its decline. To probe relation between random walk features of dFC and cognitive performance, one would need a more systematic and controlled cognitive testing setup, probing a spectrum of different cognitive functions (attention, working memory, their executive control, etc.) under different conditions within the same subject. Furthermore, it is unlikely that cognitive performance levels in specific tasks–beyond generic global assessment of performance– is associated to alterations of dFC properties at the whole brain level. Indeed, different specific tasks differentially involve alternative functional sub-networks.

To this heterogeneity of regional involvement, may correspond a parallel heterogeneity of dFC properties, which we fully ignore in this first study. It is known that variance of FC links in time (Chen et al., 2017) or even the fractal scaling properties of fMRI signals (Maxim et al., 2005; Churchill et al., 2016; Dong et al., 2018) are affected heterogeneously across brain regions. In a companion study by Lombardo et al. (2020) we make a further step forward, probing variations of cognitive performance in selective attention and other cognitive functions, before and after sleep deprivation. We do so by developing “modular” (i.e. subnetwork specific) forms of dFC speed analysis. We can thus show that, beyond the first –useful but rough– approximation of whole-brain dFC analyses, distinct subsets of functional links can evolve with different dFC speeds. Furthermore, differential modulations of modular speed can result in differential modulation of performance through different tasks. A key result of Lombardo et al. (2020) is that, when cognitive performance changes are small and more task-specific than the large overall cognitive performance losses linked to MOCA score decreases, then they cannot be captured anymore by whole-brain dFC analyses. They can however still be tracked by modular dFC analyses. The reader is referred to Lombardo et al. (2020) for more details.

Our random walk analyses may be combined with improved ways of estimating dFC streams themselves, for instance, by avoiding the use of sliding windows (Lindquist et al., 2014; Yaesoubi et al., 2018). More importantly, it would be important to investigate potential mechanisms giving rise to complex temporal structure between order and disorder in empirical dFC streams. The emergence of long-range correlations has often been associated with critical dynamics (Linkenkaer-Hansen et al., 2001; Chialvo, 2010), and resting-state is well captured by dynamic mean-field models close to a critical point (Deco et al., 2011). Various proposals have been made for the potential clinical use of dFC (McIntosh and Jirsa, 2019) and related metrics (such as DFA and multiscale entropy) quantifying the changes of flow and manifold structure across pathological conditions. More in general, long-range correlations are also generated by non-linear behavior at the “edge of chaos” (Manneville, 1980; Geisel et al., 1987), well in line with our observation that the resulting dFC lies between order and randomness, as also observed in spontaneous activity at the micro-scale (Clawson et al., 2019). Mean-field whole-brain computational models, besides providing further confirmations of the non-artefactual nature of dFC – virtual brains do not have

blood –, may allow identifying the dynamic and neurophysiological mechanisms behind its random walk properties. Generic whole-brain models are already able to qualitatively reproduce switching dFC (Hansen et al., 2015; Cabral et al., 2017a), but enhanced dynamical complexity will be required to account in silico for the rich non-linearities of empirical dFC revealed by our approach.

Future simulations might be fitted to individual subjects via automated pipelines (Schirner et al., 2015; Proix et al., 2016) to render the dFC trajectory of evolution across aging more quantitatively. Models embedding SC typical of different age classes may reproduce the slowing down and complexity loss of dFC as an emergent by-product of SC ‘disconnection’ itself (Salat, 2011). 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 (Bäckman et al., 2006) or metabolism (Arenaza-Urquijo et al., 2013).

## Author contributions

DB and VJ conceived the novel analytic methods; DB and TB performed imaging data analysis; DL and EH contributed to dFC toolbox development; AD provided expertise and toolbox for DFA and model selection; PR supervised imaging experiments; PR, AD and SC conceived the cognitive and motor experiments; SC performed the cognitive and motor experiments and processed behavioural data; RM and JZ contributed to data pre-processing; DB, TB, DL, EH, SC, AD, RM, JZ, PR and VJ wrote the paper.

## Acknowledgments

This research was supported by the Brain Network Recovery Group (through the James S. McDonnell Foundation) and by the European Union’s Horizon 2020 Framework Program for Research and Innovation under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2). DB acknowledges support from the Mission for Interdisciplinarity of the CNRS, France (Infiniti program 2017–2018, “BrainTime”), from the EU Innovative Training Network “i-CONN” (H2020 ITN 859937). DL has been funded by the Uruguayan National Agency of Research and Innovation (ANII, grant POS EXT 2015 1 123495). P.R. acknowledges the following additional funding sources: H2020 Research and Innovation Action grants VirtualBrainCloud 826421 and ERC 683049; German Research Foundation CRC 1315, CRC 936 and RI 2073/6-1; Berlin Institute of Health & Foundation Charité, Johanna Quandt Excellence Initiative. We thank: Rodrigo Sigala, Sebastian Haufe, Michael Schirner, Simon Rothmeier for help in data acquisition; and Djouya Arbabyazd, Dionysios Perdikis, Rita Sleimen-Malkoun, Annette Witt and Paul Triebkorn for discussions.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.neuroimage.2020.117156.

## References

* Allen, E.A., Damaraju, L., Plis, S.M., Erhardt, E.B., Eichele, T., Calhoun, V.D., 2012. Tracking whole-brain connectivity dynamics in the resting-state. Cerebr. Cortex. https://doi.org/10.1093/cercor/bhs352.

* Andrews-Hanna, J.R., Snyder, A.Z., Vincent, J.L., Lustig, C., Head, D., Raichle, M.E., Buckner, R.L., 2007. Disruption of large-scale brain systems in advanced aging. Neuron 56, 924–935.

* Arenaza-Urquijo, E.M., Landeau, B., La Joie, R., Mevel, K., Mézenge, F., Perrotin, A., Desgranges, B., Bartrès-Faz, D., Eustache, F., Chételat, G., 2013. Relationships between years of education and gray matter volume, metabolism and functional connectivity in healthy elders. Neuroimage 83, 450–457.

* Bäckman, L., Nyberg, L., Lindenberger, U., Li, S.-C., Farde, L., 2006. The correlative triad among aging, dopamine, and cognition: current status and future prospects. Neurosci. Biobehav. Rev. 30, 791–807.
Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Probert Smith, P.J., Woolrich, M., 2014. Fast transient networks in spontaneous human brain activity. *Elife* 3, e01867.

Baltes, B.P., 1980. Integration versus differentiation of fluid/crystallized intelligence in old age. Dev. Psychol. 16, 625–635.

Bassett, D.S., Wymbs, N.F., Porter, M.A., Mucha, P.J., Carlson, J.M., Grafton, S.T., 2011. Dynamic reconfiguration of human brain networks during learning. Proc. Natl. Acad. Sci. U.S.A. 108, 7641–7646.

Battaglia, D., Witt, A., Wolf, F., Geisel, T., 2012. Dynamic effective connectivity of inter-areal brain circuits. *PLoS Comput. Biol.* 8, e1002438.

Betzel, R.F., Byrge, L., He, Y., Goni, J., Zuo, X.-N., Sporns, O., 2014. Changes in structural and functional connectivity among resting-state networks across the human lifespan. *Neuroimage* 102P2, 345–357.

Braun, U., Schäfer, A., Walter, H., Erk, S., Romanczuk-Seiferth, N., Haddad, L., Schweiger, J.I., Grimm, O., Heinz, A., Tost, H., Meyer-Lindenberg, A., Bassett, D.S., 2015. Dynamic reconfiguration of frontal brain networks during executive cognition in humans. Proc. Natl. Acad. Sci. U.S.A. 112, 11678–11683.

Braun, U., Schaefer, A., Betzel, R.F., Tost, H., Meyer-Lindenberg, A., Bassett, D.S., 2018. From maps to multi-dimensional network mechanisms of mental disorders. *Neuron* 97 (1), 14–31.

Brockmann, D., Hufnagel, L., Geisel, T., 2006. The scaling laws of human travel. Nature 439, 462–465.

Brown, R., 1828. A brief account of microscopical observations made in the months of June, July and August, 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies. *Edinburgh New Philosophical J.* 358–371. July-September.

Cabral, J., Kringelbach, M., Deco, G., 2017a. Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: models and mechanisms. *Neuroimage* 160, 84–96.

Cabral, J., Vidaurre, D., Marques, P., Magalhães, R., Moreira, P., Soares, J., Deco, G., Sousa, N., Kringelbach, M., 2017b. Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. *Sci. Rep.* 7, 5135.

Calhoun, V.D., Miller, R., Pearson, G., Adali, T., 2014. The Chronnectome: *Time-Varying Connectivity Networks as the Next Frontier in fMRI Data Discovery*. Neuron **84**, 262–274. https://doi.org/10.1016/j.neuron.2014.10.015.

Cashdollar, N., Fukuda, K., Bocklage, A., Aurtenetxe, S., Vogel, E.K., Gazzaley, A., 2013. Prolonged disengagement from attentional capture in normal aging. *Psychol. Aging* 28, 77–86.

Cavanna, F., Vilas, M.G., Palmucci, M., Tagliazucchi, E., 2018. Dynamic functional connectivity and brain metastability during altered states of consciousness. *Neuroimage* 180, 383–395.

Clawson, W., Vicente, A., Ferraris, M., Bernard, C., Battaglia, D., Quilichini, P., 2019. Computing hubs in the hippocampus and cortex. *Sci. Adv.* 6 (5), eaax4843.

Cleveland, W.S., Sun, D.X., 2000. Internet traffic data. *J. Am. Stat. Assoc.* 95, 979–985.

Chen, J.E., Chang, C., Greicius, M.D., Glover, G.H., 2015. Introducing co-activation pattern metrics to quantify spontaneous brain network dynamics. *Neuroimage* 111, 476–488.

Chen, Y., Wang, W., Zhao, X., Sha, M., Liu, Y., Zhang, X., Ma, J., Ni, H., Ming, D., 2017. Age-related decline in the variation of dynamic functional connectivity: a resting-state analysis. *Front. Aging Neurosci.* 9, 203.

Chettouf, S., Rueda-Delgado, L.M., de Vries, R., Ritter, P., Daffertshofer, A., 2020. Are unimanual movements bilateral? *Neurosci. Biobehav. Rev.* 113, 39–50.

Chialvo, D., 2010. Emergent complex neural dynamics. *Nat. Phys.* 6, 744–750.

Churchill, N., Spring, R., Grady, C., Cimprich, B., Askren, M., Reuter-Lorenz, P., Jung, M., Peltier, S., Strother, S., Berman, M., 2016. The suppression of scale-free fMRI brain dynamics across three different sources of effort: aging, task novelty and task difficulty. *Sci. Rep.* 6, 30895.

Ciuciu, P., Varoquaux, G., Abry, P., Sadaghiani, S., Kleinschmidt, A., 2012. Scale-free and multifractal time dynamics of fMRI signals during rest and task. *Front. Physiol.* 3, 186.

Ciuciu, P., Abry, P., He, B., 2014. Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks. *Neuroimage* 95, 248–263.

Ciuciu, P., Wendt, H., Combrexelle, S., Abry, P., 2017. Spatially regularized multifractal analysis for fMRI data. In: *Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society*, pp. 3769–3772.

Clapp, W.C., Rubens, M.T., Sabharwal, J., Gazzaley, A., 2011. Deficit in switching between functional brain networks underlies the impact of multitasking on working memory in older adults. Proc. Natl. Acad. Sci. U.S.A. 108, 7212–7217.

Cohen, J.R., 2018. The behavioral and cognitive relevance of time-varying, dynamic changes in functional connectivity. *Neuroimage* 180, 515–525.

Crutchfield, J., 2011. Between order and chaos. *Nat. Phys.* 8, 17–24.

Daffertshofer, A., Peper, C.E., Frank, T.D., Beek, P.J., 2000. Spatio-temporal patterns of encephalographic signals during polyrhythmic tapping. *J Human Movement Sci.* 19, 475–498.

Damaraju, E., Allen, E.A., Belger, A., Ford, J.M., McEwen, S., Mathalon, D.H., Mueller, B.A., Pearson, G.D., Potkin, S.G., Preda, A., Turner, J.A., Vaidya, J.G., van Erp, T.G., Calhoun, V.D., 2014. Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. *Neuroimage Clin* 5, 298–308.

Davison, E.N., Turner, B.O., Schlesinger, K.J., Miller, M.B., Grafton, S.T., Bassett, D.S., Carlson, J.M., 2016. Individual differences in dynamic functional brain connectivity across the human lifespan. *PLoS Comput. Biol.* 12 (11), e1005178.

Deco, G., Jirsa, V.K., McIntosh, A.R., 2011. Emerging concepts for the dynamical organization of resting-state activity in the brain. *Nat. Rev. Neurosci.* 12, 43–56.

Desikan, R.S., Ségonne, F., Fischl, B., Quinn, B.T., Dickerson, B.C., Blacker, D., Buckner, R.L., Dale, A.M., Maguire, R.P., Hyman, B.T., Albert, M.S., Killiany, R.J.,

2006. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyrodal regions of interest. *Neuroimage* 31, 968–980.

Dingwell, J., Cusumano, J., 2010. Re-interpreting detrended fluctuation analyses of stride-to-stride variability in human walking. *Gait Posture* 32, 348–353.

Dong, J., Jing, B., Ma, X., Liu, H., Mo, X., Li, H., 2018. Hurst exponent analysis of resting-state fMRI signal complexity across the adult lifespan. *Front. Neurosci.* 12, 34.

Fink, P.W., Foo, P., Jirsa, V.K., Kelso, J.A.S., 2000a. Local and global stabilization of coordination by sensory information. *Exp. Brain Res.* 134, 9–20.

Fink, P.W., Kelso, J.A.S., Jirsa, V.K., de Guzman, G., 2000b. Recruitment of degrees of freedom stabilizes coordination. *J. Exp. Psychol. Hum. Percept. Perform.* 26, 671–692.

Fink, P.W., Kelso, J.A.S., Jirsa, V.K., 2008. Perturbation induced false starts as a test of the Jirsa-Kelso Excitator Model. *J. Mot. Behav.* 41 (2), 147–157.

Finkel, D., Reynolds, C.A., McArdle, J.J., Pedersen, N.L., 2007. Age changes in processing speed as a leading indicator of cognitive aging. *Psychol. Aging* 22, 558–568.

Gabel, A., Redner, S., 2011. Random walk picture of basketball scoring. *J. Quant. Anal. Sports* 8, 6.

Geisel, T., Zacherl, A., Radons, G., 1987. Generic 1/f noise in chaotic Hamiltonian dynamics. *Phys. Rev. Lett.* 59, 2503–2506.

Globk, B., Ponce-Alvarez, A., Gilson, M., Ritter, P., Deco, G., 2017. Resting-state Networks in empirical and simulated dynamic functional. *Neuroimage* 159, 388–402.

Goldberger, A., Amaral, L., Hausdorff, J., Ivanov, P., Peng, C., Stanley, H., 2002. Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. Unit. States Am. 99 (Suppl. 1), 2466–2472.

Golos, M., Jirsa, V., Daucé, E., 2015. Multistability in large scale models of brain activity. *PLoS Comput. Biol.* 11 (12), e1004644.

Gonzalez-Castillo, J., Hoy, C.W., Handwerker, D.A., Robinson, M.E., Buchanan, L.C., Saad, Z.S., Bandettini, P.A., 2015. Tracking ongoing cognition in individuals using brief, whole-brain functional connectivity patterns. Proc. Natl. Acad. Sci. U.S.A. 112, 8762–8767.

Gonzalez-Castillo, J., Bandettini, P.A., 2018. Task-based dynamic functional connectivity: recent findings and open questions. *Neuroimage* 180, 526–533.

Grady, C.L., Garrett, D.D., 2014. Understanding variability in the BOLD signal and why it matters for aging. *Brain Imag. Behav.* 8, 274–283.

Hansen, E.C.A., Battaglia, D., Spiegler, A., Deco, G., Jirsa, V.K., 2015. Functional connectivity dynamics: modeling the switching behavior of the resting-state. *Neuroimage* 105, 525–535.

He, B., 2014. Scale-free brain activity: past, present, and future. Trends Cognit. Sci. 18, 480–487.

He, B., 2018. Robust, transient neural dynamics during conscious perception. Trends Cognit. Sci. 22, 563–565.

Hennig, H., Fleischmann, R., Fredebohm, A., Hagmayer, Y., Nagler, J., Witt, A., Theis, F., Geisel, T., 2011. The nature and perception of fluctuations in human musical rhythms. *PLoS One* 6, e26457.

Hindriks, R., Adhikari, M.H., Murayama, Y., Ganzetti, M., Mantini, D., Logothetis, N.K., Deco, G., 2016. Can sliding-window correlations reveal dynamic functional connectivity in resting-state fMRI? *Neuroimage* 127, 242–256.

Hinton, G., van der Maaten, L., 2008. Visualizing high-dimensional data using t-SNE. *J. Mach. Learn. Res.*: Workshop and Proc. 9, 2579–2605.

Houweling, S., Daffertshofer, A., van Dijk, B.W., Beek, P.J., 2008. Neural changes induced by learning a challenging perceptual-motor task. *Neuroimage* 41, 1395–1407.

Hurst, H.E., 1951. Long-term storage capacity of reservoirs. Trans. Am. Soc. Civ. Eng. 116, 770–799.

Hutchison, R.M., Womelsdorf, T., Allen, E.A., Bandettini, P.A., Calhoun, V.D., Corbetta, M., Penna Della, S., Duyn, J.H., Glover, G.H., Gonzalez-Castillo, J., Handwerker, D.A., Keilholz, S., Kiviniemi, V., Leopold, D.A., de Pasquale, F., Sporns, O., Walter, M., Chang, C., 2013. Dynamic functional connectivity: promise, issues, and interpretations. *Neuroimage* 80, 360–378.

Hutchison, R.M., Morton, J.B., 2015. Tracking the brain's functional coupling dynamics over development. *J. Neurosci.* 35, 6849–6859.

Huys, R., Perdikis, D., Jirsa, V.K., 2014. Functional architectures and structured flows on manifolds: a dynamical framework for motor behavior. *Psychol. Rev.* 121 (3), 302–336.

Iyengar, N., Peng, C.K., Morin, R., Goldberger, A.L., Lipsitz, L.A., 1996. Age-related alterations in the fractal scaling of cardiac interbeat interval dynamics. *Am. J. Physiol. Regul. Integr. Comp. Physiol.* 271, R1078–R1084.

Jones, D.T., Vemuri, P., Murphy, M.C., Gunter, J.L., Senjem, M.L., Machulda, M.M., Przybelski, S.A., Gregg, B.E., Kantarci, K., Knopman, D.S., Boeve, B.F., Petersen, R.C., Jack, C.R., 2012. Non-stationarity in the "resting brain's" modular architecture. *PLoS One* 7, e39731.

Kantelhardt, J., Koscielny-Bunde, E., Rego, H., Havlin, S., Bunde, A., 2001. Detecting long-range correlations with detrended fluctuation analysis. *Physica A* 295, 441–454.

Kantz, H., Schreiber, T., 2004. Nonlinear Time Series Analysis. Cambridge University Press.

Kirst, C., Timme, M., Battaglia, D., 2016. Dynamic information routing in complex networks. *Nat Comms* 7, 11061. https://doi.org/10.1038/ncomms11061.

Kucyi, A., Hove, M.J., Esterman, M., Hutchison, R.M., Valera, E.M., 2017. Dynamic brain network correlates of spontaneous fluctuations in attention. *Cerebr. Cortex* 27, 1831–1840.

Lagarde, J., Kelso, J.A.S., 2006. Binding of movement, sound and touch: multimodal coordination dynamics. *Exp. Brain Res.* 173, 673–688. https://doi.org/10.1007/s00221-006-0401-0.

La Rocca, D., Zilber, N., Abry, P., Wassenhove, V., Ciuciu, P., 2018. Self-similarity and multifractality in human brain activity: a wavelet-based analysis of scale-free brain dynamics. *J. Neurosci.* Methods 309, 175–187.
Laumann, T.O., Snyder, A.Z., Mitra, A., Gordon, E.M., Gratton, C., Adeyemo, B., Gilmore, A.W., Nelson, S.M., Berg, J.J., Greene, D.J., McCarthy, J.E., Tagliazucchi, E., Laufs, H., Schlaggar, B.L., Dosenbach, N.U.F., Petersen, S.E., 2016. On the stability of BOLD fMRI correlations. Cerebr. Cortex. https://doi.org/10.1093/cercor/bhw265.

Lemaire, P., Arnaud, L., 2008. Young and older adults' strategies in complex arithmetic. Am. J. Psychol. 121, 1-16.

Leonardi, N., Van De Ville, D., 2015. On spurious and real fluctuations of dynamic functional connectivity during rest. Neuroimage 104, 430-436.

Liégeois, R., Laumann, T.O., Snyder, A.Z., Zhou, J., Yeo, B.T.T., 2017. Interpreting temporal fluctuations in resting-state functional connectivity MRI. Neuroimage 163, 437-455.

Lim, J., Teng, J., Patanaik, A., Tandi, J., Massar, S.A.A., 2018. Dynamic functional connectivity markers of objective trait mindfulness. Neuroimage 176, 193-202.

Lindquist, M.A., Xu, Y., Nebel, M.B., Caffo, B.S., 2014. Evaluating dynamic bivariate correlations in resting-state fMRI: a comparison study and a new approach. Neuroimage 101, 531-546.

Linkenkaer-Hansen, K., Nikouline, V., Palva, J., Ilmoniemi, R., 2001. Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370-1377.

Liu, X., Duyn, J.H., 2013. Time-varying functional network information extracted from brief instances of spontaneous brain activity. Proc. Natl. Acad. Sci. U.S.A. 110, 4392-4397.

Lo, A., 1991. Long-term memory in stock market prices. Econometrica 59, 1279.

Lombardo, D., Cassé-Perrot, C., Ranjeva, J.P., Le Troter, A., Guye, M., Wirsich, J., Payoux, P., Bartrés-Faz, D., Bordet, R., Richardson, J.C., Felician, O., Jirsa, V., Blin, O., Didic, M., Battaglia, D., 2020. Modular Slowing of Resting-State Dynamic Functional Connectivity as a Marker of Task-specific Cognitive Dysfunction in Sleep Deprivation. NeuroImage. https://doi.org/10.1101/2020.01.17.910810. Preprint available at: bioRxiv.

Mandelbrot, B.B., van Ness, J.W., 1968. Fractional Brownian motions, fractional noises and applications. SIAM Rev. 10, 422-437.

Mandelbrot, B.B., 1983. The Fractal Geometry of Nature. Macmillan.

Manneville, P., 1980. Intermittency, self-similarity and 1/f spectrum in dissipative dynamical systems. J. Phys. 41, 1235-1243.

Matsui, T., Murakami, T., Ohki, K., 2016. Transient neuronal coactivations embedded in globally propagating waves underlie resting-state functional connectivity. Proc. Natl. Acad. Sci. Unit. States Am. 113, 6556-6561.

Maxim, V., Sendur, L., Fadili, J., Suckling, J., Gould, R., Howard, R., Bullmore, E., 2005. Fractional Gaussian noise, functional MRI and Alzheimer's disease. Neuroimage 25, 141-158.

McIntosh, A.R., Jirsa, V.K., 2019. The hidden repertoire of brain dynamics and dysfunction. Network Neuroscience 3 (4), 994-1008. https://doi.org/10.1162/netn-a_00107.

Melozzi, F., Bergmann, E., Harris, J.A., Kahn, I., Jirsa, V., Bernard, C., 2019. Individual structural features constrain the mouse functional connectome. Proc. Natl. Acad. Sci. U.S.A. 116, 26961-26969.

Metzler, R., Jeon, J.-H., Cherstvy, A.G., Barkai, E., 2014. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys. 16, 24128-24164.

Nasreddine, Z.S., et al., 2005. The montreal cognitive assessment, MoCA: a brief screening Tool for mild cognitive impairment. J. Am. Geriatr. Soc. 53, 695-699.

Nobukawa, S., Kikuchi, M., Takahashi, T., 2019. Changes in functional connectivity dynamics with aging: a dynamical phase synchronization approach. Neuroimage 188, 357-368.

Palva, J., Zhigalov, A., Hirvonen, J., Korhonen, O., Linkenkaer-Hansen, K., Palva, S., 2013. Neuronal long-range temporal correlations and avalanche dynamics are correlated with behavioral scaling laws. Proc. Natl. Acad. Sci. U.S.A 110, 3585-3590.

Peng, C., Mietus, J., Hausdorff, J., Havlin, S., Stanley, H., Goldberger, A., 1993. Long-range anticorrelations and non-Gaussian behavior of the heartbeat. Phys. Rev. Lett. 70, 1343-1346.

Penna, T., Oliveira, P., Sartorelli, J., Goncalves, W., Pinto, R., 1995. Long-range anticorrelations and non-Gaussian behavior of a leaky faucet. Phys. Rev. 52, R2168-R2171.

Pillai, A.S., Jirsa, V.K., 2017. Symmetry breaking in space-time hierarchies shapes brain dynamics and behavior. Neuron 94 (5), 1010-1026.

Preti, M.G., Bolton, T.A., Van De Ville, D., 2017. The dynamic functional connectome: state-of-the-art and perspectives. Neuroimage 160, 41-54.

Proix, T., Spiegler, A., Schirner, M., Rothmeier, S., Ritter, P., Jirsa, V.K., 2016. How do parcellation size and short-range connectivity affect dynamics in large-scale brain network models? Neuroimage. https://doi.org/10.1016/j.neuroimage.2016.06.016.

Qin, J., Chen, S.G., Hu, D., Zeng, L.L., Fan, Y.M., Chen, X.P., et al., 2015. Predicting individual brain maturity using dynamic functional connectivity. Front. Hum. Neurosci. 9 https://doi.org/10.3389/fnhum.2015.00418. PMID: 26236224.

Radicchi, F., Baronchelli, A., Amaral, L.A.N., 2012. Rationality, irrationality and escalating behavior in lowest unique bid auctions. PloS One 7, e29910.

Salat, D.H., 2011. The declining infrastructure of the aging brain. Brain Connect. 1, 279-293.

Salthouse, T.A., 1996. The processing-speed theory of adult age differences in cognition. Psychol. Rev. 103, 403-428.

Schirner, M., Rothmeier, S., Jirsa, V.K., McIntosh, A.R., Ritter, P., 2015. An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. Neuroimage 117, 343-357.

Schlesinger, K.J., Turner, B.O., Lopez, B.A., Miller, M.B., Carlson, J.M., 20. Age-dependent changes in task-based modular organization of the human brain. Neuroimage 46, 741-762. https://doi.org/10.1016/j.neuroimage.2016.09.001.

Schöner, G., Zanone, P., Kelso, J.A.S., 1992. Learning as change of coordination dynamics: theory and experiment. J. Mot. Behav. 24, 29-48. https://doi.org/10.1080/0022895.1992.9941599.

Schurger, A., Sarigiannidis, I., Naccache, L., Sitt, J., Dehaene, S., 2015. Cortical activity is more stable when sensory stimuli are consciously perceived. Proc. Natl. Acad. Sci. U.S.A 112 (16), E2083-E2092.

Shakil, S., Lee, C.-H., Keilholz, S.D., 2016. Evaluation of sliding window correlation performance for characterizing dynamic functional connectivity and brain states. Neuroimage 133, 111-128.

Shine, J., Koyejo, O., Bell, P., Gorgolewski, K., Gilat, M., Poldrack, R., 2015. Estimation of dynamic functional connectivity using Multiplication of Temporal Derivatives. Neuroimage 122, 399-407.

Shine, J.M., Bissett, P.G., Bell, P.T., Koyejo, O., Balsters, J.H., Gorgolewski, K.J., Moodie, C.A., Poldrack, R.A., 2016. The dynamics of functional brain networks: integrated network states during cognitive task performance. Neuron 92, 544-554.

Sleimen-Malkoun, R., Temprado, J., Hong, S., 2014. Aging induced loss of complexity and dedifferentiation: consequences for coordination dynamics within and between brain, muscular and behavioral levels. Front. Aging Neurosci. 6, 140.

Tagliazucchi, E., Wegner von, F., Morzelewski, A., Brodbeck, V., Laufs, H., 2012. Dynamic BOLD functional connectivity in humans and its electrophysiological correlates. Front. Hum. Neurosci. 6, 339.

Termenon, M., Jaillard, A., Delon-Martin, C., Achard, S., 2016. Reliability of graph analysis of resting state fMRI using test-retest dataset from the Human Connectome Project. Neuroimage 142, 172-187.

Thompson, W., Fransson, P., 2016. Bursty properties revealed in large-scale brain networks with a point-based method for dynamic functional connectivity. Sci. Rep. 6 (1), 39156.

Ton, R., Daffertshofer, A., 2016. Model selection for identifying power-law scaling. Neuroimage 136, 215-226.

Van de Ville, D., Britz, J., Michel, C.M., 2010. EEG microstate sequences in healthy humans at rest reveal scale-free dynamics. Proc. Natl. Acad. Sci. U.S.A 107, 18179-18184.

Vidaurre, D., Quinn, A.J., Baker, A.P., Dupret, D., Tejero-Cantero, A., Woolrich, M.W., 2016. Spectrally resolved fast transient brain states in electrophysiological data. Neuroimage 126, 81-95.

Viswanathan, G.M., Buldyrev, S.V., Havlin, S., da Luz, M.G., Raposo, E.P., Stanley, H.E., 1999. Optimizing the success of random searches. Nature 401, 911-914.

Viviano, R., Raz, N., Yuan, P., Damoiseaux, J., 2017. Associations between dynamic functional connectivity and age, metabolic risk, and cognitive performance. Neurobiol. Aging 59, 135-143.

Witt, A., Malamud, B., 2013. Quantification of long-range persistence in geophysical time series: conventional and benchmark-based improvement techniques. Surv. Geophys. 34 (5), 541-651.

Yaesoubi, M., Miller, R.L., Calhoun, V.D., 2015. Mutually temporally independent connectivity patterns: a new framework to study the dynamics of brain connectivity at rest with application to explain group difference based on gender. Neuroimage 107, 85-94.

Yaesoubi, M., Adali, T., Calhoun, V., 2018. A window-less approach for capturing time-varying connectivity in fMRI data reveals the presence of states with variable rates of change Human. Brain Mapping 39, 1626-1636.

Yang, A.C., Huang, C.C., Yeh, H.L., Liu, M.E., Hong, G.J., Tu, P.C., et al., 2013. Complexity of spontaneous BOLD activity in default mode network is correlated with cognitive function in normal male elderly: a multiscale entropy analysis. Neurobiol. Aging 34, 428-438.

Zalesky, A., Fornito, A., Cocchi, L., Gollo, L., Breakspear, M., 2014. Time-resolved resting-state brain networks. Proc. Natl. Acad. Sci. U.S.A 111, 10341-10346.

Zalesky, A., Breakspear, M., 2015. Towards a statistical test for functional connectivity dynamics. Neuroimage 114, 466-470.

Zanone, P.G., Kelso, J.A.S., 1997. Coordination dynamics of learning and transfer: collective and component levels. J. Exp. Psychol. Hum. Percept. Perform. 23 (5), 1454-1480.

Zimmermann, J., Ritter, P., Shen, K., Rothmeier, S., Schirner, M., McIntosh, A.R., 2016. Structural architecture supports functional organization in the human aging brain at a regionwise and network level. Human Brain Mapping. https://doi.org/10.1002/hbm.23200.