Login using

You can login by using one of your existing accounts.

We will be provided with an authorization token (please note: passwords are not shared with us) and will sync your accounts for you. This means that you will not need to remember your user name and password in the future and you will be able to login with the account you choose to sync, with the click of a button.

Scaling temporal dynamics in functional MRI (fMRI) signals have been evidenced for a decade as intrinsic characteristics of ongoing brain activity (Zarahn et al., 1997). Recently, scaling properties were shown to fluctuate across brain networks and to be modulated between rest and task (He, 2011): notably, Hurst exponent, quantifying long memory, decreases under task in activating and deactivating brain regions. In most cases, such results were obtained: First, from univariate (voxelwise or regionwise) analysis, hence focusing on specific cognitive systems such as Resting-State Networks (RSNs) and raising the issue of the specificity of this scale-free dynamics modulation in RSNs. Second, using analysis tools designed to measure a single scaling exponent related to the second order statistics of the data, thus relying on models that either implicitly or explicitly assume Gaussianity and (asymptotic) self-similarity, while fMRI signals may significantly depart from those either of those two assumptions (Ciuciu et al., 2008; Wink et al., 2008). To address these issues, the present contribution elaborates on the analysis of the scaling properties of fMRI temporal dynamics by proposing two significant variations. First, scaling properties are technically investigated using the recently introduced Wavelet Leader-based Multifractal formalism (WLMF; Wendt et al., 2007). This measures a collection of scaling exponents, thus enables a richer and more versatile description of scale invariance (beyond correlation and Gaussianity), referred to as multifractality. Also, it benefits from improved estimation performance compared to tools previously used in the literature. Second, scaling properties are investigated in both RSN and non-RSN structures (e.g., artifacts), at a broader spatial scale than the voxel one, using a multivariate approach, namely the Multi-Subject Dictionary Learning (MSDL) algorithm (Varoquaux et al., 2011) that produces a set of spatial components that appear more sparse than their Independent Component Analysis (ICA) counterpart. These tools are combined and applied to a fMRI dataset comprising 12 subjects with resting-state and activation runs (Sadaghiani et al., 2009). Results stemming from those analysis confirm the already reported task-related decrease of long memory in functional networks, but also show that it occurs in artifacts, thus making this feature not specific to functional networks. Further, results indicate that most fMRI signals appear multifractal at rest except in non-cortical regions. Task-related modulation of multifractality appears only significant in functional networks and thus can be considered as the key property disentangling functional networks from artifacts. These finding are discussed in the light of the recent literature reporting scaling dynamics of EEG microstate sequences at rest and addressing non-stationarity issues in temporally independent fMRI modes.

1. Introduction

Much of what is known about brain function stems from studies in which a task or a stimulus is administered and the resulting changes in neuronal activity and behavior are measured. From the advent of human electroencephalography (EEG) to cognitive activation paradigms in functional Magnetic Resonance Imaging (fMRI), this approach proved very successful to study brain function, and more precisely functional specialization in human brain. It has relied, on one hand, on contrasting signal magnitude between different experimental conditions (Rosen et al., 1998) or task-specific hemodynamic response (HRF) shape (Dale, 1999) and, on other-hand, on statistical methods often framed within linear or bilinear modeling strategies (Friston et al., 1995; Makni et al., 2005, 2008).

Spontaneous modulations of neural activity in Blood Oxygenation Level Dependent (BOLD) fMRI signals however arise without external input or stimulus and thus depict intrinsic brain activity (Damoiseaux et al., 2006). This ongoing activity constitutes a major part of fMRI recordings and is responsible for most of brain energy consumption. It has hence been intensively studied over the last decade using various methods ranging from univariate, i.e., Seed-based linear Correlation Analysis (SCA; Biswal et al., 1995; Greicius et al., 2003), to multivariate methods such as Independent Component Analysis (ICA; Calhoun et al., 2001; Beckmann and Smith, 2004), group-level ICA (Cole et al., 2010; Varoquaux et al., 2010b), or more recent dictionary learning techniques (Varoquaux et al., 2011). All these methods have revealed that interactions between brain regions, also referred to as functional connectivity, occur through these spontaneous modulations and consistently vary between rest and task (Damoiseaux et al., 2006; Fox et al., 2007). Resting-State Network (RSN) extraction from resting-state fMRI time series is thus achieved either by thresholding the correlation matrix computed between voxels or regions (seed-based or univariate approach) or by identifying spatial maps in ICA-based algorithms that closely match RSNs such as somato-sensory systems (visual, motor, auditory), the default mode, and attentional networks (ventral and dorsal; Fox et al., 2007; Smith et al., 2009). For a recent review about the pros and cons of the SCA and ICA approaches to RSN extraction, the reader can refer to Cole et al. (2010). Once RSNs are extracted, their topological properties can be analyzed with respect to small-world or scale-free models (Chialvo, 2004; Eguiluz et al., 2005; Zemanová et al., 2006; Bullmore and Sporns, 2009).

Small-world and scale-free topology led to model brain as a complex critical system, that is as a large conglomerate of interacting components, with possibly non-linear interactions (Bak and Paczuski, 1995; Chialvo, 2010). Further, these complex systems were then regarded as potential origins for long-range correlation spatio-temporal patterns, as critical systems, i.e., complex systems driven close to their phase transitions, constitute known mechanism yielding scaling time dynamics and generic 1/f power spectral densities (see e.g., Chialvo, 2010). They however so far failed to account for the existence of possibly richer scaling properties (such as, e.g., multifractality). At a general level, scale invariance in time dynamics and scale-free property of brain topology are, in essence, totally independent properties that must not be confused one with the other. Whether or not and how these two scale-free instances are related one to the other in the fMRI context remains a difficult and largely unsolved issue, far beyond the scope of the present contribution, that concentrates instead on performing a thorough analysis of scale invariance temporal dynamics in fMRI signals.

In the existing literature, the analysis of scale invariance in fMRI signals suffers from two limitations: First, it has often been performed at the voxel or region level, thus consisting of a collection of univariate analyses, suffering from the classical bias of voxel selection or region definition. Moreover, although the fluctuation of scale-free dynamics with tissue type has been studied in Shimizu et al. (2004), Wink et al. (2008) to derive that stronger persistency occurs in gray matter and that this background activity might represent neuronal dynamics, no systematic analysis has been undertaken to disentangle the scale-free properties of RSN and non-RSN components, such as artifacts. This investigation can be better handled using multivariate or ICA-like approaches. Second, scale invariance in fMRI signals has mostly been based on spectral analysis and/or Detrended Fluctuation Analysis (c.f., e.g., Thurner et al., 2003; Stam and de Bruin, 2004; He, 2011). This amounts to considering that scaling is associated only with the correlation or the spectrum (hence with the second order statistics) of the data and thus, implicitly and sometimes even explicitly, to assuming Gaussianity and (asymptotic) self-similarity for the data (cf., e.g., Eke et al., 2002) for a survey in the fMRI context). Also, it is now well-known that such technics lack robustness to disentangle stationarity/non-stationarity versus true scaling property issues and do not allow simple extension to account for richer scaling properties such as those observed in multifractal models. It is well accepted that wavelet analysis based analysis of scaling (cf., e.g., Abry et al., 1995, 1998; Bullmore et al., 2001; Veitch and Abry, 2001; Fadili and Bullmore, 2002) yield not only better estimation performance, but also show significant practical robustness, notably to non-stationarity, while paving the way toward the analysis of scaling properties beyond the strict second order (hence beyond Gaussianity and asymptotic self-similarity).

