METHODS article

Front. Physiol., 26 July 2017

Sec. Fractal Physiology

Volume 8 - 2017 | https://doi.org/10.3389/fphys.2017.00533

Decomposing Multifractal Crossovers

  • 1. Institute of Clinical Experimental Research, Semmelweis University Budapest, Hungary

  • 2. Department of Physiology, Semmelweis University Budapest, Hungary

  • 3. Department of Radiology and Biomedical Imaging, Yale University New Haven, CT, United States

Abstract

Physiological processes—such as, the brain's resting-state electrical activity or hemodynamic fluctuations—exhibit scale-free temporal structuring. However, impacts common in biological systems such as, noise, multiple signal generators, or filtering by transport function, result in multimodal scaling that cannot be reliably assessed by standard analytical tools that assume unimodal scaling. Here, we present two methods to identify breakpoints or crossovers in multimodal multifractal scaling functions. These methods incorporate the robust iterative fitting approach of the focus-based multifractal formalism (FMF). The first approach (moment-wise scaling range adaptivity) allows for a breakpoint-based adaptive treatment that analyzes segregated scale-invariant ranges. The second method (scaling function decomposition method, SFD) is a crossover-based design aimed at decomposing signal constituents from multimodal scaling functions resulting from signal addition or co-sampling, such as, contamination by uncorrelated fractals. We demonstrated that these methods could handle multimodal, mono- or multifractal, and exact or empirical signals alike. Their precision was numerically characterized on ideal signals, and a robust performance was demonstrated on exemplary empirical signals capturing resting-state brain dynamics by near infrared spectroscopy (NIRS), electroencephalography (EEG), and blood oxygen level-dependent functional magnetic resonance imaging (fMRI-BOLD). The NIRS and fMRI-BOLD low-frequency fluctuations were dominated by a multifractal component over an underlying biologically relevant random noise, thus forming a bimodal signal. The crossover between the EEG signal components was found at the boundary between the δ and θ bands, suggesting an independent generator for the multifractal δ rhythm. The robust implementation of the SFD method should be regarded as essential in the seamless processing of large volumes of bimodal fMRI-BOLD imaging data for the topology of multifractal metrics free of the masking effect of the underlying random noise.

Introduction

Fractal and multifractal concepts focus on characterizing scale-free properties in terms of scaling exponents—such as, spectral index (β) or Hurst exponent (H; Mandelbrot, ; Eke et al., , ; Mukli et al., )—of ideal or empirical signals. The scaling is a global behavior in the case of monofractals and a local property in the case of multifractals, which requires a set of exponents to be obtained for characterization. Specifically—in addition to a range of methods operating in the frequency and time/frequency domains (Eke et al., )—in the time domain, this is achieved by analyzing a range of statistical moments (−∞ < q < + ∞) of the signal. In the monofractal case, a single qth order moment (i.e., the variance at q = 2) suffices for capturing the global roughness, H. However, for multifractals, a range of statistical moment orders are needed to obtain the generalized Hurst exponent, H(q). Submitting H(q) to the multifractal formalism yields the Hölder exponent, h, reflecting the local roughness of the process, and then the multifractal spectrum, D(h), which is essentially analogous with a histogram of local fractality in the signal. Accordingly, D(h) captures the moment-wise distribution of the singularity strength of local roughness or multifractal scaling in the temporal process (Kantelhardt et al., ; Ihlen, ; Mukli et al., ). We recently demonstrated that standard moment-based multifractal analyses were susceptible to signal inhomogeneity leading to spurious estimates of the multifractal spectrum. We resolved this issue by developing focus-based multifractal formalism (FMF), which replaced the standard—essentially monofractal—analysis for H(q) by fitting an exact multifractal to the family of moment-wise scaling functions all at once by enforcing an expected value at signal length (termed focus) as a guiding reference in the fitting procedure (Mukli et al., ). FMF explicitly relied on a previous observation on the focus (Kantelhardt et al., ) and can be related to some earlier multifractal approaches (Struzik, ; Struzik and Siebes, ).

In the pure mathematical sense of the fractal concept, scaling should be present across an infinite range of scales; a property of ideal fractals with an exact generating algorithm (Mandelbrot, ), such as, Cantor set and function (Cantor, ). Fractality can be present in a statistical sense in sampled representations of temporal processes, as it is the case with fractional Gaussian noise (fGn) and Brownian motion (fBm; Mandelbrot and Van Ness, ; Eke et al., , ). However, the estimation of fractality of even such exact fractal structures can become easily corrupted by the effect of sampling (see Figure 1), filters (Valencia et al., ), trends (Kantelhardt et al., ), shuffling (Kantelhardt et al., ), multiple fractal signal components (Thornton and Gilden, ), or other scale-dependent influences, resulting in multimodal scaling functions.

Figure 1

) within a temporal window with a lower boundary, s′, and an upper boundary set to N/2. The estimate of the generalized Hurst exponent at q = 2 [H(2)] was obtained on multifractal signals by the FMF-SSC method (Mukli et al., ) within the same range of scales. DHM fractal: synthesized by the method of Davies and Harte (DHM; Davies and Harte, ) at Htrue = 1.25; SSM fractal: synthesized by the SSM method (Saupe, ); Cantor set and function generated at appropriate weight factors p1, p2, p3 (Cantor, ; Mandelbrot, ). Estimated Hurst and generalized Hurst exponents are displayed on an extended scale within the range of [0, 2]. Signal classes based on the extended H and H(2) are shown to the right for reference. Contracting scaling range—by increasing s′—will render the information remaining in the signal for its analysis inadequate to demonstrate its scale-invariance (see the bias in reference of Htrue progressively increasing with s′). A standard treatment of scaling function multimodality typically results in contracted scaling ranges and thus results in a loss of valuable information on fractal scaling.

