Daria La Rocca

Herwig Wendt

Virginie Van Wassenhove

Philippe Ciuciu

Patrice Abry
email: patrice.abry@ens-lyon.fr

Revisiting Functional Connectivity for Infraslow Scale-Free Brain Dynamics using Complex Wavelets

Keywords: Human brain temporal dynamics, functional connectivity, infraslow, arrhythmic, scale-free, phase coupling, functional connectivity assessed from fractal dynamics, complex-wavelet, MEG data

The analysis of human brain functional networks is achieved by computing functional connectivity indices reflecting phase coupling and interactions between remote brain regions. In magneto-and electroencephalography, the most often used functional connectivity indices are constructed on Fourier-based cross spectral estimation applied to specific fast and band limited oscillatory regimes. Recently, infraslow arrhythmic fluctuations (below the 1Hz) were recognized as playing a leading role in spontaneous brain activity. The present work aims to propose to assess functional connectivity, from fractal dynamics, thus extending the assessment of functional connectivity to the infraslow arrhythmic or scale-free temporal dynamics of M/EEG-quantified brain activity. Instead of being based on Fourier analysis, new Imaginary Coherence and weighted Phase Lag indices are constructed from complex-wavelet representations. Their performance are, first, assessed on synthetic data, by means of Monte-Carlo simulations, and compared favorably against the classical Fourier-based indices. These new assessment of functional connectivity indices are, second, applied to MEG data collected on 36 individuals, both at rest and during the learning of a visual motion discrimination 1 La Rocca et al.

Functional Connectivity assessed from fractal dynamics

task. They demonstrate a higher statistical sensitivity, compared to their Fourier counterparts, in capturing significant and relevant functional interactions in the infraslow regime, and modulations from rest to task. Notably, the consistent overall increase in functional connectivity assessed from fractal dynamics from rest to task, correlated with a change in temporal dynamics, as well as with improved performance in task completion, suggests that complex-wavelet weighted Phase Lag index is the sole index able to capture brain plasticity in the infraslow scale-free regime.

INTRODUCTION

Human brain univariate temporal dynamics. The dynamics of Human brain activity can be studied non-invasively using electro-and magnetoencephalography (EEG and MEG, respectively). Interpreted as resulting from the synchronous activation of neuronal populations in specific frequency bands, these fluctuations are often analyzed as fast (10Hz and above) oscillatory rhythms, now well associated with cognitive functions, such as perception, attention or decision making (cf. e.g., [START_REF] Freeman | Mesoscopic neurodynamics: from neuron to brain[END_REF][START_REF] Jensen | Cross-frequency coupling between neuronal oscillations[END_REF])), well-described by band-limited models and well-analyzed by classical Fourier transform-based spectral analysis.

At the turn of the 21st century, the large-band infraslow activity of the brain (typically below 1 Hz), that for long had been considered as either instrumental or head-movement noises, received growing interest and has been documented as a prominent part of recorded electromagnetic brain signals and a critical component of brain activity [START_REF] Gong | Scale-invariant fluctuations of the dynamical synchronization in human brain electrical activity[END_REF][START_REF] Vanhatalo | Infraslow oscillations modulate excitability and interictal epileptic activity in the human cortex during sleep[END_REF][START_REF] Stam | Scale-free dynamics of global functional connectivity in the human brain[END_REF][START_REF] Miller | Power-law scaling in the brain surface electric potential[END_REF][START_REF] Werner | Fractals in the nervous system: conceptual implications for theoretical neuroscience[END_REF]. This large-band infraslow activity in the brain differs significantly from band-limited oscillations in the sense that it is not characterized by specific frequencies or scales of times, but rather corresponds to arrhythmic, or scale-free, temporal dynamics. While exact scale-free dynamics remains debatable [START_REF] Dehghani | Comparative power spectral analysis of simultaneous elecroencephalographic and magnetoencephalographic recordings in humans suggests non-resistive extracellular media[END_REF][START_REF] Ignaccolo | Dynamics of electroencephalogram entropy and pitfalls of scaling detection[END_REF], it has been proposed by an abundant literature (cf. eg., [START_REF] Vanhatalo | Infraslow oscillations modulate excitability and interictal epileptic activity in the human cortex during sleep[END_REF][START_REF] He | The temporal structures and functional significance of scale-free brain activity[END_REF][START_REF] Van De Ville | EEG microstate sequences in healthy humans at rest reveal scale-free dynamics[END_REF][START_REF] Dehghani | Comparative power spectral analysis of simultaneous elecroencephalographic and magnetoencephalographic recordings in humans suggests non-resistive extracellular media[END_REF][START_REF] He | Scale-free properties of the functional magnetic resonance imaging signal during rest and task[END_REF][START_REF] He | Scale-free brain activity: past, present, and future[END_REF][START_REF] Buzsáki | The log-dynamic brain: how skewed distributions affect network operations[END_REF][START_REF] Zilber | Modulation of scale-free properties of brain activity in MEG[END_REF][START_REF] Gadhoumi | Scale invariance properties of intracerebral EEG improve seizure prediction in mesial temporal lobe epilepsy[END_REF][START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF]) that infraslow macroscopic brain activity is better described by a scaling exponent (historically the power-law exponent of the Fourier spectrum, and more recently and relevantly the selfsimilarity exponent H) that relates together dynamics across a large continuum of slow time scales (or low frequencies).

While most oscillatory regimes are only observed in evoked activity, elicited by stimuli, infraslow scale-free brain temporal dynamics are persistent, observed both at rest and during task performance, or even in unconscious states (e.g., sleep stages). It was also shown that infraslow scale-free brain temporal dynamics are modulated when contrasting rest and task-related brain activity, with task inducing systematically a decrease in H and faster infraslow dynamics [START_REF] Bhattacharya | Universality in the brain while listening to music[END_REF][START_REF] Linkenkaer-Hansen | Stimulus-induced change in long-range temporal correlations and scaling behaviour of sensorimotor oscillations[END_REF][START_REF] Vanhatalo | Infraslow oscillations modulate excitability and interictal epileptic activity in the human cortex during sleep[END_REF][START_REF] Popivanov | Multifractality of decomposed EEG during imaginary and real visual-motor tracking[END_REF][START_REF] Bianco | Brain, music, and non-poisson renewal processes[END_REF][START_REF] Buiatti | Feedback modulates the temporal scale-free dynamics of brain electrical activity in a hypothesis testing task[END_REF][START_REF] He | The temporal structures and functional significance of scale-free brain activity[END_REF][START_REF] Zilber | Learning-induced modulation of scale-free properties of brain activity measured with MEG[END_REF][START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF]. Infraslow scale-free brain activity has thus been hypothesized to be functionally associated with neural excitability [START_REF] He | Scale-free brain activity: past, present, and future[END_REF].

Altered scale-free brain dynamics has also been reported in specific condition such as Alzheimer's disease for which larger selfsimilarity exponents were reported in multiple brain areas (e.g. lateral temporal lobes, insula, etc) early involved in the neurodegenerative process [START_REF] Maxim | Fractional gaussian noise, functional mri and alzheimer's disease[END_REF].

