All fractals are self-similar structures (mathematical fractals in an exact, natural fractals in a statistical sense), with their fractal dimension falling between the Euclidian and topologic dimensions (Mandelbrot, 1983; Eke et al., 2002). When self-similarity is anisotropic, the structure is referred to as self-affine; a feature, which applies to fractal time series (Mandelbrot, 1985; Barabási and Vicsek, 1991; Eke et al., 2002), too. Statistical fractals cannot be described comprehensively by descriptive statistical measures, as mean and variance, because these do depend on the scale of observation in a power law fashion: (1)μ2μ1=s2s1ε, where μ₁, μ₂ are descriptive statistical measures, and s₁, s₂ are scales within the scaling range where self-affinity is present, and ε is the power law scaling exponent. From this definition a universal scale-free measure of fractals can be derived: (2)D=-lims→0inflogNslogs.

D is called capacity dimension (Barnsley, 1988; Liebovitch and Tóth, 1989; Bassingthwaighte et al., 1994), which is related but not identical to the Hausdorff dimension (Hausdorff, 1918; Mandelbrot, 1967), s is scale and N(s) is the minimum number of circles with size s needed to cover the fractal object to quantify its capacity on the embedding dimensional space (it corresponds to μ in Eq. 1). For fractal time series, the power law scaling exponent ε is typically calculated in the time domain as the Hurst exponent (H), or in the frequency domain as the spectral index (β). H and D relate (Bassingthwaighte et al., 1994) as: (3)H=2-D.

Further, β can also be obtained from H as (H − 1)/2 for fGn and (H + 1)/2 for fBm processes (Eke et al., 2000).

While D does not vary along a monofractal time series, it is heterogeneously distributed along the length of a multifractal signal.

This phenomenon gave rise to the term “singular behavior,” as self-affinity can be expressed by differing power law scaling along a multifractal time series, X_i as: (4)Xi+Δi-Xi∝Δih(i), where h is the Hölder exponent defining the degree of singularity at time point, i. Calculating the fractal dimension for each subsets of X_i of the same h, one obtains the singularity spectrum, D(h) (Mandelbrot spectrum), which describes the distribution of singularities (Frisch and Parisi, 1985; Falconer, 1990; Turiel et al., 2006).

(5)D(h)=log(ρ(h)∕ρ(hmax))logsmin, where h_max is the Hölder exponent corresponding to maximal fractal dimension, s_min is the finest scale corresponding to Hölder trajectory, and ρ(h) is the distribution of singularities.

The singular behavior of a multifractal is a local property. Separation of the singularities can be difficult, given the finite sampling frequency of the signal of interest (Mallat, 1999). Thus, in contrast with monofractality, a direct evaluation of multifractality is a demanding task in terms of the amount of data and the computational efforts needed, which can still not guarantee precise results under all circumstance.

With the aid of different moments of appropriate measure, μ, a set of equations can be established to obtain the singularity spectrum, which is a common framework exploited by multifractal analysis methods referred to as multifractal formalism (Frisch and Parisi, 1985; Mandelbrot, 1986; Barabási and Vicsek, 1991; Muzy et al., 1993). Using a set of different moment orders, one can determine the scaling behavior of μ^q, yielding the generalized Hurst exponent, H(q) (Barunik and Kristoufek, 2010; See Figure 1): (6)μq(s)∝sq⋅H(q).

On the right side of Eq. 4 Δi corresponds to scale, s, on the right side of Eq. 6. Using the partition function – introduced in context of Wavelet Transform Modulus Maxima (WTMM) method – singularities are analyzed globally for estimating the (multi)scaling exponent (Mallat, 1999): Z(s,q)=∑k=1N(s)μiq(s)(7)τ(q)=lims→0inflogZ(s,q)logs,(8) where τ(q) can be also expressed from H(q) (Kantelhardt et al., 2002) as: (9)τ(q)=q⋅H(q)-DT, where D_T is the topological dimension, which equals 1 for time series.

The generalized fractal dimension can also describe the scale-free features of a multifractal time series: (10)D(q)=τ(q)q-1=q⋅h(q)-1q-1.