Many physical, natural, biological systems show multimodal, scale-invariant properties, for example, sunspot activities (Movahed et al., ), river water levels (Rego et al., ), human heartbeat time series (Peng et al., ; Gierałtowski et al., ), neuronal discharge dynamics (Blesic et al., ), human near infrared spectroscopy (NIRS) signals (Eke et al., ), local field potential (Bedard et al., ), human electroencephalography (EEG) signals (Gifani et al., ), rat blood oxygen level-dependent functional magnetic resonance imaging (fMRI-BOLD) data (Herman et al., ), human gait trajectories (Kuznetsov et al., ), and neonate diffuse optical tomography data (White et al., ). So far—among the numerous possible scale-dependent influences—only trends leading to artificial crossovers (Kantelhardt et al., ) have been specifically assessed within such multimodal signals. Common among these examples is that multimodality emerges, e.g., by superposition of two or more coexisting multifractal processes. Thermodynamic analogies of multifractality (Stanley and Meakin, ; Tel, ; Arneodo et al., ) have suggested that superpositioning could also lead to phase transition (Grassberger et al., ; Muzy et al., ; Arneodo et al., ; Radons and Stoop, ). Multimodality of the former case is seen in the initial step of multifractal formalism resulting in a scale-dependent impact on the scale-free pattern, while phase transition manifests itself at a later stage as a moment-dependent phenomenon captured as an inflected multiscale exponent, τ(q), at a critical q. This latter case has been detailed in the literature due to its demonstrated physiological relevance (Roux et al., ; Nicolay et al., ; Kestener et al., ); thus, our study focused on the former case of multimodality only.

As characterization of multifractality in the time domain requires assessing moment-wise scaling exponents (Kantelhardt et al., ), multifractal multimodality should likewise also be described—along with the respective scaling exponents—in terms of moment-wise breakpoints or crossovers (Ludescher et al., ; Schumann and Kantelhardt, ; Gierałtowski et al., ; Figures 2B,C). In contrast, a mere presence of breakpoints or crossovers should by no means be taken as proof of multifractality (Matic et al., ).

Figure 2

). (B) Our concept of q-wise scaling range adaptive focus-based multifractal formalism. This approach employs iterative fitting for q-dependent crossover scales by enforcing the focus (black circle) of the scaling functions found at maximum signal length. Note that the point-like q-wise breakpoints separate two adjacent SRs occupied by distinct scale-invariant components of the multifractal scaling functions. (C) Scaling function decomposition with focus-based multifractal formalism utilizing an extended version of the iterative approach shown in (B). It yields a complete decomposition of the overlapping scaling functions of the merging fractals/multifractals coexisting within the same SR.

Definition of “breakpoints” or “crossovers” appears inconsistent in the literature (Peng et al., ; Iyengar et al., ; Struzik et al., ; Kantelhardt et al., ; Ge and Leung, ; Kuznetsov et al., ; Rego et al., ), and their extension to multifractal multimodality (Figure 2A) is also lacking. In particular, recent works have either focused mostly on monofractal multimodality (Kuznetsov et al., ) or applied a monofractal approach to multimodal handling of multifractal signals (Ge and Leung, ; Rego et al., ; Figure 2A).

Previously, “breakpoints” or “crossovers” were determined by “eyeballing” or by segmented line regression (Ge and Leung, ) through scaling ranges (SRs). The latter being a typical adaptive analysis that identifies adjacent SRs with different scaling separated by a point-like breakpoint. As breakpoint is not a point-like attribute of the scaling function, some approaches omitted the transient in its vicinity (Kuznetsov et al., ) or applied fitting with multi-parametric models (Struzik et al., ). Admittedly, the benefit of these approaches is that they allow an approximation of fractal descriptors without formulating any a priori concept on signal genesis. Nevertheless, empirical signals much too often result from co-sampled signal processes, whose scaling may overlap across their respective breakpoints (Figure 2B); one such example is signal contamination by instrumental or biological noise. In these cases, the signal components obviously cross over around the apparent breakpoint (Figure 2C). This perception of multimodal signal genesis calls for an adequate concept for handling crossovers based on a genuine multifractal approach.

Accordingly, our aims were (i) to decompose the moment-wise crossover of superimposed multifractal signals based on an additive model, (ii) to validate this method, (iii) to compare this approach with an enhanced—moment-wise—version of the segmented line regression method, and (iv) to demonstrate their applicability on exemplary empirical signals.

Materials and methods

Description of the methods of signal synthesis and empirical data acquisition are followed by introduction of two adaptive multifractal analyses of bi- or multimodal signals. The first approach is based on a q-wise identification of breakpoints along each and every scaling function as a step of signal pre-treatment, and hence is referred to as the q-wise scaling range adaptive (qSRA) method. The essence of the second method is to decompose the multifractal crossovers for all scaling functions of the analysis combined, thus achieving true scaling function decomposition (SFD) in a one-pass manner. As both apply to the regression scheme of our FMF (Mukli et al., ), these will be further referred to as the qSRA-FMF and SFD-FMF methods.

The multifractal algorithms, signal synthesis, and numerical tests were implemented in Matlab (The MathWorks, Inc., Natick, MA, USA) with code written by the authors. The multifractal toolbox containing scripts described in this paper can be requested from the corresponding author.

Investigated signal populations

Synthesized monofractal time series

As described previously (Eke et al., ), statistically self-affine time series of the fGn type with exact long-range-dependent structuring were generated by the Davies and Harte method (DHM; Davies and Harte, ) at length, N, and pre-set Hurst exponent, fGnHtrue. fBm-type time series were produced by cumulatively summing (Eke et al., ) their fGn counterparts, yielding a true extended Hurst exponent of Htrue = fBmHtrue = fGnHtrue + 1. fBm-type time series were also produced directly at N and pre-set Htrue = fBmHtrue by the spectral synthesis method (SSM; Saupe, ).

Synthesized multifractal time series

Cantor sets and—by their cumulative summation—Cantor functions as examples of exact multifractal structuring were generated at pre-set weight factors. Statistical self-similar multifractals with known H(q) were synthesized for testing purposes using the generalized binomial multifractal model (Oświȩcimka et al., ; Schumann and Kantelhardt, ). The generating algorithm is an iterative process in which a stochastic binomial cascade is created at a given weight factor and at a given length. In the generalized method, the degree of correlation and the strength of multifractality can be separately tuned, the former being set with the Fourier Filtering Method; for further details see Schumann and Kantelhardt () and Mukli et al. ().

Synthesized multimodal time series