In this context, the present contribution elaborates on earlier works dedicated to the analysis of scale invariance in fMRI temporal dynamics by proposing two significant variations.

First, scale invariance dynamics is not investigated at the voxel or region spatial scale level independently. Instead, group-level resting-state networks are segmented by an exploratory multivariate decomposition approach, namely the MSDL algorithm (Varoquaux et al., 2011), detailed in Section 3: It produces both a set of spatial components and a set of times series, for each component and each subject, that conveys ongoing dynamics in functional networks but also in artifacts. As shown in Varoquaux et al. (2011), the sparsity promoting regularization involved in the MSDL algorithm enables to recover less noisy spatial maps than group-level or canonical ICA (Varoquaux et al., 2010b). This makes their interpretation easier in the context of small group of individuals. This technique is detailed in Section 3.

Second, to enable an in-depth analysis of the scaling properties of the temporal dynamics in fMRI signals, we resort to multifractal analysis, that measures not a single but a collection of scaling exponents, thus enabling a richer and more versatile description of scale invariance (beyond correlation and Gaussianity), referred to as multifractality. It is thus likely to better account for the variety and complexity of potential scaling dynamics, as already suggested in the context of fMRI in, e.g., Ciuciu et al. (2008), Wink et al. (2008). However, in contrast to Wink et al. (2008), and following the track opened in Ciuciu et al. (2008), we use a recent statistical analysis tool, the Wavelet Leader-based Multifractal formalism (WLMF; Wendt et al., 2007). This formalism benefits from better mathematical grounding and shows improved estimation performance compared to tools previously used in the literature. This framework is introduced in Section 4, after a review of the intuition, models, and methodologies underlying the definition and analysis of scaling temporal dynamics, thus, to some extend, continuing, and renewing the surveys provided in (Eke et al., 2002; Ciuciu et al., 2008).

These tools are combined together and applied to two datasets, corresponding to resting-state and activation runs. They are described in Section 2 (see also Sadaghiani et al., 2009). Modulations of scale-free and multifractal properties in space, i.e., between functional and artifactual components but also between rest and task, are statistically assessed at the group-level in Section 5.

In agreement with findings in He (2011), the results reported here confirm that fMRI signals can be modeled as stationary processes, as well as the decrease of the estimated long memory parameter under task. However, this is found to occur everywhere in the brain and not specifically in functional networks. Moreover, evidence for multifractality in resting-state fMRI signals is demonstrated except for non-cortical regions. Task-related modulations of multifractality appear only significant in functional networks and thus become the key property to disentangle functional networks from artifacts. However, in contrast to what happens for the long memory parameter, this modulation is not monotonous across the brain and varies between cortical and non-cortical regions. These results are further discussed in Section 6 in the light of recent findings related to scale-free dynamics of EEG microstate sequences and non-stationarity of functional modes. Conclusions are drawn in Section 7.

The rs-fMRI dataset we consider in this study has already been published in Sadaghiani et al. (2009). Eight hundred-twenty volumes of task-free “resting-state” data (with closed, blind-folded eyes) were acquired before getting experimental runs of 820 volumes each. These experimental runs, which have not been analyzed in Sadaghiani et al. (2009), involve an auditory detection task (run 2, motor response), and make use of a sparse supra-threshold auditory stimulus detection.

The auditory stimulus was a 500-ms noise burst with its frequency band modulated at 2 Hz (from white noise to a narrower band of 0–5 kHz and back to white noise). Inter-stimulus intervals ranged unpredictably from 20 to 40 s, with each specific interval used only once. Subjects were instructed to report as quickly and accurately as possible by a right-hand key press whenever they heard the target sound despite scanner’s background noise. Details about the definition of each subject’s auditory threshold are available in Sadaghiani et al. (2009).

2.2. Data Analysis