The singularity spectrum, D(h), can be derived from τ(q) with Legendre transform (Figure 2), via taking (11)h=τ′(q), the slope of the tangent line taken at q for τ(q), and yielding (12)D(h)=infq(qh-τ(q)), that when evaluated gives the negative of the intercept at q = 0 for the tangent line (See Figure 2).

Natural signals have a singularity spectrum over a bounded set of Hölder exponents, whose width is defined by [h_−∞, h_+∞] (Figure 3).

A combination parameter, P_c, can be calculated (definitions on Figure 3) to facilitate the separation of time series characteristics (Shimizu et al., 2004), which can aid the exploration of the physiological underpinnings, too.

A similar parameter is W (Wink et al., 2008) calculated as (14)W=W+W-.

Here we focus on widely used monofractal methods selected from those in the literature.

Three analysis methods are described here; all use different statistical moments (termed q-th order) of the selected measure to evaluate the signal’s multifractality. Despite of certain inherent drawbacks, these methods are widely used in the literature, and can obtain reliable results if their use is proper with limitations considered.

Monofractal signals of known autocorrelation (AC) structure can be synthesized based on their power law scaling. The method of Davies and Harte (1987) (DHM for short) produces an exact fGn signal using its special correlation structure, which is a consequence of the power law scaling of the related fBm signal in the time domain (Eq. 19). It is important, that different realizations can be generated with DHM at a given signal length and Hurst exponent, which consists of a statistical distribution of similarly structured and sized monofractals.

The next question is how to define meaningful end-points for the tests? For ideal monofractals with a given length and true H, Mean Square Error (MSE) is a good descriptor: it can be calculated for each set of series of known H and particular signal length, N (Eke et al., 2002). It carries a combined information about bias and variance, as MSE = bias² + variance.

Interpreting the multifaceted results of numerical experiments is a complex task. It can be facilitated if they are plotted in a properly selected set of independent variable with impact shown in intensity-coded representations (Figure 6; Eke et al., 2002). Precision index is determined as the ratio of results falling in the interval of [H_true – H_dev, H_true + H_dev], where H_dev is an arbitrarily chosen value referring to the tolerable degree of deviation.

In the monofractal testing framework, we used DHM-signals to evaluate the performance of MF-DMA (Gu and Zhou, 2010) and MF-DFA (Gu and Zhou, 2006), by the code obtained from http://rce.ecust.edu.cn/index.php/en/research/129-multifractalanalysis. It was implemented in Matlab, in accordance with Eqs 17 and 30–35. As seen in Figure 6, precision of MF-DFA and MF-DMA depends on N, H, and the order of moment.

In order to compare the methods in distinguishing multifractality, end-points should be defined reflecting the narrow or wide distribution of Hölder exponents. We select a valid endpoint Δh proposed by Grech and Pamula (2012), which is the difference of Hölder exponents corresponding to q = −15 and q = + 15 (Figure 7).