After considering various numerical testing frameworks, we chose the DHM algorithm and Cantor functions as offering the best control over the cardinal parameters in our testing. The above-listed signals represent cases of unimodality with a single SR. Multimodal synthetic (or mock) signals with multiple SRs were created by adding these unimodal fractal time series of known attributes (N, Htrue). Positioning of crossovers was controlled by setting the standard deviation (SD) ratios (i.e., focus ratio) of the signal components in addition to the differences in Htrue.

Sampled empirical time series

Human NIRS measurements using a NIRO 500 Cerebral Oxygen Monitor (Hamamatsu Photonics, Hersching, Germany) at a rate of 2 Hz were carried out to record the relative change in total hemoglobin concentration with a length of N = 16,384 data points (for details, see Eke et al., ). Human EEG signals were sampled with a length of N = 16,384 data points with eyes closed during random hand movements at 500 Hz using a Neurofax EEG System (see Acknowledgment; Pattnaik and Sarraf, ). These empirical records were acquired from healthy volunteers above the pre-frontal area. Rat fMRI-BOLD data with a length of N = 4,096 data points were obtained by using a modified 11.7 T Bruker horizontal-bore spectrometer (Bruker AVANCE, Billerica, MA, USA) using a 1 H surface coil (1.4 cm diameter) with sequential sampling gradient echo planar imaging (EPI) sequence (Hyder et al., ): field of view of 2.56 × 2.56 cm2; image matrix of 64 × 64; slice thickness of 2 mm; repetition time of 200 ms (i.e., 5 Hz of sampling frequency) and echo time of 13 ms; and voxel size of 400 × 400 × 2000 μm3 (for further details see Herman et al., ).

Multifractal analyses