Infraslow arrhythmic brain activity can be efficiently described with large-band scale-free models, such as selfsimilarprocesses (fractional Brownian motion and fractional Gaussian noise) [START_REF] Mandelbrot | Fractional Brownian motion, fractional noises and applications[END_REF]. It is also now well established and documented that, while Fourier analysis can be used to assess 1/f power law spectra at low frequencies, accurate and robust assessments of scale-free dynamics requires replacing Fourier-based spectral estimation with multiscale wavelet analysis. Interested readers are referred to [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motions[END_REF][START_REF] Muzy | Multifractal formalism for fractal signals: The structure-function approach versus the wavelet-transform modulus-maxima method[END_REF][START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF][START_REF] Kantelhardt | Fractal and multifractal time series[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF] for methodological developments, and to [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] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks[END_REF][START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF] for applications to neuroimaging data. Further, it has recently been shown that the dself-similar escription of scale-free temporal dynamics could be enriched by combining the concept of multifractality with that of selfsimilarity [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF], requiring the use of wavelet-leaders, consisting of nonlinear nonlocal transforms of wavelet coefficients, for practical analysis. The potential interest of multifractality for the analysis of fMRI and M/EEG signals has been investigated in e.g., [START_REF] Shimizu | Waveletbased multifractal analysis of fMRI time series[END_REF][START_REF] Popivanov | Effect of independent component analysis on multifractality of EEG during visual-motor task[END_REF][START_REF] Popivanov | Multifractality of decomposed EEG during imaginary and real visual-motor tracking[END_REF][START_REF] Shimizu | Neuronal response to Shepard's tones: an auditory fmri study using multifractal analysis[END_REF][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] Proekt | Scale invariance in the dynamics of spontaneous behavior[END_REF][START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF].

Human brain multivariate temporal dynamics: Functional connectivity. Remote brain regions are known to interact within large scale functional networks (e.g., the default Mode Network at rest) which mediate the information flow inside the brain integrating the activity of functionally segregated modules that are activated in particular mental states, task execution or health condition [START_REF] Power | Functional network organization of the human brain[END_REF]. These interactions (correlations, delays, phase synchronization,. . . ) between different brain regions are quantified by indices of similarity computed from signals collected in each region, and are referred to as functional connectivity. Assessing functional connectivity thus entails performing a multivariate analysis of the temporal recordings, thus complementing univariate analysis of each signal Frontiers separately. Classically, functional connectivity is assessed mostly from band-limited signals reflecting the oscillatory activity of the brain, by measures of cross (bivariate) second order statistics (correlation coefficient, cross-correlation function,. . . ). However, M/EEG measurements suffer from the so-called volume conduction effects: Linearity in Maxwell equations and electromagnetic quasi-static approximation (for the forward model below 100 Hz) induce a linear mixing of electromagnetic sources on M/EEG sensors with negligible temporal delays. Close-by EEG electrodes or SQUID MEG sensors thus redundantly capture brain activity from a given current cortical dipole hence inducing spurious correlations amongst recordings, thus precluding a relevant assessment of functional connectivity [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF][START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF][START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF]. Source-space reconstructed signals are documented to still suffer from residual volume conduction effects because of the approximate and imperfect nature of inverse problem resolutions [START_REF] Siebenhühner | Measuring large-scale synchronization with human MEG and EEG: challenges and solutions[END_REF][START_REF] Palva | Ghost interactions in MEG/EEG source space: A note of caution on inter-areal coupling measures[END_REF]. The design of indices robust to such spurious correlations has been based on measuring average phase delays, such as in the Phase Locking Value [START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF], and also naturally calls for the use of Fourier-based cross spectral estimation. Indeed, the Fourier transform, being by definition based on complex numbers, permits to incorporate automatically phases and thus delays in the assessment of functional connectivity: zero delay between correlated signals corresponds to zero phase and imaginary part but nonzero real part for the cross-Fourier spectrum (on average). Therefore, the moduli of the cross-Fourier spectrum and of the coherence function (F-COH) are affected by volume conduction effects, but their imaginary parts and phases are robust to such spurious effects and in theory depart from zero only for dependent sources with actual delays, a crucial property for assessing functional connectivity. This observation has led to the design, study and use of the Imaginary Coherence function (F-ICOH) [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] and the (weighted-)Phase Lag Index (F-wPLI) [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF] as relevant indices to assess functional connectivity for the band-limited oscillatory brain activity measured by M/EEG measurements. Interested readers are referred to e.g., [START_REF] Varela | The brainweb: phase synchronization and large-scale integration[END_REF][START_REF] Engel | Dynamic predictions: oscillations and synchrony in top-down processing[END_REF][START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF][START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF][START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF][START_REF] Siegel | Spectral fingerprints of large-scale neuronal interactions[END_REF] for thorough reviews and further details (see also Section 2.1 for definitions). Beyond second order statistics and linear correlation, higher order (nonlinear) dependencies have also been investigated using directed partial correlations; moreover, the Granger causality approach has been used to infer causal links, see [START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF] for a review.

Functional connectivity was so far mainly measured via the band-limited oscillatory activity of the brain, and has hardly been applied to characterize the infraslow arrhythmic brain activity. Preliminary attempts in that direction [START_REF] Achard | Fractal connectivity of long-memory networks[END_REF][START_REF] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks[END_REF], though based on wavelet representation, remained tied to the coherence function, hence essentially to direct correlation, and thus severely impaired by volume conduction effect in functional connectivity assessment in M/EEG. This lack of functional connectivity tools dedicated to the infraslow regime is partly due to the role infraslow arrhythmic temporal dynamics to brain activity remaining controversial, but also, and prominently, because conceptual and practical tools reconciling the modeling and analysis of both multivariate and scale-free dynamics were lacking. This situation changed recently with the theoretical definition and formal study of multivariate selfsimilarity [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF] as well as with the design and assessment of multivariate wavelet transform based practical tools [START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF]Abry and Didier, 2018a,b;Abry et al., 2019a,b), thus permitting the investigation of functional connectivity within the infraslow arrhythmic brain activity, at the core of the present work.

Goals, contributions and outline. The present work aims to revisit the analysis of functional connectivity in human brain activity in two ways: First, functional connectivity assessment will be based on the on-going (or spontaneous) infraslow arrhythmic (or scale-free) activity of the human brain, rather than on stimulusinduced band-limited oscillatory faster rhythms. This will be referred to as functional connectivity assessed from fractal dynamics (see (La Rocca et al., 2018a) for a preliminary attempt).