Extending the dichotomous model of fGn/fBm signals (introduced in context of monofractals; Mandelbrot and Ness, 1968; Eke et al., 2000) toward multifractal time series is reasonable as it can account for essential features of natural processes exhibiting local power law scaling. Description of an algorithm creating multifractional Brownian motion (mBm) and multifractional Gaussian noise (mGn) can be found here (Hosking, 1984), while implementation of such code can be found on the net (URL1: http://fraclab.saclay.inria.fr/, URL2: www.ntnu.edu/inm/geri/software). Given that these algorithms require Hölder trajectories as inputs, multifractality cannot be defined exactly on a finite set, which is a common problem of such synthesis methods. Selecting a set of meaningful trajectories is a challenging task: it should resemble those of empirical processes and meet the analytical criteria of the selected algorithms (such criteria are mentioned in Concept of Fractal Time Series Analyses).

On the contrary, iterative cascades defined with analytic functions are not influenced by the perplexity of definitions associated with multifractality outlined in the previous paragraph, given that their value at every real point of the theoretical singularity spectrum is known. Due to their simplicity, binomial cascades (Kantelhardt et al., 2002; Makowiec et al., 2012) and Devil’s staircases (Mandelbrot, 1983; Faghfouri and Kinsner, 2005) are common examples of theoretical multifractals used for testing purposes. A major drawback of this approach is that these mathematical objects do not account for features in empirical datasets, but can still be useful in comparing reported results.

The most extensive test of multifractal algorithms which used a testing framework of signals synthesized according to the model introduced by Benzi et al. (1993) was reported by Turiel et al. (2006). Briefly, it is a wavelet-based method for constructing a signal with predefined properties of multifractal structuring with explicit relation to its singularity spectrum. Since the latter can be manipulated, the features of the resulting multifractal signal could be better controlled. The philosophy of this approach is very similar to that of Davies and Harte (1987) in that a family of multifractal signals of identical singularity spectra can be generated by incorporating predefined distributions (log-Poisson or log-Normal) giving rise to controlled variability of realizations. Additionally, using log-Poisson distribution would yield multifractals with a bounded set of Hölder exponents in that being similar to those of empirical multifractals. To conclude, this testing framework should merit further investigation.

Zarahn et al. (1997) demonstrated early in a careful analysis on spatially unsmoothed empirical human fMRI BOLD data (collected under null-hypothesis conditions) that the examined datasets showed a disproportionate power at lower frequencies resembling of 1/f type noise. In spite of the very detailed analysis, these authors treated the 1/f character as a semi-quantitative feature of fMRI noise and accepted its validity over a decaying exponential model as the form of the frequency domain description of the observed intrinsic serial, or autocorrelation. The spectral index, β, however was not reported but can be reconstructed from the power slope by converting the semilog plot of power vs. frequency in their Figure 3D panel to a log-log plot compatible to |A(f ) |²∝1/f ^β model. A β value of ∼3.3 is yielded, which is exceedingly higher than the values of 0.6 < β < 1.2 reported recently for an extensive 3T dataset by He (2011). This precludes the possibility that the collected resting-state 1.5T BOLD dataset would have been of physiological origin. Our recently reported results for the rat brain with −0.5 < β < 1.5 reconfirms this assertion (Herman et al., 2011). In fact, Zarahn et al. (1997) wished to determine if the 1/f component of the noise observed in human subjects was necessarily due to physiological cause, but had to reject this hypothesis because they found no evidencing data to support this hypothesis. Zarahn et al. (1997) felt the AC structure (in the time domain, which is equivalent to the inverse power law relationship in the frequency domain) may not be the same for datasets acquired in different magnets, not to mention the impact of using various fMRI scanning schemes (Zarahn et al., 1997). Accordingly, and in light of our rat data for magnets 4, 9.4, and 11.7T, a less than optimal field strength could have led to a signal definition inadequate to capture the 1/f^β type structuring of the BOLD signal of biological origin that must have been embedded in the human datasets Zarahn et al. (1997) but got overridden by system noise. Most recently, Herman et al. (2011) and He (2011) referred to the early study of Zarahn et al. (1997) as one demonstrating the impact of system noise on fMRI data, while Fox et al. (2007) and Fox and Raichle (2007) as the first demonstration of 1/f type BOLD noise with the implication that the 1/f pattern implied fluctuations of biological origin.

Fox et al. (2007) reported on the impact of intrinsic BOLD fluctuations within cortical systems on inter-trial variability in human behavior (response time). In conjecture of the notion that the variability of human behavior often displays a specific 1/f frequency distribution with greater power at lower frequencies, they remark “This observation is interesting given that spontaneous BOLD fluctuations also show 1/f power spectrum (Figure S4). While the 1/f nature of BOLD fluctuations has been noted previously (Zarahn et al., 1997), we show that the slope is significantly between −0.5 and −1.5 (i.e., 1/f ) and that this is significantly different from the frequency distribution of BOLD fluctuations observed in a water phantom,” and in their Figure S4 conclude that “the slope of the best fit regression line (red) is −0.74, close to the −1 slope characteristic of 1/f signals.” This interpretation of the findings implies that the spontaneous BOLD fluctuations can be adequately described by the “strict” 1/f model, where the spectral index, β, in 1/f ^β – known as the “general” inverse power law model – is treated as a constant of 1, not a variable carrying information on the underlying physiology. Incidentally, studies of Gilden and coworkers (using a non-fMRI approach) have indeed demonstrated (Gilden et al., 1995; Gilden, 2001; Gilden and Hancock, 2007) that response time exhibits variations that could not be modeled by a strict 1/f spectrum but by one incorporating a varying scaling exponent (Gilden, 2009).

Scrutinizing the data of Figure S4 can offer an alternative interpretation as follows. In terms of the hardware, the use of 3T magnet must have ensured adequate signal definition for the study. In their Figure S4, spectral slopes were reported in a lumped manner, in that power at each and every frequency were averaged for the 17 human subjects first (thus creating frequency groups), and then mean slopes along with their statistical variation were plotted for the frequency groups. The mean slope of −0.74 (of thee lumped spectrum) was obtained by regression analysis. This treatment of the data implies that the |A(f )|² ∝ 1/f ^β model (Mandelbrot and Ness, 1968; Eke et al., 2000, 2002) was a priori rejected otherwise the slope should have been determined for each and every subject in the group across the range of observed frequencies and their associated power estimates (of the true spectrum) first, followed by the statistical analysis for the mean and variance within the group of 17 subjects, for the following reasons. The spectral index is found by fitting a linear model of |A(f )|² ∝ 1/f ^β across spectral estimates for a range of frequencies. In our opinion when it comes to provide the mean spectral index, it is indeed reasonable (Gilden and Hancock, 2007; Gilden, 2009) to come up with statistics on the fractal estimates for a group of time series data first by obtaining the estimates, proper. Averaging spectral estimates at any particular frequencies and assembling an average spectrum from them tend to abolish the fractal correlation structure for any particular time series and develop one for which the underlying time series is indeed missing. Because the transformation between the two treatments is not linear, the true mean slope of the scale-free analysis cannot be readily reconstructed from the reported slope of the means. Nevertheless, if we regard its value as an approximation and convert it to β, which being less than 1 warrants the use of H′ = (β_fGn + 1)/2, one would yield a value of β = 0.77 and H′ = 0.87, respectively.

A recent review by Fox and Raichle (2007) offers an impressive overview and insight of how to delineate cooperative areas (or systems) in the brain based on functional connectivity that emerges from spatial cross-correlation maps of regional fluctuating BOLD signals in the resting brain (Biswal et al., 1995). These authors place the spontaneous activity of the brain as captured in BOLD fluctuations in spatio-temporal domains of fMRI data in the focus of the review emphasizing that itis a fingerprint of a newly recognized mode of functional operation of the brain referred to as default or intrinsic mode (Fox and Raichle, 2007). They argue that the ongoing investigation of this novel aspect of the mode of brain’s operation using fractal analysis of resting-state fMRI BOLD may lead to a deeper and better understanding of the way the brain – on the expense of very high baseline energy production and consumption by glucose and oxidative metabolism – maintains a mode capable of selecting and mobilizing these systems in order to respond to a task adequately (Hyder et al., 2006). One has to add that the default or intrinsic mode of operation has been demonstrated and investigated in overwhelming proportions by connectivity analyses based on cross-correlating BOLD voxel-wise signals as opposed on AC of single voxel-wise BOLD time series.

Fox and Raichle (2007) emphasize “spontaneous BOLD follows a 1/f distribution, meaning that there is an increasing power in the low frequencies.” In their furthering on the nature of this 1/f type distribution they refer to the studies of Zarahn et al. (1997) and Fox et al. (2007) in the context it was described above (Fox et al., 2007) reaching the same conclusion, in that the characteristic model of human spontaneous BOLD is the 1/f (meaning the “strict”) model. We would like to suggest that the notion of 1/f distribution having a regression slope of close to −1 on the log-log PSD plot is somewhat misleading.

In an attempt to consolidate this issue, we suggest that the data be fitted to a model in the form of 1/f^β, where β is a variable (Eke et al., 2000, 2002) responding to states of physiology (Thurner et al., 2003; He, 2011) of characteristic topology (Thurner et al., 2003; Herman et al., 2011) in the brain, not a constant of 1. A potential advantage of this model is that by regarding β as a scaling exponent the distribution can then be described to be scale-free (or fractal).

As seen above, from the modeling point of view the issue of a reliable description of the autoregressive signal structuring of spontaneous BOLD, is fundamental and critical in resting-state. If it is done properly, it can lend a solid basis for assessing changes in the scaling properties in response to changing activity of the brain. The study of Thurner et al. (2003) was probably the first to demonstrate that spontaneous BOLD in the brain was scale-free and that the scaling exponent of inactive and active voxels during sensory stimulation differed. At the time of the publication of their study, the monofractal analytical strategy of Eke et al. (2000, 2002) based on the dichotomous fGn/fBm model of Mandelbrot and Ness (1968) did not yet reached the fMRI BOLD community, hence Thurner et al. (2003) did not rely on it, either. In this section we will demonstrate the implications of this circumstance in terms of the validity and conclusions of their study. We will do it in a detailed, didactical manner so that our reader should gain a hands-on experience with the perplexed nature of the issue.

Subtracting the mean from the raw fMRI signal precedes the analysis proper, Īx→(t), yielding Ix→(t) in Eq. 39, (39)Ix→(t)=Īx→(t)-Īx→(t)t, which step is compatible with (D)FA (Eke et al., 2000).

Subsequently, in Eq. 40, the temporal correlation function, Cx→(τ), is calculated (40)Cx→(τ)=Ix→(t)Ix→(t+τ)=1N-τ∑t=1N-τIx→(t)Ix→(t+τ).

In fact in this step of the analysis the covariance was calculated given that a division by variance was missing. Hence, it is slightly misleading to regard Eq. 40 as the temporal (or auto) correlation (see Eke et al., 2000, Eq. 2). Only, if assumed that the signal is fGn, whose variance is known to be constant over time, the covariance function can be taken as equivalent to the AC function. Because the authors have not tested and proven the signal’s class was indeed fGn (Eke et al., 2000), there is no basis for the validity of this assertion.

In Eq. 41, the signal is summed yielding Xnx→(τ), in order to eliminate problems in calculating the AC function due to noise, non-stationarity trends, etc.

This form of the signal is further referred to as “voxel-profile.”

Note, that the signal remains in this summed form for the rest of the analysis (i.e., analyzed as fBm). As a consequence, spectral analysis later in the study was applied to a summed – hence processed – signal and the results were thus reported for this and not the raw fMRI signal, which circumstance prevented reaching a clear conclusion.

Furthermore, the authors indicated that the temporal correlation function would characterize persistence. It seems the two terms (correlation vs. persistence) are used as synonyms of one another whereas they are not interchangeable terms: persistence is a property of fBm, while correlation is that of an fGn signal (Eke et al., 2000). Please note, as the raw signal has been summed, the covariance here characterized persistence that was not present in the raw fMRI signal.

In the next step (Eq. 42), the AC function is approximated by a power law function with γ as its exponent (42)Cx→(τ)~τ-γ,⋯0<γ<1.

Based on the equation of the AC function using the Hurst exponent, H, γ must be proportional to 2H (Eke et al., 2000, Eq. 15).

Subsequently, as a part of a FA of the authors (cited in their Reference 19 as unpublished results of their own), the statistics (Fx→(τ), standard deviation) was calculated for the AC function in Eq. 43 (43)Fx→(τ)=Xn+τx→-Xnx→2n1/2.

In the left side of Eq. 44, a general power law was applied to the fluctuation from Eq. 43 as Fx→(τ)~τα (44)Fx→(τ)~τα,α=1-γ/2.

(Note, as the fluctuations have not been detrended, this method is not the DFA of Peng et al., 1994 but strongly related to it).

Consider the scaling exponent, α, on the left side of Eq. 44. According to Peng et al. (1994) and Eke et al. (2002) α = H only if the raw signa l,Ix→(t), is an fGn. However, because at this point the summed raw signal, Xnx→(τ), is the object of the analysis, α and H should relate to each other as α = H + 1. Given that the signal was summed in Eq. 41 leading up to Eq. 43, and the values for “outside the brain” were reported as α ≈ 0.5, and for “inside the brain” as 0.5 < α < 1, α must have been improperly calculated because α cannot possibly yield a value of 0.5 for a summed signal given that H scales between 0 and 1 and for an fBm series α = H + 1 holds. The reported value of 0.5 < α < 1 can be regarded correct only for Ix→(t), the raw fMRI signal, which therefore had to be an fGn process. On the other hand, the reported values 2 < β < 3 are correct for the Xnx→ signal, only (for reasons given later). Hence the reported α and β values lacking an indication of their respective signal class ended up being ambiguous.

Next, consider the right side of Eq. 44, which expresses α by using γ introduced earlier. We just pointed out that the raw fMRI signal must have been an fGn with α ≡ H. Consequently, α can be substituted for H in Eq. 44 as H = 1 − γ/2 and γ expressed as (45)γ=2-2H.

The authors referring to power law decays in the correlations relate the spectral index, β, to γ as (46)β=3-γ, and further to α as β=2α+1.

Note, that these relations between β, γ, and α in principle do depend on signal class that was not reported.

Now, let us substitute γ as expressed in Eq. 45 into Eq. 46 β=3-2+2H=1+2H, then express H (47)H=β-12.

As shown by Eke et al. (2000), Figure 2; in Eke et al. (2002), Table 1, based on the dichotomous fGn/fBm model, Eq. 47 would have equivocally identified the case of an fBm signal. As pointed out earlier, the raw fMRI signal was summed before the actual fractal analysis. Consequently, the relationship β = 3 – γ ends up holding only if the raw fMRI signal was an fGn process. This is therefore the second piece of evidence suggesting that the class of the raw fMRI signal must have been fGn. Nevertheless, the relationship β = 2α + 1 could not hold concomitantly for reasons that follow. In an earlier paper of the group (Thurner et al., 2003), the authors stated “The relationship is ambiguous, however, since some authors use the formula α = 2H + 1 for all values of α, while others use α = 2H−1 for α < 1 to restrict H to range (0,1). In this paper, we avoid this confusion by considering α directly instead of H.” The fGn/fBm model (Eke et al., 2002) helps resolving this issue as neither of these relationships between α and H holds because if α is calculated with the signal class recognized and determined, the relationship between α and H is equivocally α = H_fGn and α = H_fBm + 1. Based on the fGn/fBm model, the relationship between β and α given in Eq. 46 as β = 2α + 1 needs to be revised, too, to is correct form of β = 2α − 1 (See Table 1 in Eke et al., 2002).

Thurner et al. (2003) concluded: “Outside the brain and in non-active brain regions voxel-profile activity is well described by classical Brownian motion (random walk model, α ∼ 0.5 and β ∼ 2).” Recall, the “voxel-profile” is not the raw fMRI signal (intensity signal, Ix→(t), most probably an fGn), but its summed form, Xnx→(τ), an fBm.

Our conclusion on the above analysis by Thurner et al. (2003) is as follows: (i) α was improperly calculated by the authors’ FA method because α ∼ 0.5 cannot possibly be valid for an fBm signal given that α_fBm > 1 (Peng et al., 1994), (ii) β ∼ 2 is only formally valid given that it was calculated based on Eq. 46 from an improperly calculated α and by using an arbitrary relationship between α and β. The subsequent and opposite effects of these rendered the value of β to β ∼ 2.

When the results of Thurner et al. (2003) are interpreted according to the analytical strategy of Eke et al. (2000) based on the dichotomous fGn/fBm model of Mandelbrot and Ness (1968), the reported values of Thurner et al. can be converted for their fMRI “voxel-profile” data Xnx→ to α_fBm ∼ 1.5, β_fBm ∼ 2, H_fBm ∼ 0.5 or for the raw fMRI intensity signal Ix→(t) to α_fGn ∼ 0.5, β_fGn ∼ 0, H_fGn ∼ 0.5. This interpretation of the data reported for humans by Thurner et al. (2003) is fully compatible with the current findings by He (2011) on the human and by Herman et al. (2009, 2011) on the rat brain.

Exemplary analysis on empirical BOLD data is presented on the 11.7T coronal scan shown in Figure 10 to demonstrate the inner workings of these methods when applied to empirical data, and point to potential artifacts, too (See Figure 12). For monofractal analysis, we recommend using monofractal SSC for it gives unbiased estimates across the full range of the fGn/fBm dichotomy. For this reason, the topology is well defined and not as noisy as on the PSD maps. MF-DFA, due to its inferior performance in the strongly correlated fGn range (See Figure 6 at q = 2), failed with this particular BOLD dataset. Also note, that the histograms obtained for the same datasets evaluated by these different methods do differ indicating that method’s performance were different. Proper interpretation of the data therefore assumes an in-depth understanding of the implication of method’s performance on the analysis. P_c and most certainly W seems a promising parameter to map from the BOLD temporal datasets. Their proper statistical analyses along with those of singularity spectra for different anatomical locations in the brain should be a direction of future research.