According to its indirect concept, multifractal characterization of time series is performed by sequencing through the steps of scaling, regression, and singularity analyses of the multifractal formalism (Mukli et al., ). The output of the first step is the scaling function (S) obtained in the time domain [i.e., by DFA (detrended fluctuation analysis; Peng et al., ), SSC (signal summation conversion; Eke et al., ), or DMA (detrending moving average; Gu and Zhou, )] and in the time-frequency domain (i.e., by WTMM (wavelet transform modulus maxima; Muzy et al., ) and WL (wavelet leader-based multifractal analysis methods; Jaffard, ). S is calculated from the input signal (Xi) as a function of scale (s) and moment order (q) by obtaining the power mean of the measure (μ) with q as the exponent.

where Ns stands for the number of non-overlapping windows, and v for different temporal positions within a particular signal segment of size s. For further details, see Kantelhardt et al. () and Ihlen ().

Levels of moment order were selected from −15 to 15 in increments of 1, based on (i) the findings of Grech and Pamuła () and (ii) the scaling parameter (Clauset et al., ) obtained from the thin-tailed distribution of our synthetic and empirical data. The scaling function of the signal was sampled in 100 logarithmically spaced scales, which is a much denser representation than dyadic sampling would be for a signal length N = 214. Low and high temporal scales were excluded from our analyses (Cannon et al., ; Kantelhardt et al., ), resulting in a SR of 8 ≦ s ≦ N/2. The applied data were also tested by DFA with high-order detrending, effectively excluding an “artificial crossover” (Kantelhardt et al., ). We did not use filtering, as it produces a breakpoint effectively splitting the SR (Valencia et al., ). Disregarding this scenario (Matic et al., ) may well end up in the analysis of a filtering-ridden band for multifractality instead of the range of scale-invariance of the physiological signal.

Focus-based multifractal method

The FMF (Mukli et al., ) exploits the fact that the family of q-wise scaling functions yields a convergent structure merging at the coarsest scale, s = N. At this scale, the value defined in Equation (1) becomes q-independent yielding a “focus,” S[Xi](N). Enforcing S[Xi](N) as an analytical constraint improves the extraction of q-wise regression slopes for finite-length (by excluding inherently conflicting estimates due to the merely explicit emergence of this restriction) and thus offers a robust method for multifractal analysis of empirical data (Mukli et al., ; Ali et al., ; Delignières et al., ), the very approach adopted in this work. Instead of independently performing repetitive monofractal analyses on a set of empirical scaling functions (Equation 2), FMF methods—through minimizing the residual sum of squared errors, SSE—iteratively find the best-fitting, true multifractal and assign it to the family of the evaluated scaling functions (Equation 3).

In Equation (2), x has two specific values (at x = 0 and in the case of the “focus” x = N); otherwise, it represents the scale, where the exact scale-dependent statistic is being evaluated. The case of x = N and represents the enforced constraint in Equation (3). Thus, according to FMF, a set of model (i.e., exact) scaling functions with iterated parameters—Ĥ[Xi](q) and log Ŝ[Xi](N)—are fitted all at once to the actual data set of the scaling functions. In order to obtain an overall measure of the goodness-of-fit of the FMF regression procedure, its mean squared error (MSE) was calculated according to Equation (21) of Mukli et al. ().

Moment-wise scaling range adaptivity method

The standard segmented line regression method is capable of finding breakpoints, sb, and also in the case of a superimposed signal (Figure 3) approximating crossovers, sx. Equation (4) is an adaptation of the segmented line regression method for a bimodal scaling function, where s′ could be any particular temporal scale. To capture q-dependent breakpoints, we introduced a q-wise regression algorithm, broken down into three steps of Equations (4a–c)

where indices f and n stand for different fractal processes: in our particular case uncorrelated (noise) and correlated (fractal) signals within a co-sampled arrangement, respectively. We chose noise and fractal signals as the constituents of a bimodal signal in describing our methods because this was the case for bimodal cerebral hemodynamic data reported earlier (Eke et al., ; Herman et al., ) and is used in this study as an exemplary dataset. The breakpoint for a given moment is obtained at the minimum value of the SSE(q,s′) function as the estimates, Ĥ[Xi](q) and log Ŝ[Xi](q,x) for the fractal and noise components, respectively, are being refined during the iteration process. Further away from the breakpoint, in the low range of scales, the underlying noise with low H will dominate, while in the high range of scales, the fractal component with high H will dominate. However, in between the two, where the respective function values are commensurable, their fluctuations blend into a non-fractal segment, a so-called exclusion range (ER; Figure 3). Consequently, by excluding the ER from the analysis, the error in assessing the low and high Hs will be smaller, while for the FMF analysis yielding H(q), the error will increase due to SR being contracted by the exclusion range. To find the ER, we used Equation (5) (a similar approach as Equation (20) of Ge and Leung, ) at a specific q and selected a range of scales with SSE(q,s′) lower than SSEtolerance calculated at a preset level of tolerance (0 ≤ tolerance ≤ 1)

Any moment-to-moment inconsistencies in the regression analysis will upset the expected structural aspect of multifractal scaling functions known as the “H(q)-monotonicity” (Mukli et al., ), the monotonous drop of regression slopes of H(q) from q = −∞ to q = +∞, because in subsequent steps of multifractal formalism this automatically results in “inversed” or “corrupted” multifractal spectra. To eliminate this eventuality, the analytical constraints of H[Xi](q) > H[Xi](q-k) with k > 0 were enforced on the nested iterative processes of minimization. MSE was obtained from the sum of SSEs from respective FMF analyses see Equation (3) on the noise and fractal components separated by q-wise breakpoints.

Figure 3

Scaling function decomposition method

As fluctuations from the two underlying signal components mutually contribute to each other's scaling functions near the breakpoint, they hold estimates deviating from the power-law relationship (Figure 4). Taking the exemplary case of the SSC algorithm—where the statistical measure is the standard deviation, SD—this relationship is readily seen as a realization of the Bienaymé formula stating that in the case of uncorrelated variables, the variance (SD2) of their sum equals the sum of the respective variances (Bienaymé, ). A generalization for the case of (anti)correlated signals is given in the Appendix (See Supplementary Material). Thus, the emerging scaling function (the power mean of SD) can be well estimated (Figure A1 in Supplementary Material) as the root sum square of the composing scaling functions as

where the signal (cXi) used in the calculation of S is in square brackets with c being a positive integer referring to each and every of the Nc constituent signals.

Figure 4

Earlier—for the cases of resting-state cerebral hemodynamic fluctuations—we showed that a fractally correlated signal is typically interwoven by uncorrelated noise (Eke et al., ; Herman et al., ). In this simplest case of bimodality, the above derivation can be reduced to the following direct approximation

Accordingly, instead of fitting the two constituting fractally correlated components of a bimodal scaling function separately in two distinct processes, an exact bimodal model scaling function is reconstructed from two properly fitting power-law sets, based on the rule of addition Equation (7). Performing this one-pass regression on a log-log scale, the minimization of SSE—with the generalized Hurst exponent and the focus of the scaling function being iterated—results in the best fit of the exact bimodal model as follows

A special application of this procedure is when one of the constituting components of the composite signal is known, which obviously reduces the number of tuning factors in the minimization process. Specifically, this component could be uncorrelated noise [i.e., instrument and/or biological noise (Peng et al., ; Blesic et al., ; Eke et al., ; Herman et al., )] with a fractal estimate of Ĥ ≅ Htrue = 0.5 [as will be seen later in Section Performance of the SFD Method on an Empirical Bimodal Signal with Limited Definition (fMRI-BOLD)].

SFD is not at all limited to q-wise applications, but can also be performed along with FMF. In this case, the process of minimization of the FMF analysis needs to be modified by raising the number of tuning parameters Equation (9). Thus, both the two sets (nXi and fXi) of H(q) and their associated two foci, S(N), see Equations (9a and 9b) are being simultaneously adjusted in the same iterative process Equation (9c)

Similarly to the qSRA method, H(q)-monotonicity was granted by applying the same analytical constraints. The calculation of MSE from SSE Equation (9) was as explained in Section Focus-Based Multifractal Method.

The crossover scale, sx, of the decomposed scaling functions—where the respective statistical values are in principle the same—can be determined as the common value of the equations of the two underlying regression lines

Thus, the crossover scale can be calculated as

When enforcing the respective foci of the underlying components, the best value of the crossover scale is obtained as

Characterization of methods

To assess the precision of our novel approaches in analyzing multifractal bimodal signals, estimates were compared with multifractal endpoints derived from the singularity spectrum, D(h), (Figure 5) and the results presented in the form of performance vignettes (Eke et al., ). These endpoints are hmax [the Hölder exponent, h, at the peak of D(h)] and fwhm [the full width of D(h) at half maximum] (Eke et al., ). hmax captures the degree of correlation, and fwhm can be taken as the measure of the supporting base of the singularity spectrum (see Figure 3 of Mukli et al., ).

Figure 5

For obtaining references for the estimates by subsequent SFD-FMF and qSRA-FMF methods, synthetic signal components were analyzed for their respective multifractal estimates with the FMF-DFA and FMF-SSC methods ((Mukli et al., ); Figure 5). The characterization of our methods was completed by testing their performance on empirical NIRS, EEG, and fMRI-BOLD signals.

Results

Impact of moment level on crossover scales

As seen in Figures 5A1,A2, the crossover between two components of markedly different correlation structuring is easy to detect. When H approaches Hnoise—as the true breakpoint becomes poorly defined—the bimodal signal approaches unimodal. A similar scenario is seen with the impact of moment level (Figures 5B1,B2), where the actual scale-wise distribution of crossovers will be determined by the dynamics of the H(q) of the signal components.

Performance of qSRA and SFD methods on synthetic bimodal signals

In addition to the impact of correlation and moment level, the focus has a decisive impact on how markedly a signal component dominates the bi- or multimodality of a composite signal (Figure 6). Accordingly, depending on the actual signal component, H and the component focus ratio (or SD ratio) together will impact the direction and magnitude of bias in the multifractal estimates (hmax and fwhm) for the two approaches alike. When the aim is to provide a characterization of multifractality for a bi- or multimodal multifractal signal (i.e., with hmax and fwhm, combined), the actual combination of H and the focus ratio should preferably be as close as possible to the diagonal band of low bias (Figure 6, combined).

Figure 6

). The mean of the absolute biases for a combined evaluation of the fractal and noise components was also created and displayed in contour plots with the actual coordinates of the two empirical signals (EEG, NIRS) overlayed.

An accurate multifractal output is at most partially qualified to assess performance. Determination of the proper method for a given signal is also a requirement. Lower MSE levels in the estimates of the SFD method when compared with those obtained by qSRA analysis suggests that the signals emerged as sums of two underlying scale-free processes, in which case the SFD method should be preferred. The performance of the SFD analysis was tested on the synthesized data pool used in Figure 6. The crossover-model, eventually identified by comparing the MSEs of our two methods (qSRA and SFD), showed a sensitivity of 73%.