Second, indices quantifying functional connectivity from fractal dynamics will be constructed from multivariate complex wavelet transforms, rather than from Fourier-based cross-spectral analysis. The key intuitions underlying the design of these indices are double: Based on wavelet transforms, these tools will inherit from their well-documented performance and robustness for the analysis of scale-free dynamics [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motions[END_REF][START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF][START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF][START_REF] Abry | Wavelets for the analysis, estimation, and synthesis of scaling data (Self-Similar Network Traffic and Performance Evaluation[END_REF][START_REF] Veitch | A statistical test for the time constancy of scaling exponents[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF]; Complex wavelets allow to incorporate phase information in the analysis of multivariate cross-temporal dynamics.

To that end, after a brief recall of Fourier-based spectral estimation and of the classical The proposed Complex Wavelet indices assessing functional connectivity from fractal dynamics are extensively tested on MEG data, collected on 36 individuals, both at rest and during a visual discrimination learning task. The experimental data is described in Section 3 (see also [START_REF] Zilber | Supramodal processing optimizes visual perceptual learning and plasticity[END_REF]).

Analyses of functional connectivity assessed from fractal dynamics within infraslow arrhythmic cross temporal dynamics regime, ranging from 0.1 to 1.5Hz for this data set [START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF], are reported in Section 4 and discussed in Section 5.

The proposed Complex Wavelet indices are demonstrated to have a high sensitivity in capturing significant and meaningful group-level functional connectivity assessed from fractal dynamics networks both at rest and during task performance, that present longrange spatial interactions between fronto-occipital and temporo-parietal brain regions.

Further, a significant increase in functional connectivity assessed from fractal dynamics is shown to be positively correlated with behavioral performance in the task and to be further reinforced by the training stage and thus by learning. Finally, our results suggest an interplay between temporal and spatial dynamics: Arrhythmic infraslow brain activity evolves from strongly and globally structured slow temporal dynamics for each region individually at rest, related across the brain by a clear functional network, to faster and less globally structured temporal dynamics per region, yet with significantly stronger spatial couplings across the brain, during task.

The proposed Complex Wavelet tools constitute, to the best of our knowledge, the first operational tools for a relevant assessment of functional connectivity from fractal dynamics, i.e., functional connectivity in scale-free cross-temporal dynamics for the largeband infraslow arrhythmic brain activity recorded in M/EEG. MATLAB codes, designed and implemented by ourselves, for the synthesis of multivariate scale-free synthetic data and for the computation of Complex Wavelet-based indices to assess functional connectivity from fractal dynamics, will be made publicly available at the time of publication.

METHODOLOGY: FUNCTIONAL CONNECTIVITY ASSESSMENT

Frequency domain functional connectivity assessment

The M -variate data (X m (t) m=1,...,M , t ∈ R) available for analysis are assumed to be real-valued finite power realizations of stochastic processes, with well-defined power cross-spectral density S m,m (f ). Welch periodogram constitutes one of the classical nonparametric spectral estimation procedures [START_REF] Papoulis | Signal analysis[END_REF], grounded on the use of a windowed Fourier transform. This Fourier-based estimate S (F ) of the cross-spectrum S is indeed defined as a time average of the squared-moduli of the windowed (or short-time) Fourier transform coefficients g X ( , k) = X(t)φ ,k (t)dt:

S (F ) m,m (f = ν 0 ) = k g Xm ( , k)g * X m ( , k), (1)
where φ ,k (t) = φ(t -kT 0 ) exp (-2ı ν 0 t) denotes the collection of translated and frequency-shifted templates of a reference pattern φ(t), and T 0 and ν 0 positive constants that can be arbitrarily chosen, provided that they satisfy T 0 ν 0 ≤ 1/(4π).

Straightforward calculations yield

ES (F ) m,m ( ν 0 ) = S m,m ( ν 0 -f )| φ(f )| 2 df, (2)
with φ denoting the Fourier transform of φ and E the ensemble average. This is thus showing that S (F )

m,m provides a biased estimate of S m,m (f ). The time and frequency resolutions of the functions φ ,k being uniformly controlled by the choice of the function φ, S (F ) achieves a fixed absolute-frequency resolution multivariate spectral analysis.

From S (F ) m,m (f ), three functions are classically involved in functional connectivity assessment, the modulus (F-COH), the Imaginary (F-ICOH) part of the coherence function [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF], and the weighted Phase Lag Index (F-wPLI) [START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF] (with the imaginary part of a complex number):

F-COH m,m (f ) S F m,m (f ) S F m,m (f )S F m ,m (f ) , (3)
F-ICOH m,m (f ) {S F m,m (f )} S F m,m (f )S F m ,m (f ) , (4)
F-wPLI m,m (f = ν 0 ) n j k=1 g Xm ( ,k)g * X m ( ,k) n j k=1 g Xm ( ,k)g * X m ( ,k) } . ( 5
)
To quantify functional connectivity on MEG signals, the corresponding indices are practically computed as sums of the absolute values of these functions over the range of frequencies defining the targeted band-limited oscillations. Large values (above predefined thresholds) are used as markers of functional connectivity at the individual level, usually followed by statistical testing for assessing group-level significance.

Wavelet domain functional connectivity assessment

Complex wavelet transform. The classical discrete wavelet transform relies on the use of a real-valued mother-wavelet (cf. e.g., [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]). To assess phases and delays amongst signals, it is proposed here to use a complex wavelet transform, defined as follows.

Let ψ (r) denote a real-valued oscillating and sufficiently smooth reference pattern, referred to as the mother wavelet, and constructed such that the collection of dilated and translated

templates {ψ j,k (t) = 2 -j/2 ψ(2 -j t -k)} (j,k)∈Z 2 of ψ form an orthonormal basis of L 2 (R)
(cf. e.g., [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]). From ψ (r) , an analytic complex mother-wavelet can be defined as ψ = ψ (r) + ıψ (ı) , where ψ (ı) consists of the Hilbert transform of ψ (r) . The design of a complex, invertible and analytic mother wavelet is not straightforward. In the present work, we build on the excellent approximate solution proposed in [START_REF] Kingsbury | Complex wavelets for shift invariant analysis and filtering of signals[END_REF][START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF], referred to as the dual-tree complex wavelet transform.

For a signal X, the coefficients of the dual-tree complex wavelet transform are

defined as d X (j, k) d (r) X (j, k) + ıd (ı) X (j, k), with d (r) X (j, k) ψ (r) j,k (t)X(t)dt and d (ı) X (j, k) ψ (ı) j,k X(t)dt.
Computing a dual-tree complex wavelet transform thus amounts to computing two standard Discrete Wavelet Transforms, with the two real mother-wavelets r) and ψ (ı) , respectively, independently.

ψ (
Wavelet cross spectrum and functional connectivity. It has been well-documented that the study of univariate scale-free temporal dynamics should be performed using a wavelet-based spectral estimation rather than a Fourier-based one (cf. e.g., [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motions[END_REF][START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF]Veitch andAbry, 1999, 2001)). This has recently been extended to multivariate scale-free temporal dynamics analysis and wavelet cross-spectrum estimation (cf. e.g., [START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF]La Rocca et al., 2018a;[START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF][START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]). Given a pair of signals X m , X m , the multivariate wavelet (cross-)spectrum can be defined as

S W m,m (j) 1 n j n j k=1 d Xm (j, k)d * X m (j, k) (6)
where n j ≈ N 2 j are the number of coefficients available at scale j, and * stands for complex conjugate.

It has been shown [START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF] that

ES (W ) m,m (j = ν ψ /2 j ) = S m,m (f )| ψ(f /2 j )| 2 df, (7)
with ψ denoting the Fourier transform of ψ. This indicates that S (W ) m,m (j) estimates S m,m (f = ν ψ /a j 0 ) around frequency f = ν ψ /a j 0 and achieves a fixed relative-frequency resolution multivariate spectral analysis [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF].

Eqs. ( 2) and ( 7) combined together show that Fourier-based S 2) and ( 7) and converting scales a = 2 j into frequencies as f = f 0 × f s /2 j where f s is the data sampling frequency and f 0 a constant that depends on the specific choice of the mother wavelet. Readers interested by further theoretical and practical discussions on comparing Fourier and wavelet-based spectral estimations, are referred to e.g., [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF]Veitch andAbry, 1999, 2001;[START_REF] Abry | Wavelets for the analysis, estimation, and synthesis of scaling data (Self-Similar Network Traffic and Performance Evaluation[END_REF][START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF].

Wavelet-based functional connectivity indices. From the wavelet-based estimate of the power spectrum, wavelet-based indices can be constructed to assess functional connectivity, as was the case with Fourier spectrum, mutatis mutandis:

W-COH m,m (j) S W m,m (j) S W m,m (j)S W m ,m (j) , (8)
W-ICOH m,m (j)

S W m,m (j) S W m,m (j)S W m ,m (j) , (9)
W-wPLI m,m (j)

n j k=1 d Xm (j,k)d * X m (j,k) n j k=1 d Xm (j,k)d * X m (j,k) . ( 10
)
Unlike the standard discrete wavelet transform coherence function used in, e.g., [START_REF] Whitcher | Wavelet analysis of covariance with application to atmospheric time series[END_REF][START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF], W-COH m,m (j) is complex-valued.

Functional connectivity assessed from fractal dynamics. Functional connectivity for scale-free infraslow temporal dynamics, consists in averaging the absolute values of these functions over the corresponding range of octaves j 1 ≤ j ≤ j 2 (equivalently over the range Middle: Wavelet spectra ((c), (f)), cross spectrum (d) and coherence function (e) as functions of the (log of the) scales (top row, red lines). Bottom: Comparison to Fourier spectra ((g), (j)), cross-spectrum (h) and coherence function (i) (solid black lines) after remapping scales into frequencies (bottom row). The scale-free (or arrhythmic) regime is marked by linear behaviors of the power spectra across coarse scales, 8 j 12 corresponding to low frequencies, 0.1 f 1.5 Hz, in these log-log plots.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)
of scales a = 2 j or frequencies f = f 0 /2 j ) where scale-free dynamics are observed:

1 j 2 -j 1 + 1 j 2 j=j 1 W-wPLI m,m (j) or 1 j 2 -j 1 + 1 j 2 j=j 1 W-ICOH m,m (j)
Remapping scales into frequencies, calculations inspired from those leading to Eqs. ( 2) and ( 7) permit to compare theoretically and practically W-COH, W-ICOH and W-wPLI to F-COH, F-ICOH and F-wPLI, as illustrated in Figs. 2, 3 and 5 on synthetic data. This is here critical to emphasize that functional connectivity assessed from fractal dynamics as defined and used in the present work is associated with (the statistics of) cross-temporal dynamics. It should not be confused with the so-called fractal networks, as studied in e.g., in [START_REF] Bassett | Adaptive reconfiguration of fractal small-world human brain functional networks[END_REF][START_REF] Varley | Fractal dimension of cortical functional connectivity networks & severity of disorders of consciousness[END_REF], which are related to topological (thus static) properties of a spatial graph.

Functional connectivity from fractal dynamics performance assessment

Monte Carlo numerical simulations. To assess the performance of the proposed indices aiming to quantify functional connectivity from fractal dynamics, Monte Carlo numerical simulations were conducted. They make use of synthetic bivariate fractional Brownian motion, a specific instance of the multivariate selfsimilarmodel recently introduced in [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF] and studied in (Abry and Didier, 2018a,b). Bivariate fractional Brownian motion consists of a pair of fractional Brownian motions B H 1 and B H 2 , with possibly different selfsimilarity parameters H 1 and H 2 , with pointwise correlation ρ. In addition, one component is delayed by ∆. Correlation coefficient ρ is set to range within ρ ∈ {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9} and delays range in ∆ = {0, 1, 2, 4, 8, 16, 32, 64} samples. Sample size is n = 2 14 , chosen to match the size of the infraslow regime of the MEG data (cf. Sections 3 and 4).

To model MEG data as those analyzed in Section 4 and as commonly indicated in the literature [START_REF] He | The temporal structures and functional significance of scale-free brain activity[END_REF], one needs to use both fractional Gaussian noise (fGn), the increments of fractional Brownian motion (fBm), with parameter H ranging from say 0.6 to 1 and fractional Brownian motion itself with parameters ranging from 0 to 1. Therefore, the numerical simulations conducted here were based on bivariate processes, each component being either fGn or fBm, with 0 < H < 1. For the Fourier-based spectral estimation, the classical averaged windowed periodogram estimate of the power spectral density was computed, with Hanning windows of width corresponding with the frequency bands of the complex wavelet filters, to enable relevant comparisons of the tools. For the Complex-Wavelet based estimation, q-shift complex wavelets were used, as described in [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF] and references therein, (see, e.g., [START_REF] Lina | Complex daubechies wavelets[END_REF] for an alternative choice).

Indices assessing functional connectivity from fractal dynamics (both Fourier and waveletbased) were computed as average over a range of frequencies and scales that match those of the infraslow scale-free range observed on the MEG data described and analyzed hereafter.

Performance are reported as means (and confidence intervals) computed from N = 1000 independent realizations of bivariate fractional Gaussian noise.

Spurious connectivity. To start with, we analyzed scenarii where the two components of bivariate fractional Gaussian noise were correlated but not delayed: ∆ ≡ 0. ii) As functions of j, W-ICOH and W-wPLI display different patterns that depend on ∆.

However, these patterns both show independently maximum absolute deviations from 0 at scales that vary with ∆ (second column). This was quantified for W-ICOH and used as a delay estimation procedure [START_REF] Didier | Detection and estimation of delays in bivariate self-similarity: bootstrapped complex wavelet coherence[END_REF].

iii) When a scale 2 j in relation to ∆ is chosen, both (the absolute values of) W-ICOH and W-wPLI are proportional to (the absolute value of) ρ (third column). This shows not only that W-ICOH and W-wPLI depart from 0 when delays amongst components exist, but also that the amplitude of the departure is proportional to the correlation ρ between components, a crucial property to assess quantitatively functional connectivity, clearly and originally quantified in these numerical simulations. ii) F-ICOH and F-wPLI significantly overestimate correlations at small ρ while W-COH and W-wPLI do not (third column) ;

iii) The RMSE of F-ICOH and F-wPLI becomes up to ten times larger than RMSE of W-ICOH and W-wPLI for small values of ρ (fourth column).

Functional connectivity from fractal dynamics assessment performance. In addition, Fig. 6 compares the ratio of the RMSE of W-ICOH to the RMSE of W-wPLI i) They indicate that W-COH cannot be used to assess functional connectivity as it is fooled by zero-delay (volume conduction effect) correlations, thus confirming an already documented observation for F-COH in the literature [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF][START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF].

To the converse, W-ICOH and W-wPLI (and F-ICOH and F-wPLI) are much less affected by these spurious correlations.

ii) The Complex Wavelet W-ICOH and W-wPLI can be used to assess functional connectivity for scale-free temporal dynamics.

iii) The Complex Wavelet W-ICOH and W-wPLI perform significantly better than the Fourier-based F-COH and F-wPLI, first, when the signals show very large scaling exponents β in their f -β power spectral density behavior, as is the case with fBm-like time series, or, second, when additive noise in the form of smooth and slow trends are superimposed to data with scale-free dynamics, a situation commonly observed in recordings collected from neuroimaging techniques. iv) W-ICOH and W-wPLI perform comparably with (slightly) better performance of W-wPLI when ρ or ∆ increases, or when smooth trends are superimposed to scale-free dynamics, as often the case on MEG data. This will be further discussed in Section 4.

The proposed complex wavelet based assessment of functional connectivity in infraslow arrhythmic brain activity was tested on MEG measurements, consisting of non-invasive recordings of simultaneous time-series reflecting the whole brain activity, both at rest and during the completion of a task. All details about the experimental paradigm and the task can be found in [START_REF] Zilber | Supramodal processing optimizes visual perceptual learning and plasticity[END_REF].

In short, the task was designed from a short-term learning paradigm and consisted of a visual coherence discrimination. Two sets of colored (green and red) dots were mixed and shown on a screen, each dot with random and independent movement. After a variable duration interval (0.3 to 0.6 s) of incoherent motion, a fraction of randomly chosen dots belonging to either of the two sets (also randomly chosen at each trial) followed a coherent motion during one second. Participants were asked to tell which of the red or green clouds had a coherent motion by pressing a button of the same color. Task difficulty was increased by decreasing the rate of dots in coherent motion.

