Dr Philippe Ciuciu
email: philippe.ciuciu@cea.fr

Patrice Abry

Biyu J He

Interplay between Functional Connectivity and Scale-Free Dynamics in Intrinsic fMRI Networks

Keywords: Scale-free dynamics, Cross-temporal dynamics, fMRI, Intrinsic brain activity, Task modulation

Studies employing functional connectivity-type analyses have established that spontaneous fluctuations in functional magnetic resonance imaging (fMRI) signals are organized within large-scale brain networks. Meanwhile, fMRI signals have been shown to exhibit 1/f-type power spectra -a hallmark of scale-free dynamics. We studied the interplay between functional connectivity and scale-free dynamics in fMRI signals, utilizing the fractal connectivity framework -a multivariate extension of the univariate fractional Gaussian noise model, which relies on a wavelet formulation for robust parameter estimation. We applied this framework to fMRI data acquired from healthy young adults at rest and performing a visual detection task. First, we found that scaleinvariance existed beyond univariate dynamics, being present also in bivariate crosstemporal dynamics. Second, we observed that frequencies within the scale-free range do not contribute evenly to inter-regional connectivity, with a systematically stronger contribution of the lowest frequencies, both at rest and during task. Third, in addition to a decrease of the Hurst exponent and inter-regional correlations, task performance modified cross-temporal dynamics, inducing a larger contribution of the highest frequencies within the scale-free range to global correlation. Lastly, we found that across individuals, a weaker task modulation of the frequency contribution to inter-regional Rev.#1, Q5 Rev.#1, Q5 connectivity, and thus a weaker recourse to high frequencies,was associated with bet-

ter task performance manifesting as shorter and less variable reaction times. These findings bring together two related fields that have hitherto been studied separatelyresting-state networks and scale-free dynamics, and show that scale-free dynamics of human brain activity manifest in cross-regional interactions as well.

Introduction

In recent years, functional-connectivity analysis applied to resting-state blood oxygen level-dependent (BOLD) fMRI has revealed a rich intrinsic functional architecture of brain activity, manifesting as large-scale, coherent brain networks that recapitulate the spatial patterns of task activations [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Xiong | Inter-regional connectivity to primary motor cortex revealed using mri resting state images[END_REF][START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Cordes | Frequencies contributing to functional connectivity in the cerebral cortex in "resting-state" data[END_REF][START_REF] Raichle | A default mode of brain function[END_REF][START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF][START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF]. The functional significance of fMRI resting-state networks (RSNs) has been demonstrated in various neurological and psychiatric diseases by showing that the degree of disruption of resting-state networks (RSNs) correlated with the severity of the disorder [START_REF] He | The role of impaired neuronal communication in neurological disorders[END_REF][START_REF] Zhang | Disease and the brain's dark energy[END_REF][START_REF] Zhang | Disrupted functional connectivity of the anterior cingulate cortex in cirrhotic patients without overt hepatic encephalopathy: a resting state fMRI study[END_REF]. Rev.#1, Q6 Rev.#1, Q6 Moreover, repetitive training over period of days sculpts spontaneous activity of the resting human brain, suggesting dynamic reconfiguration of RSNs [START_REF] Lewis | Learning sculpts the spontaneous activity of the resting human brain[END_REF].

Most studies assessing functional connectivity so far have used either a seed-based region-of-interest (ROI) approach, in which the time series associated with a chosen ROI is used as a regressor to identify regions of similar temporal behavior across the brain [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Raichle | A default mode of brain function[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: a network analysis of the default mode hypothesis[END_REF], or an independent component analysis (ICA) -an exploratory approach for identifying spatial regions with temporally coordinated activity [START_REF] Beckmann | Investigations into resting-state connectivity using independent component analysis[END_REF][START_REF] Calhoun | A method for making group inferences from functional MRI data using independent component analysis[END_REF][START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group PICA[END_REF]. Each approach relies on either anatomically or statistically driven a priori assumptions (see [START_REF] Cole | Advances and pitfalls in the analysis and interpretation of resting-state FMRI data[END_REF]) for a general review of the pros and cons of both approaches). The seed-based method relies largely on the computation of linear Pearson correlation coefficients between the temporal fluctuations of BOLD signal in different brain regions. The ICA approach identifies spatial components that are maximally independent, each component grouping voxels with similar temporal dynamics [START_REF] Beckmann | Investigations into resting-state connectivity using independent component analysis[END_REF][START_REF] Cole | Advances and pitfalls in the analysis and interpretation of resting-state FMRI data[END_REF].

In nature, Pearson's correlation is linear, static and global. Whether linear coupling is sufficient in describing interactions between brain regions or networks has been studied with care, using the functional integration index and mutual information [START_REF] Hlinka | Functional connectivity in resting-state fMRI: is linear correlation sufficient?[END_REF][START_REF] Marrelec | Regions, systems, and the brain: hierarchical measures of functional integration in fMRI[END_REF]. Definitive answers are still lacking, which depends on whether fMRI data are well modeled as Gaussian processes, and whether their dependence structure can be described by the sole correlation coefficient. It has been observed that depending on the spatial scale at which the correlation measure is assessed, departure from Gaussianity may be relatively minor (within-network) [START_REF] Hlinka | Functional connectivity in resting-state fMRI: is linear correlation sufficient?[END_REF] or significant (between-network) [START_REF] Marrelec | Regions, systems, and the brain: hierarchical measures of functional integration in fMRI[END_REF].

Pearson's correlation can also be considered a static measure of dependency since it does not provide practitioners with any information regarding the contributions of the different frequencies to correlation. To bridge that gap, the coherence function can be used to measure the relative contributions of the different frequencies to correlation [START_REF] Salvador | Undirected graphs of frequencydependent functional connectivity in whole brain networks[END_REF][START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF].

The global nature of Pearson's correlation prevents assessment of dependencies that vary over time. To overcome that limitation, local correlation coefficients can be computed via sliding windows [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF][START_REF] Hutchison | Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques[END_REF] to access dynamic functional connectivity in humans or animals [START_REF] Hutchison | Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques[END_REF][START_REF] Majeed | Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans[END_REF]. It is however natural and efficient to combine local (time-varying) and frequency-dependent correlation measure into time-frequency or wavelet-based measures of correlation such as the wavelet transform coherence [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]. Indeed, temporal reconfigurations of fMRI RSNs have been recently observed over typical scan durations (several minutes) using time-resolved acquisitions and a cascade of spatial and temporal ICA [START_REF] Smith | Temporally-independent functional modes of spontaneous brain activity[END_REF], or sliding-window ICA or principal component analysis (PCA) [START_REF] Kiviniemi | A sliding time-window ica reveals spatial variability of the default mode network in time[END_REF][START_REF] Leonardi | Principal components of functional connectivity: A new approach to study dynamic brain connectivity during rest[END_REF]. More recently, it has been demonstrated that the spatial signature of RSNs can be reconstructed from a few spontaneous coactivation patterns

Methods

Functional magnetic resonance imaging (fMRI) data were acquired from seventeen normal healthy young adults (9 females, age 18-27 years) at rest and performing a visual detection task. All subjects provided written informed consent. This data set has been previously published in separate studies [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] He | Spontaneous and task-evoked brain activity negatively interact[END_REF][START_REF] Fox | Intrinsic fluctuations within cortical systems account for intertrial variability in human evoked brain responses[END_REF]. Detailed methods can be found in SI Methods.

fMRI data acquisition

Blood-oxygen-level dependent (BOLD) fMRI data (4×4×4 mm 3 voxels, TE 25 ms, TR 2.16 s) were acquired in 17 normal right-handed young adults (9 females, age 18-27 years) using a 3T Siemens Allegra MR scanner. All subjects gave informed consent in accordance with guidelines set by the Human Studies Committee of Washington University in St. Louis. Each subject completed 8 fMRI runs, each 194 frames (∼ 7 min) in duration. They consisted of two alternating run types. The first run type was a resting-state study in which a white crosshair was presented in the center of a black screen. Subjects were instructed to look at the crosshair, remain still, and to not fall asleep. The second run type was a task study in which the identical crosshair was presented, but now it occasionally changed from white to dark gray for a period of 250 ms, at times unpredictable to the subjects. The subjects were instructed to press a button with their right index finger as quickly as possible when they saw the crosshair dim. Each of these button-press runs contained 20 crosshair dims time-locked to the scanner TR, with an inter-trial interval of 17.3-30.2 s. Subjects practiced this buttonpress task once in the scanner, prior to the onset of the functional scans. Anatomical MRI included a high-resolution (1×1×1.25 mm 3 ) sagittal, T 1 -weighted MP-RAGE (TR 2.1 s, TE 3.93 ms, flip angle = 7°) and a T 2 -weighted fast spin-echo scan. This data set has been previously published in separate studies [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] He | Spontaneous and task-evoked brain activity negatively interact[END_REF][START_REF] Fox | Intrinsic fluctuations within cortical systems account for intertrial variability in human evoked brain responses[END_REF]. All analyses were carried out using custom-written codes in C++ and Matlab.

fMRI data preprocessing

fMRI preprocessing steps included, i.) compensation of systematic, slice-dependent time shifts, ii.) removal of systematic odd-even slice intensity difference due to interleaved acquisition (slice-timing correction); iii.) rigid body correction for inter-frame head motion within and across runs, and iv.) intensity scaling to yield a whole-brain mode value of 1000 (with a single scaling factor for all voxels). Atlas registration was achieved by computing affine transforms connecting the fMRI run first frame (averaged over all runs after cross-run realignment) with the T 2 -and T 1 -weighted structural images [START_REF] Ojemann | Anatomic localization and quantitative analysis of gradient refocused echo-planar fMRI susceptibility artifacts[END_REF]. Our atlas representative template included MP-RAGE data from 12 normal individuals and was made to conform to the 1988 Talairach atlas [START_REF] Talairach | Co-Planar Stereotaxic Atlas of the Human Brain. 3-Dimensional Proportional System : An Approach to Cerebral Imaging[END_REF]. Data were resampled to 3×3×3 mm 3 voxels after atlas registration. The first four frames of each fMRI run were discarded in all further analyses. The fMRI time courses from each run were made zero-mean and the linear trend was removed. Lastly, the effect of head motion and its temporal derivative were removed by linear regression.