Performance of qSRA and SFD methods on high-definition empirical bimodal signals (EEG and NIRS)

High-definition empirical signals (EEG and NIRS in Figure 7) were chosen for demonstrating the optimal performance of the qSRA and SFD methods on empirical data. Both of these data sets had a combination of H and focus ratio close to the low-bias band of these methods (as seen in Figure 6, combined). The SFD method proved superior on these signals over the qSRA approach, yielding lower MSE-values and values of a magnitude lower when compared with those of the unimodal analysis (Table 1). Synthesizing the signal components based on the endpoint parameters of the SFD-FMF analysis (Table 2) yielded the same MSE when these components were added. This supports the notion that these bimodal signals could be treated as the sum of two concomitant processes, of which one could be fitted by an exact multifractal (Mukli et al., ) and the other by an exact multifractal noise (Grech and Pamuła, ). The multifractal analysis by the SFD FMF-SSC method for H(2) and hmax yielded a correlated level for the fractal and an uncorrelated level for the noise component. The multifractal measure (fwhm) was greater for the former than for the latter as shown in Table 2.

Figure 7

Table 1

Method/SignalEEGmock EEGNIRSmock NIRS
FMF-SSC0.52260.53860.13880.2776
SRqA FMF-SSC0.03530.04090.03880.0550
SFD FMF-SSC0.03110.03200.02370.0260

The goodness-of-fit statistics (MSE) of the raw (FMF-SSC) and the two adaptive FMF-SSC methods (qSRA and SFD) for the empirical signals and their numerical equivalents shown in Figure 7.

Table 2

MethodEndpointEEGNIRS
NoiseFractalNoiseFractal
SFD FMF-SSChmax0.411.830.551.26
fwhm0.180.660.240.54

The endpoint parameters of SFD-FMF-SSC analysis of exemplary bimodal empirical signals shown in Figure 7.

The scaling functions for the EEG and NIRS data sets are shown in Figure 8. The estimated crossover scale of the human EEG is 257 ms at q = 2 and in case of NIRS records is 46 s at q = 2 (Figure 8). This demonstrates that the identified moment-wise crossover scales correspond well with characteristic boundaries between the theta and delta bands of the EEG and, in the case of NIRS signals, the transient is in-between the low- (Biswal et al., ; Herman et al., ) and high-frequency fluctuations (Figure 8). In this particular case, our analyses confirmed that bimodality in the analyzed EEG and NIRS signals (recorded from the human brain cortex) should not be regarded as a scale-dependent imprint of a transfer function but as one resulting from the superposition of random noise and correlated multifractal processes (Figures 7, 8 and Tables 1, 2).

Figure 8

Performance of the SFD method on an empirical bimodal signal with limited definition (fMRI-BOLD)

Rodent fMRI-BOLD imaging data of limited definition (Eke et al., ) was chosen to demonstrate the performance of our SFD method on scans of bimodal BOLD time series data obtained across a section at the level of the bregma (Herman et al., ; Figure 9A). Criteria of 0.5 < H[fXi](2) < 2, sx within the full available SR, and 0.1 < sb/sx < 10 yielded valid estimates of crossovers with MSE:SFD-FMF < MSE:qSRA-FMF at a rate of 68%. Accordingly, this assessment led to frequent dropouts in the parametric images at locations with extreme sx (data not shown). We could, however, make the performance of the SFD method robust on limited-resolution BOLD data by exploiting the fact that a ubiquitous uncorrelated component with minor variations was present in the scaling functions throughout the section (Figure 9A) and thus could be used as a default in the analysis. Building on the maximal precision of the SFD method as demonstrated on synthetic signals (Figure 6), this default component was determined in a voxel with a mid-scale sx (Figures 9A,C, marked in red). A reliable assessment of sx(q) resulted in consistent estimates of H(q) and fwhm as seen in the parametric maps in Figure 9D. Note that, in our case, because of the impact of the uncorrelated component, lower crossover scales appeared to be associated with lower Hurst exponent values in cortical areas with the exception of the cingulate cortex, where the correlated component was very dominant. This yielded apparently inverse patterns in the H[fXi](2) and H[Xi](2) maps (Figures 9B,D). Furthermore, the crossover scale seen in the human NIRS signal (Figure 8) fell within the range of crossover scales found in the rodent fMRI-BOLD image data (Figure 9D).

Figure 9

) by the SFD-FMF method are shown along with those obtained by the FMF method for comparison. (A): Representative scaling functions for cortical (C), subcortical (S), and basal (B) voxels are shown. The uncorrelated component determined from a voxel with mid-range sx [marked in red in (A,C)] that fits very well as a ubiquitous noise component for all scaling functions. Standard, unimodal FMF-based analysis yielded the multifractal metric maps seen in (B). The concept of the SFD-FMF-based adaptive, bimodal analysis is shown in (C). The ubiquitous uncorrelated component (red) seen in (A) was used as default in decomposing the correlated (fractal) component for all scaling functions (see the representative black and blue functions with their associated foci). Parametric maps for the second moment crossover scale, Hurst exponent, and full-width-at-half-maximum of the multifractal singularity strength spectrum along with the focus of the generalized Hurst exponent function are seen in (D). Note that only multimodal analysis can reveal the real topology of multifractality in the brain (D) that cannot possibly be captured by unimodal analysis (B).

Discussion

We reported here on the SFD-FMF method as a genuinely multifractal approach to decompose the scale-free constituents of empirical bimodal signals by combining our multifractal formalism (Mukli et al., ) with the use of the Bienaymé formula (Bienaymé, ). We also developed qSRA-FMF, a moment- and FMF-based variant of the segmented line regression method (Ge and Leung, ) in order to obtain an MSE-based reference for discerning (i) additive from (ii) non-additive forms of signal genesis (see Figures 2A,B, respectively). Goodness-of-fit statistics—for given bimodal characteristics—were used as a guide in choosing between estimates obtained by the SFD-FMF or qSRA-FMF methods, as respectively valid. Accordingly, based on low MSE values, when analyzing bimodal EEG, NIRS, and fMRI-BOLD data we could demonstrate that these signals resulted from superposition. When applied to high-definition empirical signals (EEG, NIRS) with high degrees of freedom, these methods performed in a robust manner. Sub-optimally sampled physiological processes—such as, the exemplary fMRI-BOLD imaging data—however, obviously imposed limitations to the extent to which these methods could reveal bimodality in the signals; a circumstance that we could overcome by reducing the degrees of freedom in the analysis, thus yielding a robust performance on limited-definition fMRI-BOLD signals too.