We used here statistical parametric mapping (SPM5, Wellcome Department of Imaging Neuroscience, UK2. For image preprocessing (realignment, coregistration, normalization to MNI stereotactic space, spatial smoothing with a 5-mm full-width at half-maximum isotropic Gaussian kernel for single-subject and group analyses) and our own software developments for subsequent analyses. More precisely, the MSDL algorithm relies on the scikit-learn Python toolbox3 and the multifractal analysis on the WLBMF Matlab toolbox4.

3. Multivariate Decomposition of Resting-State Networks

3.1. Multi-Subject Spatial Decomposition Techniques

The fMRI signal observed in a voxel reflects many different processes, such as cardiac or respiratory noise, movement effects, scanner artifacts, or the BOLD effect that reveals the underlying neural activity of interest. We separate these different contributions making use of a recently introduced multivariate analysis technique that estimates jointly spatial maps and time series characteristic of these different processes (Varoquaux et al., 2011). Formally, this estimation procedure amounts to finding K spatial maps Vs ∈ ℝp×K and the corresponding time series Us ∈ ℝn×K, whose linear combination fits well the observed brain signals, Ys ∈ ℝn×p, of length n, measured over p voxels, for subject s:

Ys=UsVst+Es,(1)

with Es ∈ ℝn×p the subject-level noise, or residuals not explained by the model. Finding Vst enables the separation of the contributions of the different process that are mixed at the voxel level, but implies to work on spatial maps rather than on specific voxels. The number of spatial maps, K, is not chosen a priori, but selected by the procedure.

This problem can be seen as a blind source separation task in the presence of noise, and has often been tackled in fMRI using ICA, combined with principal component analysis (PCA) to reject noise (McKeown et al., 1998; Kiviniemi et al., 2003; Beckmann and Smith, 2004). In the multi-subject configuration, estimating the spatial maps on all subjects simultaneously makes it easy to relate the factors estimated across the different subjects. This can be done by concatenating the data across subject, modeling a common distribution (Calhoun et al., 2001), or by extending the data-reduction step performed in the PCA by a second level capturing inter-subject variability (Varoquaux et al., 2010b). More recently, it was proposed that the key to the success of ICA on fMRI data, is to recover sparse spatial maps (Daubechies et al., 2009; Varoquaux et al., 2010a). This hypothesis can be formulated as a sparse prior in model (1), which can then be estimated using sparse PCA or sparse dictionary learning procedures. With regards to our goal in this study, extracting time series specific to the various processes observed, a strong benefit of such procedures is that they can perform data-reduction, i.e., estimation of the residuals not explained by the model, and extraction of the relevant signals in a single step informed by our prior. On the opposite, with ICA-based procedures, the residuals are selected by the PCA step, and not the ICA step.

3.2. Multi-Subject Dictionary Learning Algorithm

In addition, Varoquaux et al. (2011) have adapted the dictionary learning procedures to a multi-subject setting, in a so-called multi-subject dictionary learning (MSDL) framework. On fMRI datasets, the procedure extracts a group-level atlas of spatial signatures of the processes observed, as well as corresponding subject-level maps, accounting for the individual specificities. They show that, with a small spatial smoothness prior added to the sparsity prior on the maps, the extracted patterns correspond to the segmentation of various structures in the signal: functional regions, blood vessels, interstitial spaces, sub-cortical structures… In these settings, the subject-level maps Vs are modeled as generated by group-level maps V ∈ ℝp×K with additional inter-subject variability that appears as residual terms, Fs ∈ ℝp×K, at the group-level:

∀s∈1,…,S,Vs=V+Fs.

The model is estimated by finding the group-level and subject-level maps that maximize the probability of observing the data at hand with the given prior. This procedure is known as a Maximum A Posteriori (MAP) estimate, and boils down to minimizing the negated log-likelihood of the model with an additional penalizing term. If the two sources of unexplained signal, i.e., subject-level residuals Es and inter-subject variability Fs are modeled as Gaussian random variates, the log-likelihood term is the sum of squares of these errors. The prior term appears as the sum of the sparsity-inducing ℓ1 norm of V, and the ℓ2-norm of the gradient of the map, enforcing the smoothness. This prior has been used previously in regression settings under the name of smooth-Lasso (Hebiri and van de Geer, 2011). Estimating the model from the data thus consists of minimizing the following criterion:

𝒥Us,Vs,V=∑s=1SYs-UsVst2+μVs-V2+λV1+VtL V∕2

where, ||V||1 is the ℓ1 norm of V, i.e., the sum the absolute values, L is the image Laplacian −VtLV is the norm of the gradient. λ is a parameter controlling the amount of prior set on the maps, and thus the amount of sparsity, that is set by Cross-Validation (CV). μ is a parameter controlling the amount of inter-subject validation, that is set by comparing intra-subject variance in the observations with inter-subject variance. For more details about the estimation procedure or the parameter setting, we refer the reader to Varoquaux et al. (2011).

3.3. Resting-State MSDL Maps

rs-fMRI runs were analyzed for S = 12 subjects, consisting of n = 820 volumes (time points) with a 3-mm isotropic resolution, corresponding to approximately p = 50000 voxels within the brain. The automatic determination rule of the number of maps exposed in Varoquaux et al. (2010a) converges to K = 42. Also, the CV procedure gives us the best CV criterion for λ = 2. The group-level maps V are shown in Figure 1. They have been manually classified in three groups: Functional (F), Artifactual (A), and Undefined (U) maps that appear color-coded in red, blue, and green, respectively. The undefined class appeared necessary to introduce some confidence measure in our classification and disambiguate well-established networks (e.g., dorsal attentional network) from inhomogeneous components mixing artifacts with neuronal regions (e.g., like in v9). The anatomo-functional description of these group-level maps and their class assignment is given in Table 1. The same rules applied for individual maps Vs. In what follows, we will denote by ℱ, 𝒜 and 𝒰 the index sets of F/A/U-maps, respectively and by Card (ℱ) = 25, Card (𝒜) = 13, and Card (𝒰) = 4 their respective size.

FIGURE 1

Figure 1. From left to right and top to bottom, group-level MSDL maps V = |v1|…|v42| inferred from the multi-subject (S = 12) resting-state fMRI dataset (Neurological convention: left is left). Functional (F), Artifactual (A), and Undefined (U) maps appear color-coded boxes in red, blue, and green, respectively. Let us denote ℱ, 𝒜, and 𝒰 the index sets of F/A/U-maps, respectively and Card (ℱ) = 25, Card (𝒜) = 13, and Card (𝒰) = 4 their respective size. Each map vk consists of loading parameters within the (−1, 1) range where positive and negative values are depicted by the hot and cold parts of the color bar.

To compare spontaneous and evoked activity, the same spatial decomposition was used on resting-state (run 1, Rest) and task-related data, which were acquired during an auditory detection task (run 2, Task). In practice, this consists of projecting the task-related fMRI data Ỹs onto the inferred spatial maps Vs by minimizing the following least square criterion, Ỹs-WsVst2, with respect to Ws. The time series solution admits a closed-form expression: Ũs=ỸsVsVstVs-1. The subsequent scale-free analysis is applied to the two sets of n × K map-level fMRI time series US = [us,1|…|us,K]t and Ũs=ũs,1…ũs,Kt in a univariate manner, that is to each time series us,k and ũs,k for Rest and Task, respectively.

4. Scale-Free: Intuition, Models, and Analyses

4.1. Intuition

In the analysis of evoked brain activity, it is common to seek correlations of BOLD signals with any a priori shape of the hemodynamic response convolved with the experimental paradigm. In the frequency domain, this amounts to seeking response energy concentration in pre-defined spectral bands, as induced for instance by periodic stimulation (e.g., flashing checkerboards). In resting-state fMRI, it is now well admitted that intrinsic brain activity is characterized by scale-free properties (Zarahn et al., 1997; He, 2011). This constitutes a major change in paradigm as it implies that brain activity is not to be analyzed via the amounts of energy it shows within specific and a priori chosen frequency bands, but instead via the fact that all frequencies are jointly contributing in an equivalent manner to its dynamics. Scale-free dynamics are usually described in the spectral domain by a power-law decrease: Let Y(t) denote the signal quantifying brain activity and ΓY (f ) its Power Spectral Density (PSD). Scale-free property is classically envisaged as:

ℳ0:ΓYf≃Cf-β,β≥0,(2)

with fm ≤ |f | ≤ fM, fM/fm ≫ 1. Such a power-law behavior over a broad range of frequencies implies that no frequency in that range plays a specific role, or equivalently, that they are all equally important. To analyze brain activity, this power-law relation thus becomes a more important feature than the energy measured at some specific frequencies. For instance, it implies that energy at frequency f1 can be deduced from energy at frequency f2 according to He (2011):

ΓYf2=ΓYf1f2f1-β.(3)

In the scale-free framework, one therefore tries to quantify brain activity by considering the scaling exponent β (or variants) as the key descriptor. Let us moreover note that the terminology scale-free is equivalent to scale invariance or simply scaling, encountered in other scientific fields, where this property has also been found to play a central role (c.f., Abry et al., 2002; Ciuciu et al., 2008).

4.2. Scale-free Models

4.2.1. From spectrum to increments

Though appealing, equations (2) and (3) do not provide practitioners with a versatile enough definition of scale-free with respect to real-world data analysis. Indeed, they concentrate only on the second order statistics and hence account neither for the marginal distribution (first order statistics) of the signal Y, nor for its higher order dynamics (or dependence structure). For instance, it does not indicate whether data are jointly Gaussian or depart, weakly, or strongly, from Gaussianity.

with H = (−α/2) = (β + 1)/2, and where fd¯¯d means equality of all joint finite dimensional distributions: i.e., (X(t+τ1)-X(t))∕τ1H and (X(t+τ2)-X(t))∕τ2H have the same joint distributions. In turn, this implies that ∀q > −1, such that 𝔼|X(t)|q < ∞:

EXt+τ-X(t)q=cqτqH,τm≤τ≤τM,or(6)EXt+τ2-Xtq=EXt+τ1-Xtqτ2τ1qH,(7)

with τm ≤ τ1,τ2 ≤ τM, which are reminiscent of equations (2) and (3).

4.2.2. Self-Similar processes with stationary increments

Equations (6) and (7) turn out to hold not only for jointly Gaussian 1/f-processes but for a much wider and better defined class, that of self-similar processes with stationary increments, referred to as H-sssi processes, and defined as, c.f., Samorodnitsky and Taqqu (1994):

ℳ1:Xtt∈ℝ=fddaHXt∕at∈ℝ,(8)

∀a > 0, H ∈ (0, 1). Essentially, it means that X cannot be distinguished (statistically) from any copy, dilated by scale factor a > 0, on condition that the amplitude axis is scaled by aH. Parameter H is referred to as the self-similarity exponent. A major practical consequence of this definition consists of the fact that equations (6) and (7) hold for all τ (resp., τ1,τ2).

The central benefit of such a definition is that it does not require the data to be Gaussian but provides both theoreticians and practitioners with a well-defined model. For analysis, fMRI data can hence be envisaged as the increment process Y(t) = X(t + τ0) − X(t) of an H-sssi process X (where τ0 is an arbitrary constant chosen to make sense with respect to physiology and data acquisition set up, e.g., τ0 = TR). This constitutes a second model to account for scale-free properties in data, that encompasses the simpler 1/f-spectrum first model.

4.2.3. Multifractal processes

In a number of situations, it has been actually observed on a variety of real-world data of very different nature (c.f., e.g., Abry et al., 2002 for reviews) that equation (6) holds over a wide range of τs, however, with scaling exponents that depart significantly from the theoretical linear behavior qH:

EXt+τ-Xtq=cqτζ(q),τm≤τ≤τM.(9)

The generic behaviors modeled by equation (9) can be considered as a practical or operational, definition of scale-free property. Let us note that, by nature, ζ(q) is necessarily a concave function of q (c.f., e.g., Wendt et al., 2007).

Scaling exponents ζ(q) that are strictly concave rule out the use of H-sssi process as models. Instead, a broader class should be used, referred to as that of multifractal processes. This is however a large and not-well-defined class of processes. For the purposes of this contribution, let us use a particular subclass of multifractal processes defined as fBm subordinated to a multiplicative Compound Poisson cascade:

ℳ2:X(t):=BHA(t),whereA(t)=∫tW(s)ds,(10)

with W(s) a multiplicative Compound Poisson cascade (or martingale), such as those defined in Barral and Mandelbrot (2002). The complete definition of these cascades has been given and studied with details elsewhere and is hence not recalled here (c.f., Bacry et al., 2001; Barral and Mandelbrot, 2002; Chainais et al., 2005). It is enough to emphasize that they rely on the choice of positive random variables whose moments of order q define the ζ(q). The process X thus defined satisfies equation (9) with strictly convex tunable scaling exponents ζ(q), has stationary increments Y, and has distributions that depart from strict jointly Gaussian laws. Such departures, that may however turn subtle and hard to detect in practice, are precisely quantified by the departure of ζ(q) from a linear behavior in q. The ζ(q) therefore convey a rich information about data X, and hence about Y, as they account for the entire dependence structure of the data, hence both to the time dynamic and distributions of data. Their accurate estimation from real-world data therefore naturally constitutes an important practical challenge discussed below.

4.3. Scale-Free Analysis

4.3.1. From spectrum to wavelet analysis

Assuming that data Y have a power-law spectrum behavior as in equation (2), it is natural to rely on spectral estimation to measure β. A classical tool in spectrum analysis is the Welch estimator that consists in splitting data Y into blocks and in averaging the squared Fourier transforms computed independently over each block. For scale-free data, it is hence expected that:

Γ^Yf=∑kY,gf,k2≃Cf-β,(11)

where the gf,k = g0(t − k)e℩2πft are translated into time and into frequency templates of a reference pattern g0(t). This relation can be further used to estimate β.

It has been shown that wavelet transforms can achieve better performance both in the analysis of scale-free properties in real-world data, and in the estimation of the corresponding scaling parameters (c.f., Abry et al., 1995, 1998; Veitch and Abry, 2001). The discrete wavelet transform (DWT) coefficients of Y are defined as:

dYj,k=∫ℝYt2-jψ02-jt-kdt≡Y,ψj,k,(12)

where the ψj,k = 2−jψ0(2−jt − k) consists of templates of a reference pattern ψ0 translated in time and dilated (by a factor a = 2j). It is referred to as the mother-wavelet: an elementary function, characterized by fast exponential decays in both the time and frequency domains, as well as by a strictly positive integer Nψ ≥ 1, the number of vanishing moments, defined as ∀k=0,1,…, Nψ − 1, ∫tℝtkψ0(t)dt ≡ 0, and ∫tℝtNψ0(t)dt ≠ 0. Note the choice of the L1-norm (as opposed to the more common L2-norm choice) that better matches scaling analysis. For further introduction to wavelet transforms, the reader is referred to, e.g., Mallat (2009).

Defining SYd(j,2)=1njΣk=1njdY(j,k)2 (with nj the number of dX(j,k) available at scale 2j), one obtains (c.f. Abry et al., 1995):

ESYdj,2=∫ℝΓYfΨ02jf2df(13)

where Ψ0 denotes the Fourier transform of ψ0. This indicates that SYd(j,2) can be read as a wavelet based estimate of the PSD and is hence referred to as the wavelet spectrum. It measures the amount of energy of Y around the frequency fj = f0/2j where f0 is a constant that depends on the explicit choice of ψ0 (for the Daubechies wavelet used here, f0 ≃ 3fs/4 with fs the sampling frequency). This correspondence between the Fourier and wavelet spectra is illustrated on fMRI signals in Figure 2. For scale-free processes satisfying equation (2), it implies:

While this formally looks like equation (11), it has been shown in detail how and why the wavelet spectrum yields better estimates of the scaling exponents β than Welch based-ones, both in terms of estimation performance and robustness to various forms of non-stationarity in data that may be confused with scale-free behaviors (Abry et al., 1995, 1998; Veitch and Abry, 2001). Notably, it was shown how wavelet analysis enables to disentangle non-stationarity, stemming from fMRI environment, from true long memory in brain activity. Also, the wavelet spectrum avoids the potentially difficult issue that consists of deciding a priori whether empirical data are better modeled by Y or X, needed by classical spectrum estimation, that can only be applied to stationary data. In a nutshell, these benefits stem from the use of the change of scale operator to design the analysis tool, that intuitively matches scale-free behavior more naturally than a frequency shift operator.

4.3.2. From 2nd to other statistical orders: wavelet leaders

As discussed in Section 4.2, analyzing in-depth scale-free properties implies investigating not only the spectrum (i.e., the second order statistics of data) but rather the entire dependence structure, i.e., the whole range of available statistical orders q. It had initially been thought that this would amount to extending the definition of SYd(j,2) to other orders q, SYd(j,q)≡1nj∑k=1njY,ψj,kq. It has however recently been shown that this approach, though intuitive and appealingly simple, fails to yield satisfactory estimation of the ζ(q). Notably, wavelet coefficients show little power in enabling practitioners to decide whether ζ(q) is a linear or strictly concave function of q. Instead, it is now well documented that the estimation of the ζ(q) should be based on Wavelet Leaders (Wendt et al., 2007).

Let us now assume that ψ0 has a compact time support and introduce the global regularity of Y, hm, defined as: hm=lim⁡inf⁡2j→0log⁡(sup⁡k|dY(j,k)|)/log⁡(2j) Therefore, hm can be estimated by a linear regression of the log of the magnitude of the largest wavelet coefficient at scales 2j versus the log of the scales 2j (Wendt et al., 2007). Let γ ≥ 0 be defined as, with ∈ > 0: γ = 0 if hm > 0, and γ = −hm + ∈ otherwise. Further, let λj,k denote the dyadic interval λj,k = [k2j, (K + 1)2j), and denote by 3λj,k the union of λj,k and its 2 closest neighbors, 3λj,k = [(k − 1)2j, (k + 2)2j). The wavelet leaders LY(γ) are defined as LY(γ)(j,k)=supλ′⊂3λj,k2γjdY(λ′). In practice, LY(γ)(j,k) simply consists of any of the largest coefficients 2γj|dY(λ′)| located at scales finer or equal to 2j and within a small time neighborhood. It is then necessary to form the so-called wavelet Leader structure functions that reproduce the scale-free properties in Y according to:

SYLj,q,γ≡1nj∑k=1njLY(γ)j,kq≃cq2jζq,γ,(14)

Moreover, for a large class of processes, one has: ζ(q,γ) = ζ(q) + γq. For all real-world data analyzed so far with WLMF, this relation is found to hold, by varying γ (c.f. Wendt et al., 2007 for a thorough discussion). This has also been verified empirically for fMRI data. Further, because it can take any concave shape, the function ζ(q,γ) is often written as a polynomial expansion (Arneodo et al., 2002): ζ(q,γ)=Σp≥1cp(γ)qp∕p!. Notably, the second order truncation ζ(q,γ)≃c1(γ)q+c2(γ)q2∕2 (with c2(γ)≤0 by concavity) can be regarded as a potentially interesting approximation that captures the crucial information regarding whether the ζ(q,γ) are linear in q (hence indicating H-sssi models) or strictly concave (hence suggesting multiplicative cascade models). Interestingly, the coefficients cp(γ) entering the polynomial expansion of ζ(q,γ) are not abstract figures but rather turn out to be quantities deeply tied to the scale-free properties of Y, as they are related to the scale dependence of the cumulants of order p ≥ 1, CY(γ)(j,p), of the random variable lnLY(γ)(j,k):

This wavelet Leader-based analysis of scale-free properties is intimately and ultimately related to multifractal analysis, the detailed introduction of which is beyond the scope of the present contribution. We restate here only its essence. Multifractal analyses describe globally the fluctuations along time of the local regularity of a signal Y(t). This local regularity is measured by the so-called Hölder exponent h(t), that essentially compares Y around time t0 against a local power-law behavior: |Y(t) − Y(t0)| ≤ |t − t0|h, |t − t0| → 0. The variations of h along time are then described globally via the multifractal spectrum, consisting of the collection of Hausdorff dimensions, 𝒟(h), of the sets of points {t,h(t) = h}. In practice, the multifractal spectrum is estimated indirectly via (a Legendre transform of) the function ζ(q). The approximation ζ(q) ≃ c1q + c1q2/2 translates into 𝒟(h) ≃ 1 − (h − c1)2/(2|c2|). For thorough and detailed introductions to multifractal analysis, the reader is referred to, e.g., Wendt et al. (2007). Examples of such multifractal spectra estimated using the WLMF from real fMRI signals are illustrated in Figure 2B. An outcome of the mathematical theory underlying multifractal analysis, of key practical importance and impact, is the following: the function 𝒟(h) theoretically constitutes a rich characterization of the scale-free properties of a signal Y and its complete and entire estimation requires the use, in equation (14), of both positive and negative order qs, concentrated left and right around 0 (Wendt et al., 2007).

5. Multifractal Analysis of MSDL Maps

5.1. Single-Subject Analysis

5.1.1. Scaling range

For analysis, orthonormal minimal-length time support Daubechies’s wavelets were used with Nψ = 3. Scale-free properties are systematically found to hold within a 4-octave range ((j1,j2) = (3,6)), corresponding to a frequency range of [0.008, 0.063] Hz5, which is hence consistent with the upper limit 0.1 Hz classically associated with the hemodynamics boundary and scaling in fMRI data (Cordes et al., 2001).

5.1.2. Fourier vs. wavelet spectra

For illustrative purposes, two time series corresponding to a functional map (k = 28, f18 in Table 1), were selected in the rest and task runs from the first subject. In Figure 2A, the Fourier spectrum estimate (log2Γ^us,k(f)) based on Welch’s averaged periodogram and its wavelet spectrum counterpart (log2Sus,kd(j,2)) are found to closely match, as predicted by equation (13). Interestingly, Figure 2A shows that the β exponent, measured within frequency range [0.008, 0.063] Hz, in equation (2) (i.e., the neg-slope of the log-spectra log2Γ^us,k(f)) decreases with task-related activity in f18. This amounts to observing lower Hurst exponent H = (β − 1)/2 in the task-related dataset: H^f18R≃0.66 and H^f18T≃0.5. As shown in the following, this decrease of self-similarity is not specific to functional maps and will be observed in artifactual and undefined maps. Following He (2011), the stationarity of fMRI signals is confirmed since we systematically observed H^kR,T<1.

5.1.3. Multifractal spectrum

For the same time series, MF spectra 𝒟(h), estimated using the WLMF tool described above, are depicted in Figure 2B. The decrease of self-similarity between rest and task is captured by a shift to the left of the position c^1 of the maximum of D(h): ((c^1)f18R,(c^1)f18T)=(0.75,0.5) It should also be noted that parameter c^1 systematically takes values that are close to those of the Hurst exponent. This is consistent with the theoretical modeling of scale-free property that establishes a clear connection between c1 and H and predicts c1 ≃ H (c.f. Wendt et al., 2007). Therefore, in the following, c1 will be referred to as the self-similarity parameter although this is a slight misnomer. Further, Figure 2B confirms the presence of multifractality in fMRI data as strictly negative c2 < 0 are almost always observed. Indeed, parameter c2 quantifies the width of 𝒟(h; as a curvature radius of 𝒟(h) around (c^1): c^2<0. Multifractality is however not specific to a given brain state since we measured ((c^2)f18R,(c^2)f18T)=(-0.07,-0.06). In this example, multifractality, as measured by the width of the multifractal spectra, is decreased from rest to task. However, opposite fluctuations will be also observed amongst F-maps.

The sole two self-similarity and multifractality parameters c1 and c2 are therefore used from now on as sufficient and relevant descriptors of the scale-free properties of fMRI signals (superscript γ is dropped for the sake of conciseness, while γ has been systematically set to γ = 2).

5.2. Group-Level Analysis

5.2.1. Group-level scale-free properties

Let ci,kj,s denote the ĉ1 and ĉ2 estimates (index i = 1:2) for differents maps (index k = 1:K), runs (index j = R,T for Rest and Task, respectively) and for different subjects (index s). The map-dependent group-level values have been computed as μi,kj=Σs=1Sĉi,kj,s∕S and sorted according to their labeling (F/A/U-maps) given in Table 1. Then, global spatial averaging of the means μi,kj has been performed so as to derive global F/A/U-average parameter estimates: μ̄i,Fj=Σk⊂ℱμi,kj∕Card(ℱ), μ̄i,Aj, and μ̄i,Uj are defined equivalently. In the same spirit, group-level multifractal attributes μ̄i,vℓj are derived for each functional network vℓ ∈ 𝒩 = {Att, DMN, Mot, N-c, Vis} such that μ̄i,nℓj=Σk∈nℓμi,kj∕Card(nℓ), ∀ℓ = 1:5, and j = (R,T). We proceed in the same way for analyzing artifact types tr ∈ 𝒯 = {Ven, WhM, Mov, Oth}, and computing μ̄i,trj for r = 1:4.

As shown in Figure 3[top], the group-averaged values of self-similarity μ1,kj lie approximately in the same range [0.55, 1], indicating long memory, for all components (F/A/U-maps). An almost systematic decrease of self-similarity is observed in the task-related dataset (δ1,k=μ1,kT-μ1,kR<0), for k ∈ ℱ∪𝒜∪ 𝒰. This trend is therefore not specific to F-maps. Moreover, the average decrease computed over F-maps is about the same as the one estimated for A and U-maps (δ¯1,F=−0.125,δ¯1,A=−0.11and δ¯1,U=−0.13). Also, the averaged standard deviations (σ¯1,FR,σ¯1,AR,andσ¯1,UR) computed over the F/A/U-maps, are close to each other (σ̄1,F∕A∕UR≈0.18) and systematically increase with the task-related activity (σ̄1,F∕A∕UT>σ̄1,F∕A∕UR).

Figure 3[bottom] illustrates that the group-averaged values of μ2,kR,T are almost all negative in the F/A/U-maps indicating multifractality in fMRI time series irrespective of the map type or brain state. Between rest to task-related situation minor changes in the A and U-maps are also observed since |δ2,k| < 0.03 for k∈𝒰∪𝒜) while we measured |δ2,k| < 0.08 for k∈ℱ (δ2,k=μ2,kT-μ2,kR). Hence, the level of multifractality does not change much between rest and task in irrelevant maps. In contrast, large changes in the multifractal parameters are observed in F-maps, while not systematically in the same direction. For instance, in cerebellum ( f4), basal ganglia ( f5), DMN ( f7) and fronto-parietal network ( f8) evoked activity induces a large increase of multifractality (δ2,k< 0) while in the auditory and attentional systems (e.g., f12 and f24, respectively), which are supposed to be involved in the auditory detection task, the converse observation holds, i.e. (δ2,k> 0. Also, it is worth noticing that the averaged standard deviations computed over the A/U-maps increase when switching from rest to task (σ̄2,AR=0.06<σ̄2,AT=0.09andσ̄2,UR=0.05<σ̄2,UT=0.085) while they remain at the same level in the F-maps: σ̄2,FR≈σ̄2,FT≈0.08.

We computed the grand means of the self-similarity parameters μ̄1,F∕A∕UR,T over the F/A/U-maps, respectively, and draw the same conclusion at this macroscopic level, as demonstrated in Figures 4A–C: the decrease of self-similarity from rest to task is not specific to functional components and only slightly fluctuates between networks and artifact types. Moreover, we did not observe any significant modification of the grand means of multifractal parameter estimates μ̄2,F∕A∕UR,T between rest and task, as illustrated in Figure 4D. This motivated deeper investigations at the network and artifact levels, especially concerning the fluctuation of multifractality induced by task. Figure 4E reveals that a major increase of multifractality (μ¯2,n4T<μ¯2,n4R) occurred only in the non-cortical regions while no major change appeared in the artifacts (μ2,tr−T≃μ2,tr−R,∀r∈𝓣) as shown in Figure 4F.

5.2.2. One-sample statistical tests

To assess the statistical significance of the multifractal parameters for the rest and task-related datasets at the group-level, we used one-sided tests associated with the following null hypotheses ∀k ∈ ℱ∪𝒜 ∪ 𝒰:

We also conducted similar tests at the macroscopic level (k ∈ 𝒩 ∪ 𝒯) by replacing μi,kj with μ¯i,kj in the null hypotheses (16). Because there is no definite proof nor evidence that MF parameter estimates ĉi,kj,s should be normally distributed across subjects, we investigated different statistics (Student-t, Wilcoxon’s signed rank (WSR) statistic). Indeed, other statistics may provide more sensitive results in presence of outliers. To account for multiple comparisons (K tests performed simultaneously) and to ensure correct specificity control (control of false positives), the Bonferroni correction was applied.

Rejecting H0,j(1,k) clearly amounts to localizing brain areas or components eliciting significant long memory or self-similarity. Rejecting H0,j(2,k) enables to discriminate multifractality from self-similarity. Similar tests involving μ¯i,F/A/Uj,μ¯i,nℓj, and μ¯i,trj in the definition of null hypotheses(16) for (i = 1, 2) were also performed.

Analysis of statistical significance of F-maps regarding H0,R(1),k showed that most components (22/25) rejected this null hypothesis at rest using T-test and thus were significantly self-similar (see blue curves in Figure 5A). The task effect then induced a loss of significance in the vast majority of components as shown in Figure 5B: only four maps ( f10, f14, f18, and f24) demonstrated a significant level of self-similarity using T-test in the task-related dataset. These maps are related to the motor, fronto-parietal, and attentional (parieto-temporal junction and IPS/FEF) networks. Two out of them are lateralized in the left hemisphere. Statistical analysis of F-maps regarding H0,R(2),k demonstrated that only six components ( f15, f17, f18, f21, f23, f24) rejected this null hypothesis at rest: see red curves in Figure 5A. The task-related modulation tends to reduce the number of significant F-maps: As depicted in Figure 5B, only 3 components survived the T-test ( f10, f15, and f19) in the task-related dataset. Interestingly, f10 and f19 are likely to be involved in the auditory detection task and the motor response since they belong to the Motor and Attentional networks. Hence, a significant level of multifractality is observed during task in components that were monofractal at rest. Besides, the level of multifractality remains significant in the ventral occipital cortex ( f15) irrespective of the brain state and that a few components in the visual ( f17), fronto-parietal ( f18), temporal ( f21), prefrontal ( f23), and attentional ( f24) networks became monofractal under the task effect.

Statistical analysis of A and U-maps regarding H0,j(1),k showed the same behavior when switching from rest to task, namely a strong decrease of the number of significant self-similar components (from 10 to 4 and 4 to 2 for A/U-maps, respectively): see blue curves in Figures 5C–F, respectively. Statistical analysis of A and U-maps regarding H0,j(2),k also demonstrated a reduction of the number of multifractal components in A/U-maps. Two artifactual components (a10 and a12) located in the white matter remained consistently multifractal in both datasets and one undefined component (u3) became significantly multifractal when switching from rest to task. In all cases, a loss of significance is observed using WSR tests (dash dotted curves) instead of T-tests (solid curves) indicating that there is no outlier in this group and thus that the Gaussian distribution hypothesis is tenable.

Then, we focused on the statistical analysis at different macroscopic scales, first by averaging all F/A and U-maps respectively so as to derive a mean behavior for F/A/U-maps. Finally, we looked at functional networks and artifact types in more details. Blue curves in Figures 6A,B report such results for the rest and task-related datasets, respectively. We still observed a significant level of self-similarity in all averaged groups (blue curves) irrespective of the brain state: H¯0,j(1,F/A/U) is systematically rejected for j = (R, T). However, we still noticed a reduction of statistical significance induced by task irrespective of the map type. More interestingly, we found at this macroscopic level that all averaged maps were multifractal at rest whereas only the functional one remained multifractal during task: see red curves in Figures 6A,B. Further, statistical analysis of functional networks defined in Table 1 was conducted to understand which network drives this effect. When comparing p-values in Figures 6C,D on functional networks, we observed that all remained significantly self-similar in both states, while the DMN is close to the significance level α = 0.05 during task (blue curves). Regarding multifractality, only the non-cortical regions appeared monofractal at rest and all networks kept a significant amount of multifractality during task. In contrast, this observation did not hold for artifacts: when looking at Figures 6E,F in detail, the signal related to ventricles became monofractal during task.

5.2.3. 2-way repeated measures ANOVA

In order to assess any significant change of self-similarity or multifractality between rest and task, we entered the subject-dependent parameter estimates (ĉi,kj,s) in several 2-way repeated measures ANOVAs involving two factors: brain state (two values: j = R, T) and map type (with varying number of values). These ANOVAs were conducted separately for assessing self-similarity (i = 1) and multifractality (i = 2) changes. First six ANOVAs (three for each parameter) were carried out by considering the F/A/U-maps as the second factor, respectively. This second factor thus took a number of values that depends on the set under study: ℱ, 𝒜, or 𝒰. Results are summarized in Table 2. Regarding the analysis of self-similarity (ĉ1,kj,sparameters), a significant brain state effect appeared in all F/A/U-maps, and a significant map effect in the F and A-sets. Significant interactions were found for the F and U-maps. This confirms that the level of self-similarity is not sufficient to disentangle functional networks from artifactual or undefined maps.

As regards ANOVAs based on ĉ2,kj,s parameters, a significant interaction for F-maps is found, thus indicating that the averaged change in multifractality between rest and task is significant for functional maps only. In summary, only F-maps exhibited significant interactions for both multifractal attributes.

Akin to the one-sample analyses above, we looked at a larger spatial scale, the functional network, and artifact type levels and performed similar ANOVAs, corresponding results are reported in Table 3. While both functional networks and artifacts demonstrate a significant change in the self-similarity parameter between rest and task, only functional networks made the map type effect significant. More importantly, the key feature for discriminating functional networks from artifacts relied on ANOVAs based on ĉ2,kj,s parameters. Indeed, a significant network effect and more importantly a significant interaction between rest and task are observed in functional networks.

5.2.4. Two-sample statistical tests

To localize which maps are responsible for statistically significant ANOVA results, we finally performed two-sample T-tests in which we tested the following null hypotheses:

H̃0(1,k):μ1,kR=μ1,kT,∀k∈ℱ∪𝒜∪𝒰H̃0(2,k):μ2,kR=μ2,kT,∀k∈ℱ∪𝒜∪𝒰.(17)

We also conducted similar tests at the macroscopic level (k ∈ 𝒩∪𝒯) by replacing μi,kj with μ̄i,kj in the null hypotheses (17). The fluctuations in self-similarity being systematically in the same direction between rest and task, we performed one-sided tests as regards the μ1,kj’s while two-sided tests were considered for the μ2,kj’s: task-related positive and negative fluctuations of μ2,kj were actually observed in Subsection 5.2.1. Figures 7A,B shows the uncorrected p-values for the F-maps and networks, respectively. We rejected H̃0(1,k) for ( f3, f4, f11, f18, f25) at a significance level set to α1 = 0.01 and H̃0(2,k) for ( f4, f7, f18) at α2 = 0.05. These components clearly explain significant results reported in Table 2 about the changes in self-similarity and multifractality that occurred in F-maps. Interestingly, among the latter, the null hypothesis was rejected because of a large increase of multifractality in ( f4, f7). In contrast, a decrease of multifractality was responsible for the rejection of H̃0(2,k) in f18. When setting α2 = α1 = 0.01, only f18 survived this threshold and thus remained the single functional component for which a significant difference of self-similarity and multifractality was found between rest and task. This component clearly drove the significant interaction reported in Table 2 for the change in multifractality in F-maps. Figure 7B also showed that the state effect reported in Table 3 on (ĉ1,kj,s) at the network level was driven by the attentional, motor, and visual systems. Last, the significant interaction reported in Table 3 on (ĉ2,kj,s) is explained by the non-cortical regions as shown in Figure 7B.

FIGURE 7

Figure 7. Uncorrected p-values associated with two-samples Student-t test performed between resting-state and task-related multifractal parameters for assessing H̃0(1,⋅) (blue curves) and H̃0(2,⋅) (red curves) on the on the functional (A), artifactual (C) and undefined maps (E), respectively. Similar tests were computed at a broader spatial scale on the functional networks (B), the artifacts (D) and the averaged map types (F), respectively. Significance levels (α1 = 0.01 and α2 = 0.05) are shown in - - and , respectively.

Figures 7C,D shows the localization of the state effects reported in Tables 2 and 3 for the changes in self-similarity that occurred in artifacts at the local and global levels. No A-map enabled to reject H̃0(1,k) at the α1 significance level but a majority of A-maps (a1:4, a6, a8, a10, a12) contributed to the significant state effect observed in Table 2. At the global artifact level, the ventricles appear as the main source of the significant state effect reported in Table 3 for the change in self-similarity. Also, no significant difference in multifractality was reported for artifacts whatever the observation level (A-maps or averaged artifacts). Similarly, Figure 7E enables us to show that u2 and u4 were the main sources of the significant state effect and interactions reported in Table 2 for the change in self-similarity. At the macroscopic level, we finally observed in Figure 7F that only the grand mean of functional maps leads to a significant modulation of self-similarity between rest and task at level α1.

6. Discussion

6.1. Results Interpretation

This study analyzed in-depth the scale-free properties of fMRI signals, using multifractal methodologies, and their modulations during rest and task both in functional networks and artifactual regions. The underlying goal was to finely characterize which properties are specific to functional networks and which modulation can be expected for these networks from task-related activity. Previous attempts in the literature (Cordes et al., 2001; Leopold et al., 2003; He, 2011) focused on functional networks without comparing results with the behavior of artifacts. The main reason comes from the fact that seed region analyses were only conducted in such studies. Hence, no comparison with vascular or ventricles-related signals was undertaken.

Our results confirmed that fMRI signals are stationary and self-similar but not specifically in functional networks. Also we showed that the amount of self-similarity significantly varies between rest and task not only in functional networks involved in our auditory detection task with a motor response (Attentional, Motor) but also quite surprisingly in the visual system and in some artifacts (ventricles) and undefined maps. This observation led us to investigate the scale-free structure of fMRI signals using richer models, namely multifractal processes, to which the WLMF toolbox is dedicated. Our statistical results demonstrate first that fMRI signals are multifractal, second that interactions between brain state and maps only occurred in F-maps and functional networks and third, that specific F-maps such as in non-cortical regions demonstrated a statistically significant fluctuation between rest and task. This result shows that the concept of multifractality permits to disentangle functional components from artifactual ones, in a robust and significant manner.

However, in contrast to self-similarity that systematically decreases with evoked activity, multifractality decreases in cortical ( f18) but increases in non-cortical ( f4, f7). Thus, task-related activity has no systematic impact with respect to increase/decrease of multifractality. Interestingly, we found a statistically non-significant trend toward a decrease of multifractality in regions primarily involved in the task ( f12, f24, f25). However, the group size of this study remains small (12 subjects only) to achieve significant results, mainly because of the between-subject variability and of the difficulty in estimating ĉ2,kj,s parameters on short time series.

Further investigations beyond the scope of this paper are necessary to find out any general trend on the direction change of multifractality with evoked activity by cross-correlating multifractal parameters with task-related activity (e.g., group-level Z-scores) and task performance. However, to derive reliable results for multifractality, a larger group of individuals will be considered and a larger number of scans will be acquired while maintaining the same scanning time: To this end, accelerated SENSE imaging will be used together with recent reconstruction algorithms so as to improve temporal resolution (Chaari et al., 2011b).

The results obtained in this contribution shows multifractal temporal dynamics in fMRI signals and thus naturally lead to question the potential origins and generative mechanisms for this departure from the more traditional long-range correlation modeling of scale invariance. A natural track to inspect consists of that of the relations between hemodynamic (fMRI) and electrical (EEG) signatures for brain activity at rest. This question has been intensively studied over the last decade (Laufs et al., 2003; Mantini et al., 2007; Britz et al., 2010; Musso et al., 2010; Van de Ville et al., 2010; Yuan et al., 2012), first by measuring cross correlations between fMRI data at rest and EEG-informed regressors derived from the convolution of the EEG power signal in five well-identified frequency bands (δ ∈ (1,4) Hz, þeta ∈ (4,7) Hz, α ∈ (8,12) Hz, β ∈ (13,30) Hz, and γ > 30 Hz) with the canonical HRF. This approach revealed the negative correlation of α-band activity with the attentional network and the positive correlation with β2-band with the default mode network (precuneus and posterior cingulate cortex; Laufs et al., 2003). Also, Mantini et al. (2007) showed that functional resting-state networks have different EEG signatures which are not specific to a given frequency band but are rather spread over several oscillations regimes (e.g., correlation between α and β power in specific RSN), a consequence of the so-called oscillation hierarchy (Buzsáki, 2006) and the of phase-amplitude cross-frequency coupling (He et al., 2010). However, none of these works enable to explain the low frequency fluctuations (<0.1 Hz) or scale-free dynamics of the fMRI signal at rest, because this phenomenon is much more widespread than oscillations.

Scale-free dynamics of brain electrical activity at rest has recently been studied (Van de Ville et al., 2010) but not directly on raw data. Instead, EEG microscates that correspond to short periods (100 ms) during which the EEG scalp topography remains quasi-stable, have been first segmented. Remarkably, it has been shown that only four different EEG microstate patterns are necessary to describe the ongoing electrical brain activity at rest (Britz et al., 2010) and that these four microscates correlate with well-known RSNs, which were classically identified from fMRI dataset alone using group-level ICA. This demonstrated that the EEG microstate with rapid fluctuations might be considered as the electrophysiological signature of intrinsic functional connectivity patterns. The investigation of scale-free dynamics was thus performed on the EEG microstate sequence to understand how fast the microscates are changing and what kind of correlation structure (short or long-range) they bring (Van de Ville et al., 2010).

The recent finding that EEG microscate sequences reveal purely monofractal dynamics (Van de Ville et al., 2010), irrespective of the data filtering, may lead to conclude that the same monofractal behavior in the fMRI signature of RSN (strongly correlate with these microstates) should be expected, if one assumes a linear and time invariant HRF model for the neurovascular coupling. However, the results obtained in the present contribution can be considered not only as evidence in favor of multifractality in fMRI data, but also as evidence that this multifractal effect is discriminant of cortical versus non-cortical regions and characteristic of functional network with respect to modulation under task.

Several factors may explain this apparent discrepancy. First, an accurate comparison of both sets of result would require a precise match of the range of scales (or frequencies) within which scale invariance is analyzed and corresponding parameters measured. Here, the selected range of frequencies corresponds to ([0.008,0.063]Hz), while the monofractal behavior of EEG microstate sequences was exhibited on a distinct frequency range, i.e. ([0.063, 3.9]Hz) in Van de Ville et al. (2010). Comparison of scaling properties requires that the same frequency range is selected but this constraint is clearly not tenable across modalities like EEG and fMRI given the fMRI sampling rate.

Second, it is indeed very unlikely that a linear and time invariant filtering may create multifractality in fMRI starting from a monofractal electrophysiological signal in EEG. The general issue of the relations between (linear and non-linear) filtering and multifractality were barely studied theoretically so far but interestingly, Abry et al. has shown that simple non-linear filter can turn mono- into multifractality. Hence, another putative origin for the apparent contradiction between our findings and those in Van de Ville et al. (2010) lies in refined descriptions of HRF model by non-linear dynamical systems (e.g., Balloon model; Buxton et al., 1998, 2004). Of course, linear and stationary approximations like the canonical HRF model (Glover, 1999) or non-parametric alternatives (Vincent et al., 2010; Chaari et al., 2011a) have been validated but only on evoked activity and considering inter-stimulus intervals larger than 3 s. For shorter ISIs, non-linear hemodynamics has turned out to be a valid property (Liu and Gao, 2000). In this context, habituation, or repetition suppression effects may occur and induce a sublinear hemodynamic response, which would modify scaling properties (Dehaene-Lambertz et al., 2006; Ciuciu et al., 2009). Hence, by modeling the sequence of EEG transient brain states as a series of short time epochs, this could induce non-linearities in the hemodynamic system that could explain the switch from purely fractal EEG microstates to multifractal signatures in the corresponding RSNs.

Third, instead of segregating EEG microstates in multiple groups based upon the maximal spatial dissimilarity between groups (Britz et al., 2010; Musso et al., 2010), a more recent analysis of joint EEG/fMRI resting-state data has revealed a larger number (thirteen) of EEG microstates that show temporal independence from each other (Yuan et al., 2012). In this latter work, all resting-state networks including visual, motor, auditory, attention, saliency and default mode networks were characterized by a specific electrophysiological signature involving several EEG microstates. This clearly indicates that the original analysis of scale-free dynamics for EEG microscates done in Van de Ville et al. (2010) should be revisited on this larger number of metastable states to disentangle whether multifractality in this larger set of microstates has been discarded due to averaging effects. It is actually clear that the sequence mixing thirteen different microstates may generate richer singularities (abrupt changes between microstates) than the ones relying on four microstates only. Fourth, the temporal signatures of EEG microstates found in Musso et al. (2010), Van de Ville et al. (2010) are correlated in time since the spatial similarity was the key factor to identify them. As a consequence, the microstate sequences is correlated too and might loose some singularities that could be found out in the microstate sequences generated by Yuan et al. (2012). Finally, the presence of multifractality in resting-state (and task-related) MEG data has been evidenced in the sensor space in Zilber et al. (2012). These findings open new research avenues: For instance, it is natural to explore whether the observed multifractal properties can be related multiplicative cascade processes, that is to one of the only practical mechanism known to generate multifractal dynamics, or to investigate whether this cascade takes place at meso or macroscopic scales, as well as to figure out how brain networks could implement such cascade mechanisms. This topic is beyond the scope of the present contribution, however the log-normal statistics of neuronal firing rate could provide us with a first clue to uncover any generative process underlying multifractal dynamics.

6.3. Stationarity vs. Non-Stationarity of the RSN Dynamics

Recent results in resting-state fMRI reveal temporally independent functional modes of spontaneous brain activity (Smith et al., 2012) and postulate the presence of temporally non-stationary modes in part of the default mode network by resorting to high temporal resolution fMRI. While stationarity receives a unique and clear definition, non-stationarity can correspond to a bunch of different situations; for example, non-stationarity might (i) refer to an apparent change over time in the correlation between two regions or (ii) refer to changes in the mean and/or variance in the time course of a functional network.

The wavelet based analysis of scaling proposed here already addresses a number of such situations. The fact that the estimated Hurst coefficient of fMRI time series remains consistently below 1 indicates that fMRI signals at hand here are better modeled as a stationary step process Y rather than as a non-stationary random walk X. Further, wavelet analysis are known to bring robustness against various forms of non-stationarities, such as smooth trends superimposed to data, to mean or variance modulation (c.f. (ii)). The multifractal analysis performed here is thus not impaired by such form of non-stationarities. This leaves open issues such as the presence of oscillations superimposed to scaling. Given that time series are very short, the use of formal stationarity test will lack power and are not likely to reject stationarity. Further, in all the analysis conducted in the present work, no evidence of non-stationarity in the fMRI time series at hand were evidenced. This is in agreement with what has been reported in He (2011) in an fMRI ROI-based analysis. Finally, previous attempts to scale-free analysis of densely sampled fMRI datasets in time (using the EVI sequence; Rabrait et al., 2008) already confirmed the validity of a the stationarity assumption; see Ciuciu et al. (2008).

7. Conclusion

We uncovered multifractal scale-free dynamics of fMRI time series over four octaves (15–125 s) both in functional networks and in artifacts. We then disentangled functional components from artifactual ones in a robust and significant manner by demonstrating that only the former gave rise to significant modulations of the multifractal attributes between rest and task-related activity. Variability in human performance scores also generally exhibits power-law distributions, whose strength (or exponent) is often modulated across conditions and tasks (Holden et al., 2011). This paves the way toward future works devoted to investigating the extent to which behavioral properties are correlated with the change of scale-free dynamics in neuroimaging time series (MEG, fMRI) acquired during multisensory learning (Seitz et al., 2007).

Conflict of Interest Statement

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

Acknowledgments

The authors thank the French National Research Agency (ANR) for its financial support to the SCHUBERT (ANR-09-JCJC-0071) young researcher project (2009-13). The authors are also grateful to the anonymous reviewers and the associate editor, Dr. BJ He, for their remarks and criticisms that helped us to improve the manuscript and enlarge the readership. Finally, we would like to warmly thank Dr. V. van Wassenhove for her careful rereading and constructive comments.