It is known that sudden head movements (like "spikes") may have a strong influence in the estimated scaling exponents. To cope with this issue, a recent approach has been proposed in [START_REF] Siegel | Statistical improvements in functional magnetic resonance imaging analyses produced by censoring high-motion data points[END_REF]. It consists of erasing segments of time series which are corrupted by very large head movements. This methodology turned out to be robust for Hurst exponent analysis at least using detrended fluctuation analysis (DFA) [START_REF] Chen | Effect of nonstationarities on detrended fluctuation analysis[END_REF]. In this study, given the short length of time series, we first investigated the presence of large head movements. As detailed in Appendix A, we found that there were very few movement spikes, such that the removal of temporal segments containing them was unnecessary. Moreover, as detailed later, we rely on a wavelet framework for scaling exponent estimation, which further brings extra robustness against non-stationarities, compared to DFA cf. e.g., [START_REF] Torres | Comparison of different methods for computing scaling parameter in the presence of trends[END_REF].

Rev.#2, Q6 Rev.#2, Q6

Definition of regions of interest (ROIs)

31 ROIs were obtained from a previous task-related functional neuroimaging studies or generated using coordinates from published fMRI studies, which included 10 pairs of homologous brain regions. These ROIs were the same as used in our previous study [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF]. Their locations in the brain are shown in Fig. 1A (mapping to brain surface was done in CARET 1 ). The regions were grouped into five cortical networks based on their known anatomical/functional properties (including attention, default-mode, motor, saliency and visual networks) and a separate group outside the neocortex including hippocampus, thalamus and cerebellum.

The anatomical locations, Talairach coordinates, references and associated networks of these ROIs are listed in Table 1. Specifically, the attention, motor, visual, thalamus and cerebellum regions were obtained from functional studies conducted in [START_REF] He | Breakdown of functional connectivity in frontoparietal networks underlies behavioral deficits in spatial neglect[END_REF]. The default network regions were obtained from task-deactivation patterns from a meta-analysis of nine PET studies, which originally unveiled the default network [START_REF] Shulman | Common blood flow changes across visual tasks: II. Decreases in cerebral cortex[END_REF]. To generate these ROIs, following methods described in [START_REF] He | Breakdown of functional connectivity in frontoparietal networks underlies behavioral deficits in spatial neglect[END_REF], the activation or deactivation Z-score maps were subjected to an automatic peak search, peaks closer than 10 mm were consolidated by averaging their coordinates, and ROIs were defined around peaks by thresholding the map to yield regions of approximately 905 mm 3 , similar size as the coordinates-derived ROIs described below.

The dosolateral prefrontal cortex (DLPFC), part of the frontoparietal attention network, and the saliency (also called "core task-control") network regions were obtained from published coordinates in three studies [START_REF] Dosenbach | A core system for the implementation of task sets[END_REF][START_REF] Seeley | Dissociable intrinsic connectivity networks for salience processing and executive control[END_REF][START_REF] Vincent | Evidence for a frontoparietal control system revealed by intrinsic functional connectivity[END_REF]. The coordinates for Broca's area and the hippocampal formation (HF) were obtained from [START_REF] Embick | Mapping syntax using imaging: problems and prospects for the study of neurolinguistic computation[END_REF] and [START_REF] Vincent | Coherent spontaneous activity identifies a hippocampal-parietal memory network[END_REF], respectively. In cases where coordinates from multiple studies were obtained for one ROI, such as the R DLPFC and R TPJ (Table 1), the center-of-mass of these coordinates were used. A 6-mm-radius sphere ROI centered at these coordinates was created for each region. All regions used in the present study have been investigated in seed-based functional connectivity analyses applied to resting-state fMRI data by the author and, for ROIs in the attention and saliency networks as well as the HF, also in previous published studies [START_REF] He | Breakdown of functional connectivity in frontoparietal networks underlies behavioral deficits in spatial neglect[END_REF][START_REF] Seeley | Dissociable intrinsic connectivity networks for salience processing and executive control[END_REF][START_REF] Vincent | Evidence for a frontoparietal control system revealed by intrinsic functional connectivity[END_REF][START_REF] Vincent | Coherent spontaneous activity identifies a hippocampal-parietal memory network[END_REF][START_REF] Dosenbach | Distinct brain networks for adaptive and stable task control in humans[END_REF] and have yielded networks consistent with those reported in the literature [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF][START_REF] Fox | Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems[END_REF].

Table 1: Anatomical information and references for each ROI. Note that the ROI appearance order below defines the order of ROI entries in the following functional connectivity matrices such as Fig. 1B.

Network

ROI Anatomical location Talairach coordinates

Attention (AN) vIPS [59] (L and R) ventral intraparietal sulcus -24, -69, 30 30, -80, 16

pIPS [59] (L and R) posterior intraparietal sulcus -25, -63, 47 23 -65 48 R TPJ R temporoparietal junction 49, -50, 28

MT [59] (L and R) middle temporal region -43, -70 -3 42 -68 -6

FEF [59] (L and R) frontal eye field -26, -9, 48 32, -9, 48 R DLPFC [61, 62, 63] R dorsolateral prefrontal cortex 43, 22, 34 Default Mode (DMN) AG [60] (L and R) angular gyrus -51, -54, 30 45, -66, 27 SFG [60] (L and R) superior frontal cortex -15, 33, 48 18, 27, 48 PCC [60] posterior cingulate cortex -6, -45, 33 MPF [60] medial prefrontal cortex -6, 51, -9 FP [60] frontopolar cortex -3, 45, 36

Motor (MN)

L SII [59] L second somatosensory area -57, -27, 21 L motor [59] L primary motor cortex -39, -27, 48 Broca [64] Broca's area -42, 13, 14

Non-Neocortical (NC) Thalamus citeHe07b (L and R) thalamus -15, -21, 6 9, -18, 9 R Cerebellum [59] R Cerebellum 21, -54, -21 HF [65] (L and R) hippocampal formation -21, -25, -14 23, -23, -14