Physiological significance

Complex dynamics in biological systems—like that of the brain—have recently become the focus of intensive research as they represent an essential attribute for normal functioning (Bullmore et al., ). Staying with the example of the brain, multifractality is regarded as one of the main facets of complex, scale-free dynamics emerging from the underlying immense neuronal networks (Bullmore et al., ; Bullmore and Sporns, ). Multifractal characterization of scaling is inherently complex. It cannot be applied to empirical data without defining the application criteria in terms of the properties of the empirical signals in which resting-state brain dynamics are captured (i.e., EEG, MEG, NIRS, fMRI-BOLD; Eke et al., ). One such important property is multi- or bimodality that was the subject of this study.

The standard moment-based analyses of multifractal behavior, operating on the basis of an assumed unimodal model, estimates the scaling exponents within a single SR. This approach, however, will lead to erroneous estimates if unimodality does not hold. Indeed, it has been shown that EEG, NIRS, and fMRI-BOLD signals (Eke et al., ; Gifani et al., ; Herman et al., ) were in fact multimodal, a signal property that, therefore, must be taken into consideration in their multifractal characterization. To this end, our SFD method provides a means for decomposing the signal components of a bimodal signal with each part having its own set of single SRs needed to meet the aforementioned criteria for scale-invariance.

Beyond obtaining correct estimates for the scaling exponents, an understanding of the signal genesis in reference to the underlying physiological factors should be the subject of future research. Accordingly, in this work, we were motivated to develop the multifractal signal decomposition methods as needed and likely useful instruments to study multimodal signal genesis, in particular in the case of hemodynamic signals—such as, NIRS or fMRI-BOLD—that are widely used in brain connectivity research (Biswal et al., ; Mesquita et al., ). Both could be modeled as the convoluted product of ongoing regional neuronal activity (EEG) and the regional hemodynamic response function (HRF; Liu et al., ). As HRF is equivalent to low-pass filtering, it should render the hemodynamic signal multimodal with a breakpoint in its scaling function. Our qSRA and SFD methods can readily separate a breakpoint manifesting as signal convolution from a crossover resulting from signal superpositioning; we found the latter to be the case with the empirical signals studied. Despite its fundamental importance in signal genesis, this issue has not yet been addressed in the field of physiology or in particular in that of brain dynamics. Thus, our results on multifractal crossovers in the exemplary resting-state NIRS and fMRI-BOLD signals should be regarded as not only relevant but also promising in that our qSRA and SFD methods can reveal various facets of hemodynamic signal genesis in the brain.

Crossover scales

Our FMF formalism (Mukli et al., ) as implemented in SFD-FMF offers an explicit framework to deal with the crossover of empirical multifractals. Incorporating the focus in the regression scheme of the SFD-FMF method allows for a robust estimation of crossovers in empirical signals. Omitting the use of the focus in guiding the regression process of H(q)—as is the case with standard methods like MF-DFA (Kantelhardt et al., ; Ihlen, )—would potentially upset the moment-wise order of their independently obtained regression slopes. Ensuring such order of regression slopes in H(q) is essential in meeting the application criteria of the Legendre transformation incorporated in the multifractal formalism (Frisch and Parisi, ; Halsey et al., ; Bacry et al., ) and thus preventing outcomes with inversed or corrupted singularity spectra (Mukli et al., ; Delignières et al., ). To this end, Gierałtowski et al. () attempted to circumvent the use of the Legendre transformation, known to be an inherently critical step in the formalism, with the multiscale assessment of H(q) (termed multiscale multifractal analysis, MMA). Nevertheless, a “reversed” H(q,s)—i.e., H increasing with q—as shown in Figure 17 of Gierałtowski et al. () does evidently fail in meeting the above-mentioned criterion for a monotonously declining H(q,s) with increasing qs. This also explains the artifacts in the singularity spectra seen in their Figure 6. While their interpretation of the change of the average slope of the scaling function profile at some scale, s, would not suffice for multifractal analysis complete with singularity spectrum, it can still approximate breakpoints and/or crossover scales, but without discerning these phenomena. Nevertheless, there are other effective approaches to circumvent the difficulties associated with the use of the Legendere in obtaining D(h) (Jensen et al., ; Chhabra et al., ).

Impact of component focus ratio and temporal correlation

Additive random or correlated noise readily upsets multifractal analysis as demonstrated by Ludescher et al. (). In particular, when random noise with increasing amplitudes is added to a multifractal signal, the q-wise scaling function obtained by standard tools such as, MF-DFA will exhibit crossovers gradually shifting to larger scales. Accordingly, H(q) evaluated within the range of chosen scales becomes dominated by that of the added noise component. While the aim of these authors was not to provide a solution for handling noise contamination, they made an important contribution demonstrating that added noise propagates across scales ultimately leading to spurious results in multifractal analysis.

As for the impact of the component focus ratio (Figure 6), from the above-mentioned geometrical properties of multifractal scaling functions and the relationships shown in Figure 5A2 it follows that the crossover scale is low when both H and the component focus ratio are low (Figure 6, vignettes in lower left corner). Conversely, it is high when both H and the component focus ratio are high (Figure 6, vignettes in upper right corner). In between these extremes, a diagonal band of low bias due to the impact of mid-range crossover scales in the data is seen (Figure 6, bottom row) where the presence of merging scale-free patterns can be statistically confirmed (Clauset et al., ). When—due to the actual representation of scaling in the empirical data—the analysis is not performed within this optimized range, the signal definition should be improved: crossover scales that are much too low require an increased sampling rate, those that are much too high call for a longer signal to be collected (Eke et al., ).