The experiment was organized as interleaved MEG blocks alternating rest and task measurements: It started with a 5-minute rest recording (REST Brain activity was recorded via MEG modality, in a magnetically shielded room using a 306 MEG system (Neuromag Elekta LTD, Helsinki). MEG signals originally sampled at 2 kHz were downsampled at 448 Hz, and preprocessed to remove external and internal interferences, in accordance with accepted guidelines for MEG research [START_REF] Gross | Good practice for conducting and reporting MEG research[END_REF]. Signal Space Separation (SSS) was applied with MaxFilter to remove exogenous artifacts and noisy sensors [START_REF] Taulu | Spatiotemporal signal space separation method for rejecting nearby interference in MEG measurements[END_REF]. Ocular and cardiac artifacts (eye blinks and heart beats) were removed using Independent Component Analysis (ICA) on raw signals. ICA were fitted to raw MEG signals, and sources matching the ECG and EOG were automatically found and removed before signals reconstruction, following the procedure described in [START_REF] Jas | Autoreject: Automated artifact rejection for MEG and EEG data[END_REF]. Source localization from MEG signals was used to reconstruct source cortical activity using the mne analyze tools within MNE [START_REF] Gramfort | Meg and eeg data analysis with mne-python[END_REF]. Details regarding the source localization technique are reported in [START_REF] Zilber | Supramodal processing optimizes visual perceptual learning and plasticity[END_REF]. Finally, following analyses reported in [START_REF] Zilber | Supramodal processing optimizes visual perceptual learning and plasticity[END_REF][START_REF] Rocca | Emergence of β and γ networks following multisensory training[END_REF], 28 cortical regions-of-interest (ROIs), recruited in task performance (including frontal, somato-sensory, temporal, parietal and occipital areas) were retained for the analysis of functional connectivity in infraslow temporal dynamics.

FUNCTIONAL CONNECTIVITY ASSESSED FROM FRACTAL DYNAMICS IN INFRASLOW ARRHYTHMIC MEG-RECORDED BRAIN ACTIVITY

Infraslow scale/frequency range and functional connectivity from fractal dynamics assessment methodology

Infraslow scale/frequency range. Following the systematic inspections of the wavelet spectra and cross-spectra reported in (La [START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF] for the same MEG data, the scale-free range of scales is set uniformly for the 28 times series and across the 36 participants, for the analysis of arrhythmic functional connectivity to 8 j 12 , thus corresponding to frequencies in 0.1 f 1.5Hz or equivalently to time scales ranging roughly from 1 to 10s. This scale-free regime is illustrated in Fig. 1 Tests. These indices were filtered at the group-level (N = 36), using a recently introduced network density threshold method, the Efficiency Cost Optimization [START_REF] Vico Fallani | A topological criterion for filtering information in complex brain networks[END_REF], thus yielding group-level 28 × 28 fractal dynamics-based functional connectivity matrices across the brain for each experimental condition independently. See also [START_REF] Rocca | Emergence of β and γ networks following multisensory training[END_REF] for further details on the use of such technique.

To investigate significant differences in infraslow functional connectivity between two different experimental conditions (e.g., TASK i -REST i ), independently for each chosen index, a group-level paired t-test was performed, with a demanding preset significance level: p < 0.01. The false discovery rate (FDR) procedure was used to correct p-values for multiple comparisons across the 28 × 27/2 possible connections.

Comparisons against Fourier-based indices. To compare Fourier-based F-ICOH and F-wPLI to Complex Wavelet-based W-ICOH to W-wPLI, Fourier-based spectral estimation was conducted using Welch Periodogram procedures (as described in Section 2.1), using a windowed Fourier transform with a Hanning-type window of duration 80s.

Figure 8. Fractal dynamics-based functional connectivity assessment (W-wPLI) differences between REST i and TASK i and between REST i and TASK f . The increase in functional connectivity assessed from fractal dynamics from rest to task is strengthened with training, i.e. from TASK i to TASK f , and emerged between several intra-or interhemispheric pairs of regions (Frontal polar/IPS, ITC/MT, FEF/pSTS) involved in task performance. i) The connectivity networks yielded by W-COH predominantly display short-range and inter-hemispheric interactions throughout the cortex and most notably amongst frontal regions on one hand and temporo-occipital regions on other hand, both for REST i and TASK i .

Fractal dynamics-based functional connectivity networks

ii) The connectivity networks yielded by W-ICOH and W-wPLI display similar structures, dominated by long-range spatial interactions, that differ significantly from those of the networks produced by W-COH, dominated by shorter-range spatial interactions. These differences in network structures can be quantified using the Average Degree, i.e., the average number of connections per node, as a network structure metrics. For REST i , the Average Degrees for the graphs obtained by W-COH, W-ICOH and W-wPLI are respectively of 0.95(±0.37), 0.21(±0.24) and 0.44(±0.52). Medians in distributions of number of links per node differ significantly between W-COH and W-ICOH (p < 10 -11 ) or between W-COH and W-wPLI (p < 10 -6 ). The same holds for TASK i , with average degrees of respectively 1.0(±0.49), 0.25(±0.24) and 0.52(±0.50), and significance of respectively p < 10 -8 and p < 10 -3 .

iii) While yielding comparable networks, W-wPLI and W-ICOH differ insofar as the former yields larger connectivity indices than the latter. In addition, connectivity networks using W-wPLI or W-ICOH differ in structure, however much less than when comparing W-wPLI vs. W-COH or W-ICOH vs. W-COH. Indeed, for REST i the Average Degrees of W-wPLI and W-ICOH are respectively of 0.44(±0.52) and 0.21(±0.24), yielding a quantifiable difference (p = 0.04), and for TASK i the Average Degrees of W-wPLI and W-ICOH are respectively of 0.52(±0.50) and 0.25(±0.24), yielding a clearer difference (p = 0.01). iv) When comparing TASK i versus REST i , W-wPLI and W-ICOH both indicate an increase in functional connectivity during task performance. This increase in functional connectivity assessed from fractal dynamics highlights interactions between regions recruited in the achievement of the task, notably fronto-temporal couplings (between the right ventro-lateral prefrontal cortex (vlPFC) and inferior temporal cortex (ITC)), interactions linking temporal regions (anterior superior temporal sulcus (aSTS) and auditory cortex) with the intra-parietal sulcus (IPS), motor-occipital couplings between the left frontal BA6 (including premotor and supplementary motor regions) and primary visual areas (V1/V2). Interaction between the key region hMT+, sensitive to visual motion, and the associative area, pSTS, is also significant in the left hemisphere.

Focusing on the W-wPLI index only, Fig. 8 shows the additional comparisons of the post-training task TASK f to the initial rest REST i , which, compared to the contrast TASK i -REST i (cf. Fig. 7 bottom left plot), indicates first that functional interactions in infraslow temporal dynamics are globally strengthened by the training and second that new intra-and inter-hemispheric couplings emerged with training involving much more the parieto-occipito-temporal network (IPS, primary visual cortex and anterior STS). We also noticed new interactions between the left fronto-polar region and the left IPS, between the right frontal eye fields (FEF) and the pSTS and finally between the BA6 complex and hMT+ region.

Figure 9. selfsimilarity (H) differences between REST i and TASK i and between REST i and TASK f . The decrease in selfsimilarity from rest to task is strengthened with training, i.e. from TASK i to TASK f , and more heavily in the parieto-occipital (hMT+, visual cortices, V1/V2/V4) regions involved in task performance. Note that a value of H was computed per cortical label here. See [START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF] for methodological details.

Functional connectivity assessed from fractal dynamics and selfsimilarity