Saliency (SN)
R FI [62] R frontoinsular cortex 36, 21, -6 dACC [61] dorsal anterior cingulate cortex -1, 10, 46 Visual (VN) vRetino [59] (L and R) ventral retinal region -15, -75, -9 15, -75, -9 dRetino [59] (L and R) dorsal retinal region -6, -75, 9 9, -75, 12

Scale-free temporal dynamics modeling

Scale-free temporal dynamics is now a commonly observed property in brain activity [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF]. To account for scale-free temporal dynamics in brain dynamics, these references propose quasi-exclusively, either implicitly or explicitly, to use fractional Gaussian noise (fGn), the celebrated model put forward by Mandelbrot [START_REF] Mendelbrot | Fractional brownian motions, fractional noises and applications[END_REF] and massively used in many other scientific fields (see [START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF] for a review). In essence, fGn assumes that data have Gaussian marginal distributions and a power-law type spectral behavior, across a large range of frequencies:

Γ X (f ) ≈ C|f | -α , f m |f | f M , f M /f m 1, with 0 < α < 1. ( 1
)
where α is related to the Hurst parameter H as α = 2H -1. This model is relevant when analyzing brain activity measured from univariate time series, each associated with a given region of interest. However, to assess functional connectivity, it is needed that a collection of time series each associated with a different region of interest, are studied jointly (or simultaneously), to measure for instance how they correlate one to another. It is thus natural to make use of a model inspired from the multivariate extension of fGn (mfGn), proposed e.g., in [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF][START_REF] Amblard | Identification of the multivariate fractional brownian motion[END_REF]. In essence, this model assumes joint Gaussianity for the time series and power-law behaviors both for the auto-and cross-spectra, across a large range of frequencies.

For the sake of simplicity, the model is stated here in the bivariate case only, with a straightforward multivariate extension. Let X and Y denote two time series associated with two brain regions. Their auto-and cross-spectra are defined as, for

f m |f | f M , f M /f m 1: Γ X (f ) = ω X |f | -α X , Γ Y (f ) = ω Y |f | -α Y , Γ XY (f ) = ω XY |f | -α XY (2)
with parameters α X = 2H X -1 and α Y = 2H Y -1 confined to the range [0, 1]. Two important notes are now in order. First, mfGn theoretically further imposes that

α XY = (α X + α Y )/2 = H X + H Y -1.
In the present work, we allow α XY to depart from (H X + H Y ) -1 this is why the model used here is not strictly mfGn but rather inspired from. Second, the theoretical definitions of both fGn and mfGn, actually imply that their spectra exhibit power-law behavior in the limit of low frequencies: |f | → 0. Practically, however, power law behaviors are often assumed to hold across a large but finite range of frequencies:

f m |f | f M (with possibly f m > 0)
. We will stick to that standard practice, while, with a slight abuse of language, continuing to refer to these processes as fGn and mfGn. While in the univariate case, this has little impact on the actual use of fGn as a model for real data, this is of much larger importance in the multivariate setting as this allows theoretically both positive and negative departures of

γ XY = α XY -(H X + H Y -1)
from 0, thus providing us with a significant gain in versatility for analyzing spontaneous brain activity.

Coherence function

In the classical assessment of linear dependencies, the correlation coefficient ρ XY is used. It consists of a static quantity that conveys no information related to the way the different frequencies contribute to the global correlation of X and Y . To overcome that limitation, the coherence function can naturally be used. It consists of a (sort of) frequency-dependent correlation coefficient [START_REF] Marple | Digital spectral analysis with applications[END_REF]:

C XY (f ) = |Γ XY (f )| Γ X (f )Γ Y (f ) . (3)
By definition, it takes values between -1

C XY (f ) 1 and ω XY / √ ω X ω Y is pro- portional to ρ XY . When |C XY (f )| = |ρ XY |,

Fractal connectivity

When X and Y follow the bivariate correlated scale-free model defined in Eq. ( 2), the coherence function becomes, for

f m |f | f M : C XY (f ) = ω XY √ ω X ω Y |f | -γ XY , where γ XY = α XY -(H X + H Y ) + 1, ω XY √ ω X ω Y ∝ ρ XY . (4)
Fractal connectivity is theoretically defined as the case where C XY (f ) reduces to a (non-zero) constant over the range

f m |f | f M , i.e.,
ρ XY = 0 and γ XY ≡ 0.

(

) 5
The intuition underlying fractal connectivity is that, for scale-free data, all frequencies are contributing to the correlation (and hence to functional connectivity) in an equivalent manner, or in a mfGn-type compatible manner. In that case, the coherence function does not bring any extra information compared to the sole static correlation coefficient. Conversely, assuming ρ XY = 0, γ XY > 0 (resp., γ XY < 0) indicates that low frequencies contribute more (resp., less) to correlation of X and Y than do high frequencies. Therefore estimating γ XY and hence the coherence function brings complementary information related to the way temporal dynamics contribute to functional connectivity, compared to the sole static correlation coefficient ρ XY .

Note that, in that context, low and high frequencies are defined in a relative manner: First, the range of frequencies, f m |f | f M , where scale-free properties are observed is estimated. Second, low and high refer respectively to the lower and upper sub-ranges of that scale-free range of frequencies. Two interesting limit behaviors are worth being described:

1. When |ρ XY | → 1 then necessarily γ XY → 0 because 0 |C XY (f )| 1, ∀f ;
2. When |ρ XY | → 0 then γ XY is ill-defined as the cross-spectrum is identically zeros, and γ XY is thus observed to be estimated with large variance [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF].

The Wavelet estimation framework

It is now well-documented that the analysis of real-world data with scale-free properties can be conducted in a theoretically well grounded and and practically robust and efficient manner using wavelet coefficients, cf. e.g., [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF][START_REF] Abry | Wavelets, spectrum estimation and 1/f processes[END_REF][START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF][START_REF] Abry | Multiscale network traffic analysis, modeling, and inference using wavelets, multifractals, and cascades[END_REF][START_REF] Torres | Comparison of different methods for computing scaling parameter in the presence of trends[END_REF] in general contexts and [START_REF] Ciuciu | Log wavelet leaders cumulant based multifractal analysis of EVI fMRI time series: evidence of scaling in ongoing and evoked brain activity[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF][START_REF] Bullmore | Colored noise and computational inference in neurophysiological (fMRI) time series analysis: resampling methods in time and wavelet domains[END_REF][START_REF] Fadili | Wavelet-generalized least squares: A new BLU estimator of linear regression models with 1/f errors[END_REF] for brain activity analysis. Therefore, the fractal connectivity model proposed here is recast into a wavelet framework. This is detailed in Appendix B, which also discusses scaling parameter estimation.

Rev.#2, Q8 Rev.#2, Q8

Correction for multiple comparisons

Unless otherwise mentioned, all reported p-values for the statistical tests of cor-Rev.#1, Q10 Rev.#1, Q10 relation, self-similarity or fractal connectivity model where Bonferroni corrected for multiple comparisons. To this end, we divided the type-I error rate α = 0.05 by the number of simultaneous comparisons due to the 21 × 20/2 = 210 distinct region pairs. To control the family-wise error rate at level α (equivalent to 1.3 in a -log 10 p val scale as shown in Fig. 6A-B for instance), we computed the corrected p-values by applying the following rule to uncorrected p-values p val-corr = min(1, 210 × p val-uncorr ).

Results

Correlation-based functional connectivity analysis

The fMRI dataset comprised of seventeen right-handed subjects who were scanned at rest and during a visual detection task [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] He | Spontaneous and task-evoked brain activity negatively interact[END_REF][START_REF] Fox | Intrinsic fluctuations within cortical systems account for intertrial variability in human evoked brain responses[END_REF]. FMRI data was preprocessed before extracting signals from 31 ROIs obtained from previous task-related functional neuroimaging studies or generated using coordinates from published fMRI studies, which covered five brain networks [Attention, Default-mode (DMN), Motor, Saliency and Visual] as well as several non-neocortical ROIs (thalamus, cerebellum and hippocampus) (see SI). Projection of regions onto the cortical surface is shown in Fig. 1A. In following analyses, 10 pairs of homologous regions were each averaged together. This step was performed because we observed in preliminary analyses (see also [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF]) that the time series of these homologous regions are highly correlated with each other and that their statistical properties are very similar. Thus, we averaged across homologous reigons in the same RSN in order to enhance statistical independence between investigated brain regions.