There are cases when the superposition of two fractal components yields a composite signal with a crossover falling outside the observed range of scales (Figure 10). The multifractal spectrum in this case is typically asymmetric (Drożdż and Oświęcimka, ). While under these conditions the crossover is not directly accessible to our SFD-based analysis, our additive model still allows for its characterization and offers an explanation for the asymmetry in D(h). This way, a composite process yielding asymmetric D(h) can be modeled too. Asymmetric D(h) can also be interpreted as a phase transition (Grassberger et al., ; Muzy et al., ; Arneodo et al., ; Radons and Stoop, ) based on the deep analogy that exists between the multifractal formalism and equilibrium statistical thermodynamics (Stanley and Meakin, ; Tel, ; Arneodo et al., ). According to this model, the superimposed partition functions under and above a critical q and the signals mutually perturb each other to a slight extent (Arneodo et al., ); thus, the partition function, τ(q) and D(h) are always dominated by a single component under which condition the Bienaymé formula yields similar results in describing the superposition ().

Figure 10

Impact of moment level

In a scaling function representation of empirical temporal multifractality, the crossover scale for the chosen smallest negative moment is the largest and it becomes the smallest at the largest positive moment (see for example Figure 5B2). This moment-wise distribution of crossover scales emerges from the geometrical underpinnings of FMF (see Figure 2C) and the way H(q)-dependence is formulated in Equation (12) yielding the crossover scale itself. As crossover scales and breakpoints are similar manifestations of scaling, breakpoints should also be captured in a moment-wise manner.

The significance of the breakpoint in the analysis of bi- or multimodal signals has already been recognized in the literature (Peng et al., ; Kantelhardt et al., ; Eke et al., ; Herman et al., ; Ge and Leung, ). However, disregarding their moment-wise nature (Ge and Leung, ) distorts the acquired H(q)s with breakpoints falling outside the proper SR. This generates uneven error propagation across the moments in the case of standard multifractal methods (Ludescher et al., ), or introduces a slight but global error when FMF methods are applied. To the best of our knowledge, the qSRA- and SFD-FMF methods should be regarded as first attempts to carry out genuinely adaptive multifractal analyses on bi- or multimodal signals in a properly designed, moment-wise manner. They can readily be combined with various fractal and multifractal tools as their adaptive step, thus opening new possibilities for future applications.

Significance of the fGn-fBm framework

Mono- and multifractal analyses alike have been shown to benefit from the fGn-fBm fractal signal model of Mandelbrot and Van Ness () as implemented by Eke et al. (). Despite its simplicity, this model captures the most fundamental properties of scale-free signals as realizations of fGn and fBm processes. These two signal classes are mutually convertible via cumulative summation or differencing of their subsequent values leading to an increase or decrease in their extended H by 1, respectively. As seen in Figure 11, methods differ in the range in which they yield minimal-bias estimates of extended H across the fGn-fBm framework. A multimodal signal may well contain components with H falling above or below the minimal-bias range. Hence, signal classification by methods such as, SSC (Eke et al., ) should be incorporated as a first step in the analysis. The actual implementation is by repetitive signal conversion until a difference in two successive H estimates is found of ~1. All of our signals and signal components proved one or other of the two classes of the fGn-fBm dichotomy. This way of handling signal classes in fractal analyses should reconcile issues seen with other adaptive approaches (Kuznetsov et al., ).

Figure 11

Comparing overall performances and limitations of qSRA and SFD methods

The precision of the qSRA method increases with the level of tolerance, which in turn results in contracted SRs. This tends to weaken the estimates of H(q) due to falling short of securing wide enough scale-invariance (as seen in Figure 1). This effect is altogether eliminated by the SFD method, which makes use of all the data of the merging signal components.

Our SFD method was validated against synthetic signals (Figures 6, 7, and Table 1). It outperformed the qSRA method for low-scale crossovers as the latter was shown to be susceptible to increased fluctuations typically seen in the large-scale region with limited number of available non-overlapping windows (Cannon et al., ). Hence, qSRA tends to locate the breakpoint for spuriously high scales (Figure 6). Multimodal scaling functions may deviate from the exact segmented line regression model near the breakpoint (Kuznetsov et al., ). This phenomenon always appears in superimposed fractals or multifractals, and can be modeled by added scaling functions (Figure 12).

Figure 12

Significance of the design concept

The SFD approach is built around the notion that the multimodality emerges from the superposition of multiple and typically scale-free signal components. Multimodal multifractal scaling functions can also be produced by non-fractal generators like the infinitely divisible cascades (Chainais, ), filters, etc. Only specific generator models can adequately treat scale-dependent dynamics in the multimodal processes. As such, causal models are typically unavailable; our qSRA method can offer a phenomenological solution. Thus, where a causal generator model is sought, analysis by SFD should be preferred over that by qSRA. As the number of modeling parameters are the same in both cases, a smaller level of goodness-of-fit statistics (i.e., MSE, SSE) can indicate which of the two seems a realistic model for a particular signal (for a demonstration see Figures 7, 9). Performance of our methods and their discriminating power evidently depend on empirical signal definition (length, sampling, non-fractal contamination, etc.) as seen in the cases of our exemplary applications to high-definition EEG and NIRS, and limited-definition fMRI-BOLD signals.

Performance of qSRA and SFD methods on empirical signals

Human EEG and NIRS signals

The crossover between the EEG signal components was found at the boundary between the δ and θ bands (Figure 8) of EEG classification. An independent δ and θ rhythm has already been proposed due to the significant interregional gap in synchrony (Mormann et al., ; Fetterhoff et al., ). Our current findings also prove the presence of an independent multifractal δ rhythm. The NIRS signal was also found to be the sum of temporally correlated, scale-free fluctuations dominated by low frequencies (i.e., low-frequency fluctuations) and uncorrelated [but not instrumental (Eke et al., )] noise.

Rodent fMRI-BOLD imaging data