In (La [START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF], selfsimilarity was systematically quantified by waveletbased measurements of the selfsimilarity exponent H and a global decrease from rest to task was observed over the whole brain (see Fig. 4E in (La Rocca et al., 2018b)). This result, obtained from 24 participants, is here strengthened by using 36 subjects. Fig. 9 reports a decrease in H not only between REST i and TASK i but also between REST i and TASK f . Further, Fig. 9 shows a strengthening of the decrease in H from TASK i to TASK f in the parieto-occipital regions involved in task performance, notably the bilateral hMT+ regions, the visual cortices including V1/V1 and V4 for the visual color detection.

Interestingly, after training these regions are also more strongly coupled with others during task performance (TASK f vs REST i ).

To investigate a potential training-induced relation between the decrease in selfsimilarity and the increase in W-wPLI, ∆H = H TASF f -H REST i and ∆ W-wPLI = W-wPLI TASF f -W-wPLI REST i were averaged across the whole brain for each subject. Corresponding averages are shown in Fig. 10 which interestingly suggests a significant (p = 0.05) anticorrelation

Complex-Wavelet assessment

Averaging (the absolute values) of F-wPLI across a range of frequencies that match the range of scales associated with the infraslow scale-free scaling range, permits to compare Fourier-assessed functional connectivity from fractal dynamics. Fig. 12 reports the density networks obtained from F-wPLI for REST i and TASK i , showing significant differences with those obtained using W-wPLI. The network topography associated with the F-wPLI index are denser compared to W-wPLI. Indeed, using the Average Degree, used as a graph structure metrics, it was found that for REST i , the Average Degrees of W-wPLI and F-wPLI are respectively of 0.44(±0.52) and 1.62(±1.11), yielding a very significant difference, assessed by a p-value below 6 × 10 -6 , and for TASK i , the Average Degrees of W-wPLI and F-wPLI are respectively of 0.52(±0.50) and 1.65(±1.21), yielding also a significant difference, assessed by a p-value of 5 × 10 -5 . Further, the number of significant interactions with F-wPLI is more balanced between the two hemispheres during REST i in contrast to W-wPLI, which captures more couplings in the right one. Also, the resting-state W-wPLI-based network configuration is more dominated by fronto-occipital couplings whereas the F-wPLI-based shows a larger number of inter-hemispheric interactions. During the pre-training task TASK i , W-wPLI and F-wPLI network topographies both show similar connections but also strong differences: the former is more dominated by fronto-parietooccipital couplings with a hub role played by the visual cortices, while the latter does not strongly differ from the F-wPLI network found during REST i . Finally and more importantly, no statistically significant difference in F-wPLI TASK i -F-wPLI REST i can be evidenced (see Fig. 12-bottom), while a significant increase in W-wPLI was found from REST i to TASK i between fronto-parieto-occipital regions that are involved in task performance (see Fig. 9- top). The coupling between V4 and MT in the right hemisphere reflects the color-motion binding, while the significant interactions involving the anterior STS, IPS and vlPFC are likely due to their role in multisensory processing. Hence, the W-wPLI index provides much more meaningful information when contrasting rest to task brain activity.

DISCUSSION

Functional connectivity from fractal dynamics assessment. At the methodological level, the results presented in Section 4 clearly showed that W-COH fails to characterized correctly functional connectivity, in clear agreement with the numerical simulations reported in Section 2.3 on synthetic data fGn/fBm and with results reported in the literature (cf. [START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF][START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF]).

More interestingly, compared to W-ICOH, W-wPLI was observed to more accurately quantify functional connectivity assessment from fractal dynamics, both at rest and during task in MEG data, as well as to better highlight relevant changes in functional connectivity assessed from fractal dynamics between rest and task. This is in agreement with previously reported results, showing that for band-limited oscillatory activities, F-wPLI was a better index to assess functional connectivity than F-ICOH. This was attributed to the denominator of F-wPLI being different from that of F-ICOH and less sensitive to (residual) volume conduction effects [START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF][START_REF] Vinck | An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias[END_REF]. These arguments straightforwardly extend to W-wPLI and W-ICOH thus likely explaining the enhanced ability of W-wPLI to assess functional connectivity from fractal dynamics compared to W-ICOH. Interestingly, the numerical simulations conducted in Section 2.3 on synthetic fGn/fBm data showed only a moderate superiority of W-wPLI over W-ICOH to quantify functional connectivity from fractal dynamics, except for slightly improved estimation (RMSE) performance. This suggests that fGn/fBm, even with delays, correlations and possible additive trends, are not rich enough models to account for all the difficulties encountered in modeling real MEG data. This is calling for richer modeling, potentially involving multifractality. This will be further explored.

The benefits of using wavelet-based (multiscale) tools to analyze scale-free temporal dynamics and estimate the corresponding scaling exponent compared to classical Fourierbased spectral estimation have been abundantly documented elsewhere (cf. e.g., [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF]Veitch andAbry, 1999, 2001;[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] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF]). First, they provide better (unbiased and controlled variance) estimates of H;

Second, by tuning the so-called number of vanishing moments of the mother wavelet [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], wavelet-based spectral estimation is robust to additive smooth slow trends in data which are, to the converse, strongly altering Fourier-based spectral estimation. These benefits are straightforwardly inherited by the wavelet-based indices for assessing functional connectivity from fractal dynamics. This was clearly evidenced by the numerical simulations reported in Section 2.3 showing the robustness to trends and the better performance for large scaling exponents of Complex Wavelet-based indices over Fourier-based ones.

Functional connectivity assessed from fractal dynamics in time relates to long-range spatial interactions. On MEG data, functional connectivity in the infraslow arrhythmic regime assessed by W-COH, i.e., based on direct correlation, was observed to yield mostly spatial short-range connectivity networks across the brain, notably with spurious shortrange functional intra-and inter-hemispheric interactions, visible between frontal regions both at rest and during task. This is likely a consequence of residual common source effects, strongly biasing the real part of thecoherence function, and thus yielding spurious connectivity measures, in agreement with results reported in [START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF]. In contrast, functional connectivity assessed by W-ICOH and W-wPLI indices, i.e., based on phase coupling, did not show such short-range links, but rather functional connectivity patterns dominated by long-range spatial interactions. This yields the first major result of the present work: Functional connectivity pertaining to the large-band infraslow arrhythmic temporal dynamics, (from 1 to 10s, or equivalently from 0.1 to 1Hz), reveals long-range spatial interactions, evidencing notably couplings between frontal, parietal and occipital brain regions. Functional connectivity assessed from fractal dynamics thus permits to quantify phase couplings and interactions associated with large lags. This departs from functional connectivity networks produced by the analysis of band-limited oscillatory temporal dynamics, that pertains to the fast (high frequency) brain activity and thus focuses on short time delays.

Functional connectivity assessed from fractal dynamics increases during task performance and with training. Compared to F-wPLI, W-wPLI showed an enhanced statistical sensitivity as it revealed a positively engaged parieto-temporo-occipital network in infraslow temporal dynamics when contrasting rest to pre-training activities. This network comprises previously identified key brain regions (e.g. hMT+, ITC, vlPFC, pSTS) during task performance. Interestingly, such regions also consistently identified as recruited by task when using standard temporal or spectral data analysis [START_REF] Zilber | Supramodal processing optimizes visual perceptual learning and plasticity[END_REF][START_REF] Rocca | Emergence of β and γ networks following multisensory training[END_REF]. However, W-wPLI was the only index further showing that functional connectivity assessed from fractal dynamics actually increased during task performance in these regions.

A second key result consists of the observation of the strengthening of this functional connectivity from fractal dynamics based functional network with training, i.e. when contrasting rest to post-training activity. It shows the rising of new key couplings between frontal and parieto-temporal cortices that suggest that some cortical representations of the visual detection and decision making process may emerge even at slow time scales (1 s to 10 s) and may be used as a substrate for facilitating faster dynamics in oscillatory regimes.