Rev.#1, Q13 Rev.#1, Q13
Correlation coefficients ρ XY were estimated from the rest and task dataset separately for all pairs of regions, and were Fisher z-transformed (Z XY ) for statistical testing. Group-level means of correlation coefficients were computed at rest (ρ R XY ) and during task (ρ T XY ) and mapped onto the significant one-sample t-tests ZR XY = 0 and ZT XY = 0, respectively. Comparing Fig. 1B vs Fig. 1C, it can be seen that correlation remains high during task between regions belonging to a same network, whereas lower correlations are observed for between-network region pairs -specifically, between DMN and other networks, a result likely attributable to their respective taskdeactivation and activation [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]. Fig. 1D shows region pairs with a significant change in correlation between rest and task (paired t-test across subjects; H 0 : ZR XY -ZT XY = 0). Significant (p < 0.01) differences between rest and task were observed mainly for pairs Group-averaged inter-regional correlation matrix at rest (p < 0.05, Bonferroni corrected). Regions are grouped by network to ease visualization. Details of the ROIs are provided in Table 1. Middle (C): Group-averaged inter-regional correlation matrix during the visual detection task (p < 0.05, Bonferroni corrected). Bottom (D): Difference of the correlation coefficients between rest and task (thresholded at p < 0.01, uncorrected, two-sample t-test for rest vs. task). The ROIs are grouped by networks whose names are given in Table 1; these networks correspond to the diagonal triangles surrounded by white dashed lines.

of regions consisting of one region located in DMN and the other in the decoupled networks (Attention, Visual, Motor, Saliency and the thalamus). Moreover, in all these region pairs, correlation was higher under rest than task (ρ R XY -ρT XY > 0).

Scale-free univariate analysis

We first applied the wavelet spectrum estimation framework to investigate univariate temporal dynamics of fMRI signals. As an example, Figs. [START_REF] Xiong | Inter-regional connectivity to primary motor cortex revealed using mri resting state images[END_REF][START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]A-B shows the superimposition of the power spectrum estimated by the standard Welch-Periodogram and by wavelet method (log 2 S X (2 j ) vs. log 2 2 j ) for two regions [posterior intraparietal sulcus (pIPS) and middle temporal area (MT)] located in the dorsal attention network (DAN) at rest (Fig. 2A-B) and during task (Fig. 3A-B), respectively. The match between the classical power spectrum and wavelet spectrum confirms that wavelet coefficients can serve as an efficient estimator for the spectrum. Both the classical and wavelet spectra exhibited power-law scaling behavior over the range of 0.01 < f < 0.1 Hz (corresponding to 3.3 < j < 6.6 with f = 2 -j ). Indeed, applying the goodness-of-fit assessment procedure described in Appendix C shows that, out of the 357 analysed time series (17 subjects ×21 regions), rejection of the null hypothesis that the wavelet auto-spectrum is well described by a power law in that range of scales, occurs for only 16% and 14% of cases, at rest and during task respectively. This confirms earlier finding in [42, p. 13788]. Also, similar results were obtained across Rev.#1, Q17

Rev.#1, Q17 Rev.#2, Q5 Rev.#2, Q5
ROIs and subjects, consistent with earlier reports [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF]. Henceforth, this range of frequencies is referred to as the scaling range.

Hurst exponents were then estimated from the wavelet spectrum for each region and individual, separately for the rest (H R ) and task (H T ) dataset. Then, group-level means HR and HT were computed in each ROI to assess the overall effect. When averaged across ROIs within a network (Fig. 4), DMN exhibited the strongest long memory ( HR , HT ) = (0.91, 0.86), followed by the saliency ( HR , HT ) = (0.9, 0.83), attention ( HR , HT ) = (0.9, 0.83) and visual ( HR , HT ) = (0.86, 0.77) networks. The non-neocortical regions ( HR , HT ) = (0.78, 0.73) and the motor network ( HR , HT ) = (0.77, 0.72) exhibited the weakest Hurst exponents. These results are consistent with those previously obtained from the same data using detrended fluctuation analysis (DFA) [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF]. Furthermore, the same network-level ordering was maintained during task while Hurst exponents in all ROIs systematically decreased from rest to task (comparing Fig. 4A and Fig. 4B).

We further estimated Hurst exponents for each ROI and each subject, at rest and during task, using six different estimators: three based on spectral estimation (Direct FFT, Welch-Periodogram and Whittle), two relying on time domain representation (DFA and increments), one constructed on wavelet coefficients. Because it relies on a maximum likelihood principle, Whittle estimator theoretically yields the best estimates for Gaussian data, whereas the wavelet-based estimator has been observed to show significant robustness against additive non-stationary smooth and non-smooth trends. Fig. 5 shows that while the six group-averaged estimates take slightly different values, they are observed to systematically belong to the long memory range 0.5 < H < 1. Moreover, Hurst exponents during task are systematically smaller compared with those at rest, regardless the estimation method: HT < HR . Using paired t-tests, we investigated the statistical significance of this effect and showed that the reduction of self-similarity from rest to task was found significant in the visual network by 3 estimation methods (DFA, increments and wavelets, all results were Bonferroni-corrected for multiple comparisons). The wavelet method additionally uncovered a significant change in the saliency network.These observations confirm, in a robust Rev.#1, Q18 Rev.#1, Q18 manner, earlier reports of long memory at rest and of a decrease in long memory during task [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF]. A decrease in the Hurst exponent during task implies an increased contribution of high frequencies to the temporal dynamics of fMRI signals.

Scale-free multivariate analysis

To assess cross-regional temporal dynamics, we applied wavelet-based multivariate estimation to pairs of regions. For illustration, the Welch-Periodogram estimation of the cross-spectrum between MT and pIPS at rest and during task are illustrated in Figs. [START_REF] Xiong | Inter-regional connectivity to primary motor cortex revealed using mri resting state images[END_REF][START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]C, respectively, and superimposed with the wavelet estimation of the cross-spectrum (log 2 S XY (2 j ) vs. log 2 2 j ). The cross-spectra exhibited powerlaw scaling behavior, within a range of frequencies that matched the scaling range of univariate power spectra. Similar observations were obtained for almost all other region pairs as illustrated in Fig. 6A-B at rest and during task. Indeed, applying the Rev.#1, Q2 Rev.#1, Q2 Goodness-of-Fit assessment procedure described in Appendix C shows that, out of the 3570 = 17 × 21 × 20/2 analysed pairs of regions, the hypothesis that the wavelet cross-spectrum is well described by a power law was rejected only for 14% and 10% of cases, at rest and during task respectively. These observations reveal that scale-free properties are not only characteristics of the univariate fMRI temporal dynamics but also of the multivariate cross-temporal dynamics, within the same scaling range: 0.01 < f < 0.1 Hz. Scale-invariance in univariate temporal dynamics implies that no frequency (in the scaling range) plays a dominant role in the temporal dynamics. Scale-invariance in the multivariate crosstemporal dynamics suggests that synchronization between different brain regions does not rely on a specific frequency, but instead on the intertwining of all frequencies within the scaling range. These findings raise two questions: First, do all frequencies contribute in a balanced manner to inter-regional correlation? Second, does task performance modify scale-free cross-temporal dynamics?

To address these questions, we normalized the cross-spectrum by the auto-power spectra to derive the coherence spectrum. The classical and wavelet-based coherence functions between MT and pIPS for rest and task are shown in Figs. 2-3D, respectively. Because coherence is estimated as the ratio of estimated quantities, it is necessarily noisier. For robustness, we estimate the scaling exponent of the coherence function from the auto-and cross-spectra (γ XY = α XY -(H X + H Y ) + 1, where α XY is the power-law exponent of the cross-spectrum, and H X and H Y are the Hurst exponents of the individual time series), rather than from the coherence function directly.