Resting-state brain dynamics as captured in fMRI-BOLD fluctuations is powered by ongoing neurodynamics spreading across the functional connections of a fractally organized anatomical network of an immense neuronal pool (Bullmore et al., ; Werner, ). Previously, we demonstrated that the resting-state fMRI-BOLD signal in an animal (Herman et al., ) and the fNIRS signal in a human model (Eke et al., ), alike, were bimodal temporal fractals. Thus, subject to future but likely advances in fMRI-BOLD technology, an improved signal definition may likely reveal bimodality in human fMRI-BOLD signals too. However, due to limitations inherent to the current fMRI-BOLD technology (Eke et al., ), their scale-free character—especially in cases of bimodality—is hard to capture (see Figure 12). Accordingly, while the SFD analysis yielded valid estimates of crossovers in a sizable fraction of the scanned voxels, still—depending on the correlation level and component focus ratio—numerous crossover estimates—due to sub-optimal data acquisition—proved invalid at extreme scales. This severely upset the topology of the multifractal parametric maps and prompted us to look for justifiable grounds to reduce the degree of freedom in finding breakpoints and/or crossovers in an attempt to improve the performance of the analysis. We did indeed find that the uncorrelated components of the scaling functions could be taken as default throughout the voxels of the scan. This uncorrelated component was determined at a site (Figure 9) where the precision of the method has been shown to be maximal in tests on synthetic signals (Figure 6). This procedure effectively compensated for the limited definition in the BOLD signals as seen in the enhanced topology of the parametric maps (Figure 9D).

Conclusions and future perspectives

The issue of bimodality presents a major challenge when it comes to multifractal analysis of complex biological signals. We reported a novel approach (SFD-FMF method) as a genuinely multifractal tool to decompose the scale-free components of empirical bimodal signals by combining our multifractal formalism (Mukli et al., ) with the use of the Bienaymé formula (Bienaymé, ). We also developed a moment- and FMF-based variant of the segmented line regression algorithm (qSRA-FMF method) to discern additive from non-additive forms of signal genesis based on respective goodness-of-fit statistics. When applied to high-definition empirical signals (EEG, NIRS), these methods performed in a robust manner. The performance on sub-optimally sampled physiological signals is expected to be weaker; a circumstance that we could overcome by reducing the degree of freedom of the analysis, thus restoring a robust performance of the SFD-FMF and qSRA-FMF methods on a limited-definition fMRI-BOLD imaging dataset too. These methods offer the means to identify signal generators in physiological processes. Most importantly, they open ways to characterize the topology of key multifractal metrics in the brain emerging from its complex network dynamics.

Statements

Author contributions

ZN developed the method and wrote the manuscript. PM performed numerical tests and analyzed empirical datasets for demonstration purposes. PH provided fMRI BOLD scans for demonstrational purposes. AE helped developing and writing the manuscript and provided conceptual guidance in the study.

Acknowledgments

The authors acknowledge the use of the EEG data provided as a downloadable file last accessed on February 17, 2015, at https://sites.google.com/site/projectbci by A. Yazin of National University of Sciences and Technology, Islamabad, Pakistan. The authors declare no conflict of interest.

Conflict of interest

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

Supplementary material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2017.00533/full#supplementary-material

    Symbols and definitions

  • ∧

    estimated value

  • A

    amplitude

  • β

    spectral index

  • ΔH15

    – the difference between the H(−15) and H(15) values

  • DFA

    detrended fluctuation analysis

  • DHM

    davies and harte method

  • DMA

    detrending moving average

  • D(h)

    spectrum of singularity strength (or singularity spectrum for short)

  • EEG

    electroencephalography

  • ER

    exclusion range

  • f

    frequency (Hz)

  • fBm

    fractional Brownian motion (non-stationary signal)

  • fGn

    fractional Gaussian noise (stationary signal)

  • FMF

    focus-based multifractal formalism (an approach using a focus-based regression scheme)

  • fMRI

    functional magnetic resonance imaging

  • fNIRS

    functional near-infrared spectroscopy

  • f

    measure describing the fractal signal component dominating over the higher scales

  • fwhm

    full width of the singularity spectrum, D(h), at half of its maximum

  • H

    extended Hurst exponent (fGn: 0<H<1; fBm: 1<H<2; summed fBm: 2<H<3)

  • Htrue

    known value of the Hurst exponent in numerical syntheses of time series

  • H(q)

    generalized Hurst exponent

  • h

    Hölder exponent

  • hmax

    the value of h at the peak position of D(h)

  • HRF

    hemodynamic response function

  • m

    detrending order (in DFA)

  • MF

    multifractal (an approach using a standard regression scheme)

  • MSE

    mean squared error

  • μ

    measure

  • N

    the length of time series (in data points)

  • n

    measure describing the fractal (or noise) dominant over the lower scales

  • Nc

    number of constituent signals

  • NIRS

    near-infrared spectroscopy

  • Ns

    number of non-overlapping segments

  • q

    statistical moment order used in multifractal analysis (moment for short)

  • qSRA

    moment-wise scaling range adaptivity (method)

  • s

    temporal scale

  • s′

    scaling boundary (possible breakpoint)

  • sb

    breakpoint

  • sx

    crossover scale

  • SD

    standard deviation

  • S[Xi](q,s)

    scaling function value at a given q and s calculated from signal Xi

  • S[Xi](N)

    the focus of the scaling function for signal Xi

  • SR

    scaling range

  • SFD

    scaling function decomposition (method)

  • SSC

    signal summation conversion (method)

  • SSE

    sum of squared error

  • SSM

    spectral synthesis method

  • v

    the order of non-overlapping segments v = 1,…,Ns

  • x

    the scale selected for obtaining an actual value of a scaling function

  • Xi

    time series (signal), where i = 1,…,N

  • WL

    wavelet leader (method)

  • WTMM

    wavelet transfer modulus maxima (method)

References

Summary

Keywords

multifractality, focus-based multifractal analyses, multimodality, breakpoint, crossover, NIRS, EEG, fMRI-BOLD

Citation

Nagy Z, Mukli P, Herman P and Eke A (2017) Decomposing Multifractal Crossovers. Front. Physiol. 8:533. doi: 10.3389/fphys.2017.00533

Received

24 April 2017

Accepted

10 July 2017

Published

26 July 2017

Volume

8 - 2017

Edited by

Zbigniew R. Struzik, University of Tokyo, Japan

Reviewed by

Stanislaw Drozdz, Institute of Nuclear Physics (PAN), Poland; Alain Arneodo, University of Bordeaux 1, France

Updates

Copyright

© 2017 Nagy, Mukli, Herman and Eke.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Andras Eke

This article was submitted to Fractal Physiology, a section of the journal Frontiers in Physiology

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Share on WeChat

Scan with WeChat to share this article

Article metrics