Such increased functional connectivity assessed from fractal dynamics is a hallmark of brain plasticity induced by the training stage.

The third finding of this study is the positive correlation between the increase in functional connectivity assessed from fractal dynamics and task performance when contrasting preto post-training brain activity. This suggests that the consolidated network eases task completion for each individual, experiencing averaged increase in functional couplings within the infraslow regime.

Functional connectivity from fractal dynamics and selfsimilarity quantifying an interplay between temporal and spatial dynamics. Finally, the increase in functional connectivity assessed from fractal dynamics was shown to be correlated with a decrease in selfsimilarity from rest to task. These results on functional connectivity assessment from fractal dynamics, combined with the univariate (regionwise) analysis of scale-free temporal dynamics of the same data (La [START_REF] Rocca | Selfsimilarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF], lead to the following global picture for the large-band arrhythmic infraslow temporal dynamics of brain activity.

At rest, each region displays a globally very structured and slow activity in time

CONCLUSIONS

In this work, we have introduced the notion of functional connectivity assessment from fractal dynamics for MEG data, defined as functional connectivity associated with the large-band infraslow (typically below the Hz) arrhythmics (scale-free) cross temporal dynamics, in contradistinction with the classical functional connectivities associated with the band-limited rapid oscillatory rhythms (α-, β-, γbands).

It has been argued and demonstrated that complex wavelet (multiscale) based analyses permit to construct indices to assess functional connectivity from fractal dynamics that inherit from the theoretical and practical benefits of wavelet representations for scale-free (cross-temporal) dynamics analysis, notably in terms of robustness to trends and large selfsimilarity parameters H. It was confirmed that wPLI outperforms ICOH, as commonly observed and that COH is not suited for functional connectivity assessment.

While Fourier-based tools are natural to use to assess functional connectivity in band-limited rapid oscillatory rhythms, it was shown, using simulated synthetic data and mostly on MEG data, that the assessment of functional connectivity for large-band slow scale-free cross-temporal dynamics is better achieved by complex wavelet based indices. Therefore, Fourier and complex wavelet based spectral estimation must be regarded as complementary, rather than as mutually exclusive, tools. Routines (in MATLAB) to synthesize (correlated and delayed) bivariate fractional Gaussian noise, to perform Fourier and complex-wavelet based analysis and to compute indices quantifying functional connectivity from fractal dynamics, on synthetic or MEG data, have been developed by ourselves and will be made publicly available at the time of publication.

Such tools could further be used to examine the relevance of functional connectivity assessed from fractal dynamics in the context of network physiology, and of networks of networks, relating brain activity to other physiological functions (heart rate, respiration, sleep, ocular and motor systems,. . . ), cf. e.g., [START_REF] Bartsch | Network physiology: how organ systems dynamically interact[END_REF][START_REF] Bartsch | Network physiology: how organ systems dynamically interact[END_REF][START_REF] Liu | Plasticity of brain wave network interactions and evolution across physiologic states[END_REF][START_REF] Catrambone | Quantifying functional links between brain and heartbeat dynamics in the multifractal domain: A preliminary analysis[END_REF].

  Fourier-based functional connectivity indices (F-ICOH and F-wPLI) in Section 2.1, Complex wavelet transforms and the corresponding Complex Wavelet-based functional connectivity indices (W-ICOH and W-wPLI) are defined in Section 2.2. The performance of several Complex Wavelet-based functional connectivity indices proposed here are compared ones against the others, and against their corresponding Fourier counterparts, by means of Monte Carlo numerical simulations, involving a large number of independent drawings of synthetic signals, sampled from stochastic processes commonly used to model scale-free temporal dynamics, multivariate fractional Brownian motions and multivariate fractional Gaussian noises (cf. Section 2.3). Several scenarii (different temporal dynamics, connectivity networks, additive trends) are investigated to assess the interest and relevance of the proposed Complex Wavelet indices (W-ICOH and W-wPLI) compared to Fourierbased ones, in terms of estimation performance and robustness to trends.

  constitute two biased estimates of the power spectral density S m,m , that can be compared theoretically and practically, as illustrated in Fig.1. Interested readers are referred to[START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF][START_REF] Abry | Multivariate scale-free temporal dynamics: From spectral (fourier) to fractal (wavelet) analysis[END_REF] for further discussions. As an illustration, the wavelet spectra and cross-spectrum of the two MEG signals displayed in Fig. 1 (a)-(b) are shown in Fig. 1(c)-(f) and compared to Fourier spectra and cross spectrum (cf. Fig. 1(g)-(j)), using Eqs. (

Figure 1 .

 1 Figure 1. Fourier versus wavelet spectral estimation on actual source-reconstructed MEG time series. Top: Two source-reconstructed MEG time series X 1 (a) and X 2 (b).Middle: Wavelet spectra ((c), (f)), cross spectrum (d) and coherence function (e) as functions of the (log of the) scales (top row, red lines). Bottom: Comparison to Fourier spectra ((g), (j)), cross-spectrum (h) and coherence function (i) (solid black lines) after remapping scales into frequencies (bottom row). The scale-free (or arrhythmic) regime is marked by linear behaviors of the power spectra across coarse scales, 8 j 12 corresponding to low frequencies, 0.1 f 1.5 Hz, in these log-log plots.

Figure 2 .

 2 Figure 2. Complex Wavelet-based connectivity on synthetic bivariate fractional Gaussian noise with correlation but no delay. W-COH (a), W-ICOH (b) and W-wPLI (c) as function of octaves j and correlation coefficient ρ.As it should, W-COH correctly assesses correlations with no delays and thus departs from 0 at all scales. W-COH would hence lead to incorrectly assessing functional connectivity. In contrast, W-ICOH and W-wPLI show averages values of 0 at all scales, and across all correlation levels, thus leading to assess no connectivity, as expected for non delayed components.

Fig 3 and

 3 Fig 3 and Fig 4 show that: i) Both W-ICOH and W-wPLI do depart from 0 across j and ρ when ∆ = 0 (left column).

Figure 3 .

 3 Figure 3. Complex Wavelet-based connectivity on synthetic bivariate fractional Gaussian noise with correlation and delay ∆ = 8. Top row: W-ICOH results. Bottom row: W-wPLI results. From left to right: W-ICOH (a) and W-wPLI (e) as functions of octaves j and correlation coefficient ρ; W-ICOH (b) and W-wPLI (f) as functions of octaves j, for a given ρ; W-ICOH (c) and W-wPLI (g) as functions of ρ for given octaves j; Ratio of the RMSE of F-ICOH to W-ICOH (d) and ratio of RMSE of F-wPLI to W-wPLI (h), averaged across scales 3 ≤ j ≤ 7, and color-coded in red as functions of delay ∆ and correlation coefficient ρ. A ratio larger than the value of 1 (made explicit to ease comparisons by horizontal blue plans) indicates poorer performance for Fourier-based estimates compared to wavelet-based ones. Synthetic data consists of bivariate fGn with H 1 = 0.7 and H 2 = 0.8.

Figure 4 .Figure 5 .

 45 Figure 4. Complex Wavelet-based connectivity on synthetic bivariate fractional Brownian motion with correlation and delay ∆ = 8. Top row: W-ICOH results. Bottom row: W-wPLI results. From left to right: W-ICOH (a) and W-wPLI (e) as functions of octaves j and correlation coefficient ρ; W-ICOH (b) and W-wPLI (f) as functions of octaves j, for a given ρ; W-ICOH (c) and W-wPLI (g) as functions of ρ for given octaves j; Ratio of the RMSE of F-ICOH to the RMSE of W-ICOH (d) and ratio of the RMSE of F-wPLI to the RMSE of W-wPLI (h), averaged across scales 3 ≤ j ≤ 7, and color-coded in red as functions of delay ∆ and correlation coefficient ρ. A ratio larger than the value of 1 (made explicit to ease comparisons by horizontal blue plans) indicates poorer performance for Fourier-based estimates compared to wavelet-based ones. Synthetic data consists of bivariate fBm with H 1 = 0.7 and H 2 = 0.8.