Using the wavelet-based framework, we estimated exponents α XY and γ XY for all subjects and all region pairs, both at rest and during task. We then computed grouplevel means ᾱXY and γXY . Fig. 6 reports the Bonferroni corrected p-values for the statistical test associated with the null hypothesis ᾱXY = 0, which was rejected for most region-pairs, both at rest (Fig. 6A) and during task (Fig. 6B). This result suggests that cross-temporal dynamics in most region-pairs exhibit a non-zero scaling exponent. Fig. 7A-B reports the Bonferroni corrected p-values for the test against the null hypothesis γXY = 0 at rest and during task respectively. Several conclusions can be reached. First, the null hypothesis γXY = 0 was rejected only for a few region-pairs, 18 at rest and 19 during task, out of 210 pairs (where, due to multiple comparisons, by chance 10 out of 210 may be rejected). Low rejection rate may stem from two reasons: Bonferroni correction for multiple comparisons yields a conservative result; and the statistical power for a test against γXY = 0 has been shown to decrease when |ρ XY | decreases [START_REF] Wendt | Testing fractal connectivity in multivariate long memory processes[END_REF].

However, for region-pairs where the null hypothesis (γ XY = 0) was rejected, a couple of interesting conclusions can be drawn. First, as shown in Fig. 7A-B, whenever departure from 0 was significant, γXY was found to be positive; this was the case both at rest (Fig. 7A) and during task (Fig. 7B). Hence, for these pairs of regions, scale-free properties observed in the cross-spectra convey significant extra information beyond those carried by the auto-spectra: The lowest frequencies within the scaling range contribute substantially more to inter-regional correlation than the highest frequencies, consistent with earlier reports [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Sasai | Frequency-specific functional connectivity in the brain during resting state revealed by NIRS[END_REF][START_REF] Shim | Frequency distribution of causal connectivity in rat sensorimotor network: resting-state fMRI analyses[END_REF].

Second, while a priori ρ XY and γ XY are two independent parameters, it is worth noting that region pairs where γXY = 0 was rejected also tended to show large ρXY both at rest and during task; compare Fig. 1B-C to Fig. 7A-B. In addition, region pairs exhibiting high correlations (see Fig. 8A-B, arbitrary threshold of ρR,T XY > 0.5) systematically showed γR,T XY > 0. These observations indicate that a strong correlation involves dominant contributions of low frequencies to cross-temporal dynamics. Importantly, these observations are not a trivial effect of dependencies, as theoretically the limit ρ XY → 1 imposes γ XY → 0 (i.e., total correlation requires a balanced contribution from all frequencies): This effect can be seen in Fig. 8C-D, where as ρR,T XY approaches 1, γR,T XY approaches 0. To further assess the evolution of γXY from rest to task, a paired t-test across all subjects was applied to every region pair (H 0 : γR XY = γT XY ). Four region pairs showed a significant (p < 0.05, Bonferroni corrected) difference between rest and task (Fig. 9A, R cerebellum -pIPS, hippocampus -R dorsolateral prefrontal cortex, thalamus -L motor cortex and posterior cingulate cortex -superior frontal gyrus). For all significant changes, we observed γR XY > γT XY (Fig. 9B). In 3 out of the 4 regions pairs where the change of γR XY was significant, we observed γT XY > 0 which means that although task significantly modulated γXY and made the contribution of high frequencies to correlation more important, the global correlation remains dominated by the contribution of low frequencies. In contrast, in the R cerebellum -pIPS pair, the switch to γT XY < 0 occurred during task, making the contribution of high frequencies to correlation more significant than that of low frequencies. Altogether, the more balanced frequency-range contribution to correlation during task suggests that cross-dynamics under task involves a larger contribution of the high frequencies, and is driven more by univariate temporal dynamics, since γ T XY is converging to zero. ZT XY and for H 0 : γR XY = γT XY , respectively), we found only one region pair (hippocampus -R dorsolateral prefrontal cortex) showing a significant change in both linear correlation ZXY and coherence scaling exponent γXY , with both quantities decreasing from rest to task.

Scale-free modulation and behavior performance

Lastly, we investigated the relationship between scale-free cross-temporal dynamics, as measured by cross-coherence scaling exponent γ XY and behavioral performance as measured by reaction time (RT). Specifically, we assessed whether across subjects, γ T XY predicts an individual's response speed (mean of RT across trials) and response reliability (s.d. of RT across trials). Table 2A shows that, in 4 region pairs, γ T XY Table 2: Correlation between scale-free cross-temporal dynamics and task performance. A: Acrosssubject correlation between the cross-coherence scaling exponent during task and the s.d. of RT. B: Acrosssubject correlation between difference in cross-coherence scaling exponent between rest and task and the mean of RT. Only significant region-pairs are shown (p < 0.05, FDR-corrected). Group-average values γT XY (A) and γR XY -γT XY (B) are reported in the first column. For full anatomical details of each ROI see Table 1. is significantly correlated with the standard deviation of the recorded RT (p < 0.05, FDR-corrected). In all region pairs, the correlation coefficient was negative, suggesting that the larger the γ T XY , the more reliable the subject's response was across trials. Table 2B further indicates that γ R XY -γ T XY is significantly and positively correlated to the mean RT for three region pairs. Thus, the smaller the difference γ R XY -γ T XY , the shorter the mean RT. By contrast, no significant correlation was found between ρ T XY and the mean or s.d. of RT, after correcting for multiple comparisons.

Group

Recalling that a positive γ T XY indicates that low-frequency range contributes more to correlations than high-frequency range, these observations suggest that a larger contribution of high frequencies to correlation corresponds to poorer behavioral performance. Conversely, the less a subject needs to mobilize high frequencies to accomplish the task, the better his/her performance. Comparing Table 2A and 2B, one can notice the specific role played by the pIPS -medial prefrontal cortex (MPF) pair, a link between the attention and DMN networks, which correlates with both the mean and the standard deviation of the reaction time. In addition, the significant positive correlation between mean RT and the γ R XY -γ T XY measure in the frontal eye field (FEF) -ventral primary visual cortex (vRetino) pathway indicates that the stability of the scale-free cross-temporal dynamics in this region pair predicts the speed of task execution.

Altogether, these observations indicate that the extent to which task modulates inter-regional correlation and the balance between different frequencies' contribution to global correlation is negatively correlated with behavioral performance. Thus, across subjects, the less the temporal dynamics of cross-correlation are altered during task, the better the performance.

Discussion

Using multivariate analyses of fMRI signals within the framework of scale-free dynamics, the present work sheds light on several characteristics of brain temporal dynamics.

First, scale-invariance is an important aspect of brain dynamics, being observed not only in univariate analysis of each region independently, but also in the crossregional temporal dynamics. This observation links functional connectivity and scalefree dynamics in fMRI signals -two hitherto separately studied topics. In terms of neurophysiological interpretation, our findings have shown that the communication between distant brain regions, which is captured by the cross-spectrum, is brought by information/energy exchange over a range of frequencies in a scale invariant manner, without any dominant frequency between 0.01 and 0.1 Hz.

Rev.#1, Q1, Q7 Rev.#1, Q1, Q7
Second, the scale-free cross-temporal dynamics do not follow the fractal connectivity model in many region-pairs (Fig. 7), suggesting that: i) fMRI signals are not well described by the simple and direct multivariate extension of the classical fGntype univariate models; ii) estimation of the cross-spectrum brings extra information beyond that contained in the auto-spectra; iii) low frequencies (close to 0.01 Hz) contribute more to functional connectivity, and thus to communication between regions, than high frequencies (close to 0.1 Hz). Several previous studies have also reported Rev.#1, Q1, Q7 Rev.#1, Q1, Q7 a stronger contribution of low frequencies to inter-regional correlation [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Sasai | Frequency-specific functional connectivity in the brain during resting state revealed by NIRS[END_REF][START_REF] Shim | Frequency distribution of causal connectivity in rat sensorimotor network: resting-state fMRI analyses[END_REF] by comparing the power in different frequency bands, yet without investigating the frequency-range balance quantified by scaling parameters. Our study thus significantly Rev.#2, Q2 Rev.#2, Q2 extends these previous reports by exposing the scale-invariance thereof and demonstrating a link between functional connectivity and scale-free dynamics. In particular, we observed that the scaling exponent measured from the cross-spectrum does not simply consist of an average of those estimated on auto-spectra, since the coherence scaling exponent (γ XY ) departs from zero. Also, we systematically observed that this deviation occurs in the same direction: γ XY > 0, hence cross-temporal dynamics are driven by low frequencies. On a technical note, the lack of scale-invariance in prior re-Rev.#2, Q2 Rev.#2, Q2 ports was likely due to a methodological difference from our study, as the methods used there were not tailored to investigating long-range temporal behavior (e.g., Ref. [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF] focused on the time-evolution of inter-regional correlations and Ref. [START_REF] Shim | Frequency distribution of causal connectivity in rat sensorimotor network: resting-state fMRI analyses[END_REF] employed an autoregressive model which accounted only for short-range correlations).

Third, in addition to a decrease in linear correlation and Hurst exponent reported in previous work [START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF], task induces a decrease of the cross-coherence scaling exponent γ XY . In other words, cross-temporal dynamics are closer to fractal connectivity under task. From a neurophysiological perspective, this means that the temporal dynamics of communication between brain regions are altered under task performance in terms of frequency content.In additoin, these results are consistent with a previous Rev.#1, Q7 Rev.#1, Q7 study showing that increased attention induces a decrease of coherence between neuronal populations in macaques specifically in the low frequencies [START_REF] Mitchell | Spatial attention decorrelates intrinsic activity fluctuations in macaque area V4[END_REF]Fig. 4]. Note that a priori the variations of these three parameters are independent; thus, the finding that they occur jointly in brain activity is nontrivial. These findings indicate that the decrease in correlation is accompanied by a stronger mobilization of high frequencies within the scaling range in both the univariate temporal dynamics and in the multivariate cross-temporal dynamics, the latter corresponding to a more balanced contribution of all frequencies to correlation (i.e., functional connectivity).

The current task was a very simple visual detection task, which engaged the visual, motor, saliency and attention systems as shown for instance in [42, Fig. 6, p. 13192]. This simple task already impacted functional connectivity between task-positive (attention, motor, saliency, visual) and task-negative (default mode) networks, as illustrated in Fig. 1D. The use of Bonferroni correction made these results very specific and reliable. The recourse to an alternative approach like FDR for addressing multiple comparisons would probably have shown larger functional connectivity differences between rest and task. Nonetheless, more demanding cognitive tasks, especially those with a learning component, may demonstrate a larger change in functional connectivity [START_REF] Lewis | Learning sculpts the spontaneous activity of the resting human brain[END_REF]. In the same spirit but focusing on scaling exponents instead of correlation measures, other contributions [START_REF] Zilber | Modulation of scale-free properties of brain activity in MEG[END_REF][START_REF] Zilber | Learning-induced modulation of scalefree properties of brain activity measured with MEG[END_REF] have shown that the multifractal properties of MEG source reconstructed time series continuously evolve with perceptual learning in the task-related networks associated with a visual discrimination task.

Rev.#1, Q3, Q14 Rev.#1, Q3, Q14 Fourth, we observed that across subjects, a larger increase of the high-frequency contribution to cross-temporal dynamics under task was associated with worse behavioral performance. Hence, a strong modulation of cross-temporal dynamics may indicate difficulty in performing the task, consistent with the idea that ongoing fluctuations captured by low-frequency functional connectivity are important for behavioral processes [START_REF] Sadaghiani | The relation of ongoing brain activity, evoked neural responses, and cognition[END_REF]. However, we also outlined that these findings may results from attention effects, as originally observed in macaques [START_REF] Mitchell | Spatial attention decorrelates intrinsic activity fluctuations in macaque area V4[END_REF]. To further investigate such issues and disentangle attention from operative effects in the recourse to high frequencies, future work will be devoted to the analysis of another existing MEG dataset [START_REF] Zilber | Supramodal plasticity optimizes perceptual learning[END_REF] for which complementary eye tracker recordings will permit to probe attention through measurements of ocular saccades in cunjonction with behavioral performance.

Rev.#2, Q3 Rev.#2, Q3
On a methodological note, we made use of 28 minutes of resting-state fMRI acquisition which is a relatively large amount of data compared to a typical resting-state fMRI experiments. However, the dataset was split in 4 alternating blocks of resting-state and task-related scans of 7 min each. Thus, we computed the scaling parameter estimates for each block individually, and then averaged the results over the 4 blocks for each condition (rest or task). This averaging increased the robustness of our analysis, compared with a single block analysis.

Rev.#1, Q8 Rev.#1, Q8

The use of the wavelet framework in the present work allows significant robustness in estimating univariate or multivariate scale-free dynamics in fMRI data, especially with respect to short sample size or the presence of slow superimposed trends [START_REF] Wendt | Testing fractal connectivity in multivariate long memory processes[END_REF][START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF][START_REF] Torres | Comparison of different methods for computing scaling parameter in the presence of trends[END_REF]. In the present study, results were averaged across rest and task runs, respectively, which assumes that they were stable across time. This might hide a source of variability as suggested by recent studies on fMRI dynamic functional connectivity [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF][START_REF] Hutchison | Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques[END_REF][START_REF] Majeed | Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans[END_REF][START_REF] Smith | Temporally-independent functional modes of spontaneous brain activity[END_REF][START_REF] Kiviniemi | A sliding time-window ica reveals spatial variability of the default mode network in time[END_REF][START_REF] Leonardi | Principal components of functional connectivity: A new approach to study dynamic brain connectivity during rest[END_REF][START_REF] Liu | Time-varying functional network information extracted from brief instances of spontaneous brain activity[END_REF] and a previous MEG study showing that the amount of self-similarity might change over time [START_REF] Zilber | Learning-induced modulation of scalefree properties of brain activity measured with MEG[END_REF]. To address this issue, it would be informative for future studies to make use of high temporal resolution fMRI [START_REF] Boyacioglu | Generalized iNverse imaging (GIN): Ultrafast fMRI with physiological noise correction[END_REF][START_REF] Feinberg | Multiplexed echo planar imaging for sub-second whole brain FMRI and fast diffusion imaging[END_REF][START_REF] Moeller | Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI[END_REF][START_REF] Rabrait | High temporal resolution functional MRI using parallel echo volume imaging[END_REF], and to explore scale-free cross-temporal dynamics examined herein in a time-dependent manner.

Recent studies have also shown that functional connectivity can change over time [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF][START_REF] Hutchison | Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques[END_REF][START_REF] Majeed | Spatiotemporal dynamics of low frequency BOLD fluctuations in rats and humans[END_REF][START_REF] Smith | Temporally-independent functional modes of spontaneous brain activity[END_REF][START_REF] Kiviniemi | A sliding time-window ica reveals spatial variability of the default mode network in time[END_REF][START_REF] Leonardi | Principal components of functional connectivity: A new approach to study dynamic brain connectivity during rest[END_REF][START_REF] Liu | Time-varying functional network information extracted from brief instances of spontaneous brain activity[END_REF] suggesting the presence of nonstationarity in fMRI correlation structures. Future work should investigate whether the frequency-range balance of inter-regional connectivity is changing concomitantly at critical time points where functional connectivity is reconfigured (e.g., [START_REF] Liu | Time-varying functional network information extracted from brief instances of spontaneous brain activity[END_REF]) or whether these phenomena occur independently. This issue of nonstationarity further points out the need for continuation of the present work beyond a second-order stationary framework. First, potential time evolutions should be investigated (following, e.g., approach in [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]), yet tailored to scale-free dynamics. Second, extensions of this multivariate framework towards scalefree dynamics at higher statistical orders (referred to as multifractal properties) in fMRI signals [see Refs [START_REF] Ciuciu | Log wavelet leaders cumulant based multifractal analysis of EVI fMRI time series: evidence of scaling in ongoing and evoked brain activity[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF] for univariate applications] should be investigated. This effort should help elucidating whether dependencies beyond correlation and second-order statistics play an active role in the dynamical reconfiguration of functional connectivity.

Conclusion

Rev.#2, Q1 Rev.#2, Q1