Figure 6 .

 6 Figure 6. Ratio of the RMSE of W-ICOH to the RMSE of W-wPLI, averaged across scales 3 ≤ j ≤ 7, as functions of delay ∆ and correlation coefficient ρ, for the synthetic data in Figs. 3, 4 and 5. Horizontal blue plans indicate the constant level 1 to ease reading.

  i ), followed by a 12-minute pre-training block (TASK i ); this was followed by 4 successive 5-minute long individualized training blocks. Another 5-min resting-state block (REST f ) was recorded prior to a final 12-minute post-training block (REST f ), consisting of the same visual coherence discrimination task as in TASK i . During TASK i and TASK f , the motion coherence discrimination accuracy of each participant was assessed. Pre-training and post-training behavioral thresholds were computed for each participant as the visual coherence level associated with 75 % of correct responses (hit rate). During REST blocks, participants were instructed to keep eyes opened, and were not following any other explicit instruction, thus permitting the analysis of spontaneous fluctuations of brain activity from MEG recordings. For the experiment, 36 healthy participants (mean age: 22.1 +/-2.2) were recruited. All participants were right-handed, had normal hearing and normal or corrected-to-normal vision. Before the experiment, all participants provided a written informed consent in accordance with the Declaration of Helsinki (2008) and the local Ethics Committee on Human Research at NeuroSpin (Gif-sur-Yvette, France).

  for arbitrarily chosen MEG signals shown in Fig. 1(a)-(b).Experimental conditions. Infraslow functional connectivity was assessed for several experimental conditions: resting-state (REST i ), pre-training (TASK i ) and post-training (TASK f ) tasks, thus enabling us to assess changes in functional interactions from rest to task and modulations related to learning.Functional connectivity from fractal dynamics indices. Three proposed complex wavelet based indices were then computed to assess infraslow functional connectivity by averaging across octaves corresponding to the scale-free regime, 8 ≤ j ≤ 12, the functions W-COH(j), W-ICOH(j) and W-wPLI(j), resulting in 3 sets of 28 × 28 × 36 indices.

Figure 7 .

 7 Figure 7. Functional connectivity assessment from fractal dynamics: Group-level functional connectivity in infraslow MEG-source reconstructed brain dynamics. Filtered 28 × 28 connectivity networks measured from Complex Wavelet based W-wPLI (left), W-ICOH (middle) and W-COH (right), for REST i (top row) and pre-training TASK i (center row). The red color intensity codes for the values of the connectivity indices (ranging from 0 to 1 by construction). Functional connectivity differences between conditions TASK i and REST i when assessed significant by a group level FDR corrected t-test are displayed in bottom row. Color codes for the TASK i -REST i differences in the values of indices from blue (negative) to red (positive), thus indicating that only increases in functional connectivity are observed from REST i to TASK i .

Fig. 7

 7 Fig. 7 reports the 28 × 28 thresholded connectivity networks yielded by the Complex Wavelet based indices defined in Section 2, W-wPLI (left), W-ICOH (middle) and W-COH (right), for two different experimental conditions REST i (top row) and pre-training TASK i (center row). Further, Fig. 7(bottom row) reports the FDR corrected statistically significant differences between indices measured during TASK i and REST i . Fig. 7 leads to the following observations:

Figure 10 .

 10 Figure 10. Decrease of selfsimilarity vs. increase in functional connectivity assessed from fractal dynamics from rest to task. ∆H = H TASK f -H REST i as a function of ∆ W-wPLI = W-wPLI TASK f -W-wPLI REST i , averaged across the whole brain for each of the 36 participants (each marked as a dot), shows that the decrease of selfsimilarity correlates negatively (r = -0.33, p = 0.05) with the increase of functional connectivity assessed from fractal dynamics.

Figure 11 .

 11 Figure 11. Functional connectivity assessment from fractal dynamics vs. Task Performance. Individual performance in post-training task shows significant (p = 0.01) positive correlation (r = 0.45) with the difference in functional connectivity assessed from fractal dynamics from pre-to post-training, i.e., W-wPLI TASK f -W-wPLI TASK i . Each particpant is represented as a dot, outliers are marked with a ×.

Figure 12 .

 12 Figure 12. Fourier-based wPLI estimator in the scale-free regime. No significant difference between F-wPLI TASK i and F-wPLI REST i in arrhythmic regime can be found.

(

  large selfsimilarity exponent H, hence strong temporal autocorrelation) with no transient structures (no burstiness, no multifractality, (La Rocca et al., 2018b)). The regions are connected across the brain by a clear spatial structure, that of functional connectivity assessed from fractal dynamics, constructed on measures of infraslow arrhythmic interactions. During task performance, temporal dynamics in each region independently become less globally structured and faster (decrease in H hence globally less correlated) with transient dynamical structures for regions involved in the task (burstiness and multifractality (La Rocca et al., 2018b)). These changes in regionwise temporal dynamics are accompanied by stronger functional connectivity assessed from fractal dynamics, i.e., by stronger spatial structures connecting regions. This permits to conjecture an interplay between temporal and spatial dynamics for the large-band infraslow arrhythmic brain activity: A decrease in global temporal structures induces faster and transient temporal dynamics and is associated with an increase in spatial structures and interactions between remote brain regions. Interestingly, these modulations are further strengthened with training, i.e. when contrasting the post-training to the restingstate activity in comparison with the pre-training vs. rest contrast. Overall, such modulations of brain spatio-temporal dynamics can be conjectured as a hallmark of brain plasticity.

  Complex wavelet based analyses of functional connectivity assessment from fractal dynamics conducted on MEG data recorded on 36 participants at rest and during a visual discrimination task with individualized training, yielded several key conclusions. First, large-band infraslow arrhythmic cross temporal dynamics can be associated with long-range (fronto-temporo-occipital) spatial interactions. Second, functional connectivity from fractal dynamics increases during task performance (in a set of brain regions consistent with those evidenced by other analyses performed on the same data with different tools) and is strengthened with training. Interestingly, a larger overall fractal dynamics-based functional connectivity increase correlates with a better task performance (larger hit rate). Third, the increase in spatial structure (quantified by the increase in functional connectivity assessed from fractal dynamics) is accompanied by changes in temporal structures, combining a decrease in the global temporal correlations (quantified by a decrease in the selfsimilarity index) and the increased occurrence of local transient structures (quantified by an increase in multifractality). These spatiotemporal modulations are reinforced with intensive and individualized training to the task.

This is a provisional file, not the final typeset article

Frontiers

CONFLICT OF INTEREST STATEMENT

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

AUTHOR CONTRIBUTIONS

Frontiers Functional Connectivity assessed from fractal dynamics

The original experimental design and access to MEG data was provided by Virginie van Wassenhove.

The methodological question studied here was framed and conceptualized by Patrice Abry and Philippe Ciuciu.

The data analysis tool design, implementation and performance assessment and interpretation were performed by Herwig Wendt and Patrice Abry.

MEG data analysis, results production and interpretation were performed by Daria La Rocca and Philippe Ciuciu.

The article was written by Patrice Abry and Philippe Ciuciu.

FUNDING

Work supported by ANR-16-CE33-0020 MultiFracs, France.