In conclusion, by showing that scale-free temporal dynamics manifest in the communication between brain regions, our results provide a bridge between two related, but so-far separate, fields -resting-state networks and scale-free dynamics, which have respectively studied spontaneous brain activity in the spatial and temporal domain. In particular, we observed that the lowest frequencies contributed more to inter-regional communication under both rest and task, but interestingly, this effect was ameliorated under task performance, with different frequencies contributing more equally to interregional correlation. Furthermore, we found that the degree to which task performance modulated the scaling behavior of cross-regional temporal dynamics was correlated, across subjects, with behavioral performance, such that smaller task modulation was accompanied by faster and more consistent reaction times. These results should inspire future studies of the interplay between scale-free brain dynamics and large-scale brain networks. where the R x , R y and R z matrices match the pitch, roll and yaw movements, respectively. Of course, composing rotations is not a commutative operation so that any alternative composition will deliver distinct results. For each TR in each rest and task run, we computed the corresponding R G using Eq. (A.1). Hence, we extracted its spectral norm as follows:

R G 2 = λ max (R t G R G ) (A.2)
where λ max (A) stands for the largest eigenvalue of matrix A. Since R t G R G is positive semi-definite, all its eigenvalues are bounded below by zero. Clearly, Eq. (A.2) summarizes in a single scalar the largest rotation direction. By computing this value for all TRs and then extracting the meanµ r and standard deviation (σ r ) over all volumes in a given run, we can therefore define the associated thresholds α 1 = µ r + / -2σ r , α 2 = µ r + / -3σ r for detecting rotation movement spikes. As done before for the translation movement, we repeated this procedure over all rest and task-related runs for each individual so as to average the number of spikes per individual over the complete fMRI session. The grand average number of rotation spikes was eventually computed by averaging over the 17 subjects who underwent the study. The outcomes of our spike detection procedure wrt rotation movement are summarized in Table 3[Bottom]. They show that the presence of spikes is even more negligible compared to what we found for translation (1.4% wrt α 1 threshold and less than 1% wrt α 2 threshold both in rest and task related runs).

Appendix B. The Wavelet estimation framework

Therefore, the wavelet-based coherence function can be estimated as:

CW XY (2 j ) = S XY (2 j ) S X (2 j )S Y (2 j ) . (B.10)
When X, Y follow the bivariate model defined in Eq. ( 2), one obtains: S X (2 j ) ω X 2 j α X , S Y (2 j ) ω Y 2 j α Y , S XY (2 j ) ω XY 2 j α XY and CW XY (2 j ) γ 0 2 j γ XY .

Fig. 10, for synthetic data, and Figs. 23, for the real data analyzed here, display log S X , log S Y , log S XY and log CW XY as functions of log 2 2 j = j and thus illustrate the corresponding power law behaviors for the auto and cross (wavelet) spectra and (wavelet) cross coherence. Following [START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF] or [START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF], estimation of the scaling parameters α X (hence Ĥx ), α Y (hence Ĥy ) and α XY stems from linear regressions performed in these log-coordinate plots, across the scaling range (j 1 = log 2 f 0 /f M ; j 2 = log 2 f 0 /f m , tuned to match the frequency range f m |f | f M . Further, the estimate of γ XY is obtained as γ XY = α XY -( Ĥx + Ĥy ) + 1, and not as a linear regression of log CW XY against log 2 2 j = j, cf. [START_REF] Wendt | Testing fractal connectivity in multivariate long memory processes[END_REF] for details.

Varying N ψ , this wavelet framework provides practitioners with an efficient and robust tool to estimate γ XY and γ 0 (and thus ρ XY ) on real-world data and is systematically used in the present work to produce the results.

Appendix B.4. Illustration on bivariate fractional Gaussian noise

The wavelet estimation framework is illustrated by application to synthetic bivariate fractional Gaussian noise, synthesized using the theoretical procedure and practical codes devised in [START_REF] Helgason | Fast and exact synthesis of stationary multivariate Gaussian time series using circulant embedding[END_REF], [START_REF] Helgason | Synthesis of multivariate stationary series with prescribed marginal distributions and covariance using circulant matrix embedding[END_REF]. Parameters for synthetic data are chosen to match as closely as can be those estimated from the real data used to produce Fig. 2, and with same sample size. It permits to observe that estimated scaling from real data, in Fig. 2, are visually as convincing as those in Fig. 10 obtained from synthetic bivariate fractional Gaussian noise, with same parameters, known theoretically to have true scaling behaviors. Scaling relevance is further assessed by means of statistical tests as described in Appendix C below. Gaussian noise (mfGn) process synthetised with the same scaling exponent γ XY = as the one estimated on real fMRI data at rest; see Fig. 2. A-B: Superimposition of the auto-power spectrum estimated by Welch-Periodogram (black) and wavelet method (red) for two of its components X and Y . C: Superimposition of the cross-power spectrum between X and Y estimated by Welch (black) and wavelet method (red) in log-log coordinates. D: Superimposition of the coherence function between X and Y estimated by Welch (black) and wavelet (red) method. All exponents are estimated from linear regression based on the wavelet estimate (dashed red lines) in the scaling range defined by vertical dashed black lines (-6.6 < log 2 f < -3.3) -1000 copies of synthetic bivariate fractional Gaussian noises, with same parameters as data, were simulated (according to the procedure theoretically devised in [START_REF] Helgason | Fast and exact synthesis of stationary multivariate Gaussian time series using circulant embedding[END_REF][START_REF] Helgason | Synthesis of multivariate stationary series with prescribed marginal distributions and covariance using circulant matrix embedding[END_REF] ) and then for each copy, we performed estimation of H X , H Y , α XY , and computed the classical χ 2 goodness-of-fit quantity for both the auto-wavelet spectra and the cross-wavelet spectra;

-the p-value, corresponding to the test aiming at rejecting the null hypothesis that true data have same auto-and cross-wavelet spectra as bivariate fractional Gaussian noise, with same parameters and same sample size, is computed as the percentage of times the χ 2 goodness-of-fit value computed from synthetic data exceeds that computed from real data;

-the null hypothesis is rejected at the 0.05 level.

Figure 1 :

 1 Figure 1: Networks definition and correlation structure. Top (A): ROIs mapped onto the cortical surface, with each color denoting a different network. Middle (B):Group-averaged inter-regional correlation matrix at rest (p < 0.05, Bonferroni corrected). Regions are grouped by network to ease visualization. Details of the ROIs are provided in Table1. Middle (C): Group-averaged inter-regional correlation matrix during the visual detection task (p < 0.05, Bonferroni corrected). Bottom (D): Difference of the correlation coefficients between rest and task (thresholded at p < 0.01, uncorrected, two-sample t-test for rest vs. task). The ROIs are grouped by networks whose names are given in Table1; these networks correspond to the diagonal triangles surrounded by white dashed lines.

Figure 2 :

 2 Figure 2: Univariate (A-B) and multivariate (C-D) scale-free properties of fMRI signals at rest -example from a single subject. A-B: Superimposition of the auto-power spectrum estimated by Welch-Periodogram (black) and wavelet method (red) for pIPS and MT. C: Superimposition of the cross-power spectrum between pIPS and MT estimated by Welch (black) and wavelet method (red) in log-log coordinates. D: Superimposition of the coherence function between pIPS and MT estimated by Welch (black) and wavelet (red) method. All exponents are estimated from linear regression based on the wavelet estimate (dashed red lines) in the scaling range defined by vertical dashed black lines (-6.6 < log 2 f < -3.3).

  Rev.#2, Q5 Rev.#2, Q5 Rev.#1, Q17 Rev.#1, Q17

Figure 3 :

 3 Figure 3: Univariate (A-B) and multivariate (C-D) scale-free properties of fMRI signals during task (same subject as in Fig. 2. A-B: Superimposition of the auto-power spectrum estimated by Welch-Periodogram (black) and wavelet method (red) for pIPS and MT. C: Superimposition of the cross-power spectrum between pIPS and MT estimated by Welch (black) and wavelet method (red) in log-log coordinates. D: Superimposition of the coherence function between pIPS and MT estimated by Welch (black) and wavelet (red) method. All exponents are estimated from linear regression based on the wavelet estimate (dashed red lines) in the scaling range defined by vertical dashed black lines (-6.6 < log 2 f < -3.3).

Figure 4 :

 4 Figure 4: Hurst exponents from scale-free univariate analysis. Group-averaged region-wise (bottom row) and network-average (top row) wavelet-based estimates of Hurst exponent at rest (A) and during task (B). The color scale is the same for A and B, to illustrate the decrease of Hurst exponents from rest to task.

Figure 5 :

 5 Figure 5: Group-level network-average Hurst exponent at rest (black) and during task (red) using six different estimation methods. A: Standard periodogram (squared fast Fourier transform normalized by the number of samples) and subsequent linear regression on the log-log power spectrum plot. B: Welchbased (i.e., averaged across overlapping windows) periodogram followed by linear regression on the loglog power spectrum plot. C: Whittle estimator, which consisted of a maximum likelihood estimator of power spectrum under fractional Gaussian noise model. D: Detrended Fluctuation Analysis (DFA). E: Timedomain increment-based Hurst exponent estimate. F: Wavelet spectrum estimate, where the Hurst exponent was estimated from a linear regression in the log-log plot where the log along the x-axis involves scales instead of frequencies. X-axis labels indicate networks: A, Attention; D, Default-mode; M, Motor; N, Nonneocortical; S, Saliency; V, Visual. The shaded areas outline the significant differences of Hurst exponents between rest and task. Significant differences are indicated by *-marks (p-val<0.05, Bonferroni corrected).

Figure 6 :

 6 Figure 6: Cross-spectrum scaling exponents. Group-average values of the scaling exponent of the crossspectrum at rest (A) and during task (B), i.e. ᾱRXY and ᾱT XY , respectively. Only region pairs where the scaling exponent significantly departed from 0 (p < 0.05, Bonferroni corrected, one-sample t-test) are shown. Thus, scaling is significant for all cross-spectra associated with "hot squares".

Figure 7 :

 7 Figure 7: Cross-coherence scaling exponents. Group-average values of the scaling exponent of the crosscoherence function at rest (A) and during task (B), i.e. γR XY and γT XY , respectively. Only region pairs where the scaling exponent significantly departed from 0 (p < 0.05, Bonferroni corrected, one-sample ttest) are shown. So, scale invariance on the cross-temporal dynamics has a more complex structure than just averaging those coming from the ROI-based univariate time series in all "hot squares".

Figure 8 :

 8 Figure 8: Relationship between the cross-coherence scaling exponent and the linear Pearson correlation. A-B: For region pairs with a relatively large Pearson correlation (ρ XY > 0.5), the scaling exponent of the cross-coherence function is plotted as color, for rest (A) and task (B), respectively. C-D: Cross-coherence scaling exponent plotted against the linear correlation coefficient for all region-pairs at rest (C) and during task (D). Within-and between-network region pairs are shown in red and blue, respectively.

Figure 9 :

 9 Figure 9: Comparison of cross-coherence scaling exponent between rest and task. A: Significant differences in cross-coherence scaling exponent between rest and task (p < 0.05, Bonferroni corrected). B: For the four significant region-pairs in A, cross-coherence scaling exponent is shown for rest and task, respectively. Only one region-pair was within the same network (PCC -SFG, within the DMN). C: Grand average (with +/-standard deviations of the mean) of wavelet-based coherence functions between the Thalamus and Lmotor regions at rest (black trace) and during task (red trace).

Figure 10 :

 10 Figure 10: Univariate (A-B) and multivariate (C-D) scale-free properties of a multivariate fractional Gaussian noise (mfGn) process synthetised with the same scaling exponent γ XY = as the one estimated on real fMRI data at rest; see Fig.2. A-B: Superimposition of the auto-power spectrum estimated by Welch-Periodogram (black) and wavelet method (red) for two of its components X and Y . C: Superimposition of the cross-power spectrum between X and Y estimated by Welch (black) and wavelet method (red) in log-log coordinates. D: Superimposition of the coherence function between X and Y estimated by Welch (black) and wavelet (red) method. All exponents are estimated from linear regression based on the wavelet estimate (dashed red lines) in the scaling range defined by vertical dashed black lines (-6.6 < log 2 f < -3.3)

  ∀f , all frequencies are equally or equivalently contributing to global correlation. Conversely, frequencies such that |C XY (f )| > |ρ XY | contributes more to global correlation compared to frequencies where |C XY (f )| < |ρ XY |.

Table 3 :

 3 Results of our spike detection procedure with respect to translation and rotation movements. Grand average results are reported separately for rest and task runs after averaging first over all runs at the subject-level and then over all the individuals.

		Rest		Task
		% spikes (α 1 ) % spikes (α 2 ) % spikes (α 1 ) % spikes (α 2 )
	Translation	3.6%	0.5%	3.4%	0.8%
	Rotation	1.4%	0.7%	1.4%	0.6%

http://brainvis.wustl.edu/wiki/index.php/Caret:About

Rev.#1, Q4 Rev.#1, Q4 Rev.#2, Q4 Rev.#2, Q4

Appendix A. Spike detection procedure

We analyzed the movement parameter estimates by looking at translation and rotation separately. As regards translation, for each TR, run and individual, we computed the vector norm defined by the 3 translation parameters. Then, we computed the mean translation (µ t ) by averaging over all TRs in a given run. We extracted similarly the corresponding standard deviation (σ t ). Hereafter, we identified the number of TRs for which the translation movement exceeds α 1 = µ t + / -2σ t and α 2 = µ t + / -3σ t . We repeated this procedure over all rest and task-related runs for each individual so as to average the number of spikes per individual over the complete fMRI session. The grand average number of translation spikes was eventually computed by averaging over the 17 subjects who underwent the study. The outcomes of our spike detection procedure with respect to (wrt) translation movement are summarized in Table 3[Top]. They show that the presence of spikes is very negligible (less than 1% wrt α 2 threshold both in rest and task related runs).

As regards rotation, the problem is more complex. The difficulty lies in how to collapse the 3 rotation parameters in a single quantity describing a global 3D rotation. We decided to compute the global rotation matrix R G as follows:

Let ψ 0 (t) denote a reference oscillating function with narrow supports in both time and frequency domains, referred to as the mother wavelet. It is characterized by its number of vanishing moment, a strictly positive integer N ψ defined as:

The discrete wavelet transform (DWT) coefficients of X are defined as:

Scale 2 j qualitatively corresponds to the inverse of the frequency, 2 j ∼ f 0 /f , where f 0 is a constant that depends on the choice of ψ 0 (t). For further details, readers are referred to e.g., [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Appendix B.2. Wavelet coherence function

Let X, Y denote a bivariate second order stationary processes. It has been shown that [START_REF] Abry | Wavelets, spectrum estimation and 1/f processes[END_REF][START_REF] Abry | Multiscale network traffic analysis, modeling, and inference using wavelets, multifractals, and cascades[END_REF]:

where Ψ 0 stands for the Fourier transform of ψ 0 and E for the mathematical expectation. The quantities Ed X (j, k) 2 , Ed Y (j, k) 2 and Ed X (j, k)d Y (j, k) can thus be read as the auto-and cross-wavelet spectra, measuring the frequency contents of data around frequency f = f 0 2 -j . Following [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF], a wavelet-based coherence function can now be introduced as:

When X, Y follow the bivariate model defined in Eq. ( 2), it yields:

Appendix B.3. Estimation procedure Following [START_REF] Abry | Wavelets, spectrum estimation and 1/f processes[END_REF][START_REF] Abry | Multiscale network traffic analysis, modeling, and inference using wavelets, multifractals, and cascades[END_REF][START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF], relevant estimators for the auto-and cross-wavelet spectra can be defined as time averages of the (squared) wavelet coefficients at scale 2 j :

Appendix C. Goodness-of-Fit test for multivariate scaling Following the methodologies outlined in [START_REF] Clauset | Power law distributions in empirical data[END_REF][START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF], we have implemented the following Goodness-of-Fit assessment procedure for each subject and each pair of regions:

Rev.#2, Q5 Rev.#2, Q5 -estimation of the scaling and correlation parameters (H X , H Y , α XY and ρ) is performed;

-computation of the classical χ 2 (sum of squared errors to the best fitted linear model) goodness-of-fit quantities for both the auto-and cross wavelet spectra is done;