Impact of Healthy Aging on Multifractal Hemodynamic Fluctuations in the Human Prefrontal Cortex

Fluctuations in resting-state cerebral hemodynamics show scale-free behavior over two distinct scaling ranges. Changes in such bimodal (multi) fractal pattern give insight to altered cerebrovascular or neural function. Our main goal was to assess the distribution of local scale-free properties characterizing cerebral hemodynamics and to disentangle the influence of aging on these multifractal parameters. To this end, we obtained extended resting-state records (N = 2¹⁴) of oxyhemoglobin (HbO), deoxyhemoglobin (HbR) and total hemoglobin (HbT) concentration time series with continuous-wave near-infrared spectroscopy technology from the brain cortex. 52 healthy volunteers were enrolled in this study: 24 young (30.6 ± 8.2 years), and 28 elderly (60.5 ± 12.0 years) subjects. Using screening tests on power-law, multifractal noise, and shuffled data sets we evaluated the presence of true multifractal hemodynamics reflecting long-range correlation (LRC). Subsequently, scaling-range adaptive bimodal signal summation conversion (SSC) was performed based on standard deviation (σ) of signal windows across a range of temporal scales (s). Building on moments of different order (q) of the measure, σ(s), multifractal SSC yielded generalized Hurst exponent function, H(q), and singularity spectrum, D(h) separately for a fast and slow component (the latter dominating the highest temporal scales). Parameters were calculated reflecting the estimated measure at s = N (focus), degree of LRC [Hurst exponent, H(2) and maximal Hölder exponent, h_max] and measuring strength of multifractality [full-width-half-maximum of D(h) and ΔH₁₅ = H(−15)−H(15)]. Correlation-based signal improvement (CBSI) enhanced our signal in terms of interpreting changes due to neural activity or local/systemic hemodynamic influences. We characterized the HbO-HbR relationship with the aid of fractal scale-wise correlation coefficient, r_σ(s) and SSC-based multifractal covariance analysis. In the majority of subjects, cerebral hemodynamic fluctuations proved bimodal multifractal. In case of slow component of raw HbT, h_max, and Ĥ(2) were lower in the young group explained by a significantly increased r_σ(s) among elderly at high temporal scales. Regarding the fast component of CBSI-pretreated HbT and that of HbO-HbR covariance, h_max, and focus were decreased in the elderly group. These observations suggest an attenuation of neurovascular coupling reflected by a decreased autocorrelation of the neuronal component concomitant with an accompanying increased autocorrelation of the non-neuronal component in the elderly group.

Introduction

Scale-free dynamics is an ubiquitous property of physiological processes (West, 1991; Eke et al., 2000; Ivanov et al., 2001) such as low frequency fluctuations of cerebral hemodynamics (Fox and Raichle, 2007; Herman et al., 2009; Pierro et al., 2012) and neural activity (Linkenkaer-Hansen et al., 2001; Ivanov et al., 2009; He, 2014). Scale-free dynamics is a hallmark of complexity viewed as an emergent property of biological systems composed of numerous elements with a network of stochastic (typically weak) connections amongst them (Csermely, 2006). Several human studies investigated the scale-free phenomenon of functional brain imaging signals by using mono-(Eke and Hermán, 1999; Thurner et al., 2003; Maxim et al., 2005; Eke et al., 2006; Khoa and Nakagawa, 2008; Bullmore et al., 2009; He et al., 2010; He, 2011; Herman et al., 2011) and multifractal analysis (Shimizu et al., 2004; Wink et al., 2008; Ciuciu et al., 2012; Quang Dang Khoa and Van Toi, 2012). Monofractal analysis reveals global, long-range correlation (LRC) structuring in terms of the influence of past events in the process on future ones (Bassingthwaighte et al., 1994; Eke et al., 2000, 2002). Multifractal analysis yields a distribution of fractality measures (Barabási et al., 1991; Stanley et al., 1999; Kantelhardt et al., 2002; Ihlen and Vereijken, 2010; Mukli et al., 2015) that enables a more detailed characterization of local temporal scaling behavior provided that fluctuations at wide range of temporal scales are sufficiently represented in the sampled physiological process (Eke et al., 2012). The estimation of these complexity parameters is essentially based on a power-law model fitted to the appropriate statistics of the data, which is reliable only if sample size is large enough [at the order of hundreds, Eke et al., 2002; Clauset et al., 2009]. Such statistics usually shows power-law scaling—indicating LRC—within a bounded interval of temporal scales usually termed as scaling range (SR; Caccia et al., 1997). In addition to the finite representation of the dynamics the lower and upper boundary of SR are determined by the signal genesis of the underlying physiological process. Nevertheless, numerous examples of empirical data exhibit multiple SR indicating multimodal scaling, see examples cited in Nagy et al. (2017). Multimodality has also been of concern in case of cerebral hemodynamics and typically present itself with two (case of bimodality) or even more apparent SRs in which the statistical measure of fractal analysis scales with a different exponent (Nagy et al., 2017).

Application of a possible bimodal scale-free model on resting-state hemodynamic signals requires a measurement technology, which assures that the process is sampled at high enough rate in long enough records. Near-infrared spectroscopy (NIRS) is an emergent imaging technology which readily captures cerebrocortical resting-state hemodynamic fluctuations at a cm spatial resolution and at high sampling frequency with no particular limitations on signal length (Jöbsis, 1977; Chance et al., 2007; Fox and Raichle, 2007). In case of continuous wave near-infrared spectroscopy (cwNIRS), the measured intensity signals are determined by the relative tissue concentration of total hemoglobin (HbT) and its constituents: oxy- and deoxyhemoglobin (HbO and HbR, respectively; Cope et al., 1988).

By now the physiological underpinnings of the functional NIRS (fNIRS) signal has been elucidated (Jöbsis, 1977; Kocsis et al., 2006a; Tachtsidis et al., 2008). As to its dynamics, oscillations of cerebral hemodynamics has been characterized by spectral analysis of NIRS signals (Elwell et al., 1999). Monofractal pattern of NIRS spectral data was first reported by Eke and Herman for the human brain cortex (Eke and Hermán, 1999). Later, multifractal properties of fNIRS were also demonstrated (Quang Dang Khoa and Van Toi, 2012). These pioneering reports understandably focused on signal analysis and not making an attempt to identify the underlying physiological mechanisms shaping the observed complex patterns. As to the nature of local hemodynamic fluctuations, they are primarily elicited by neural activity via neurovascular coupling (NVC; Devor et al., 2003; Drake and Iadecola, 2007) but the hemodynamic response is also modulated by endothelial mechanism (Li et al., 2013; Chen et al., 2014). In addition, systemic hemodynamic effects should be considered (Yamada et al., 2012; Scholkmann et al., 2014) for an enhanced interpretation of results obtained from resting-state records from subjects with similar age. Apart from non-biological noise and motion artifacts, resting-state NIRS (rsNIRS) signal bears the fingerprint about systemic hemodynamics such as cardiac cycle and respiration (Tian et al., 2009; Li et al., 2013). Separation of the functional and the systemic components became of considerable interest and various approaches have been developed to address this issue (Scholkmann et al., 2014). Correlation between the fluctuations of oxy- and deoxyhemoglobin concentration is the basis of signal improvement presented by Cui et al. (2010) and the technique proposed by Yamada et al. (2012). Under certain assumptions—that usually hold in resting state, neural activity renders the relationship between HbO and HbR fluctuations more anticorrelated while fluctuations of systemic origin in the resting state (Cui et al., 2010; Scholkmann et al., 2014) has a correlated influence on hemoglobin chromophores.

The nonstationary fractal character of HbT (Eke et al., 2006) implies that its constituents, HbO, and HbR, also exhibit non-stationary dynamics. Consequently, their relationship should be explored in terms of a non-stationary characterization of correlation. Therefore, the HbO-HbR relationship was studied with the aid of scale-wise fractal cross-correlation coefficient and a novel adaptation of multifractal covariance analysis. The former is essentially a measure applicable to nonstationary time series building on the correlation of scale-wise mean variances obtained separately for HbO and HbR signal (Podobnik and Stanley, 2008; Zhou, 2008), while the latter is a multifractal approach examining the scaling properties—and its moments of various order—of HbO-HbR covariance similarly to the analyses described in Refs (Kristoufek, 2011; Jiang et al., 2016; Zhao et al., 2017).

It has been shown that aging (Goldberger et al., 2002; Lipsitz, 2003) and various diseases (Ivanov et al., 1999; Goldberger et al., 2002; Maxim et al., 2005) affect complexity parameters and the impact of other factors, such as gender (Ni et al., 2014), have also been recognized. This study contributes to this accumulating body of knowledge on the influence of aging on the complexity of physiological processes. In general, the contraction of homeodynamic space is an essential attribute of an aging biological system meaning that the dynamics of physiological processes in an elderly person is typically confined to a restricted state space (Rattan, 2014). Taking the cardiovascular system as an example the well-known dependence of heart rate variability (Beckers et al., 2006; Vandeput et al., 2012) on age can be attributed to a decline in autonomic modulation (Nunes Amaral et al., 2001; Lipsitz, 2003; Tulppo et al., 2005; Silva et al., 2017) reducing the adaptational reserve of regulatory mechanisms (Goldberger et al., 2002). Such reports on aging and altered complexity parameters are typically based on demonstrating coincidences, while there is a palpable need to establish a causal relationship for the changing complexity.

Accordingly, our goal was to provide plausible explanation for the physiological mechanism of observed age-related alterations based on parameters of complexity obtained from measures of brain hemodynamics captured by rsNIRS. Firstly, the measured signals were evaluated for the presence of true LRC-type multifractality. Multifractal parameters obtained for young and elderly volunteers were compared to assess the impact of aging on the complexity of cerebral hemodynamics. In addition, we extended our analysis incorporating the oxy- and deoxyhemoglobin signals in order to explore the age-dependent alterations in their relationship. The influence of age was characterized by exploring the strength of causal link between these measures of the coupled fluctuations of oxy- and deoxyhemoglobin and multifractal parameters of cerebral hemodynamics.

Methods

Extended records of fluctuating rsNIRS signals from the human brain prefrontal cortex (PFC) analyzed for their multifractality in this work have been collected in a previous study of the group reporting on the monofractal serial correlation in these signals (Eke et al., 2006).

Near Infrared Spectroscopy

According to the principle of cwNIRS, backscattered light intensities were measured at wavelengths of 775, 830, 849, and 907 nm by a NIRO 500 Cerebral Oxygen Monitor (Hamamatsu Photonics, Hersching, Germany), a single-channel instrument. The mean penetration depth of near infrared photons for the 4 cm interoptode distance of this device can be estimated at approximately 2 cm (Firbank et al., 1998), thus our NIRS optode sampled the brain cortex (Chance, 1994). Based on the modified Beer-Lambert law (Kocsis et al., 2006b) the relative tissue concentrations of HbO and HbR were calculated along with HbT obtained as their sum. While the fluctuation of HbT reflects on cerebral blood volume (CBV) dynamics, that of the other two chromophores and their relationship are dependent on oxygenation, too.

Data Collection

HbO, HbR, and HbT data were dumped via the RS232 port of the NIRO instrument into a computer file at a rate of 2 Hz. Extended records of HbO, HbR, and HbT samples were created and processed for each subject in length of N = 2¹⁴ proven to be adequate for fractal analysis (Eke et al., 2002). The source and detector fibers were secured in a rubber pad. The optode was mounted just under the hairline over the forehead. The cranium was shielded from ambient light. Instrument noise was determined by placing the optode over a slab of “mock” brain, whose scattering (μ_s = 10.96 1/mm) and absorption (μ_a = 0.099 1/mm) coefficients were matched to that of the human brain (courtesy of Prof. Britton Chance, University of Pennsylvania, Philadelphia, U.S.A.). The power of instrument noise was found negligible when compared to the power of resting-state fluctuations in vivo (Eke et al., 2006).

Subjects

Following approval by the Local Research Ethics Committee of Semmelweis University and having obtained informed written consent, 52 healthy volunteers with no neurological disorders (28 women, 24 men) participated in the study. To evaluate the effect of age and gender, subjects were assigned into groups of young females (Y_F, n = 9, age < 45 years), young males (Y_M, n = 13, age < 45 years) elderly females (E_F, n = 19, age≥45 years) and elderly males (E_M, n = 10, age≥45 years). The rsNIRS measurements were carried out in a comfortable sitting position in a session slightly exceeding 2.5 h as previously described in (Eke et al., 2006).

Data Preprocessing

Multifractal analyses were performed on raw signals and following correlation-based signal improvement (CBSI; Cui et al., 2010; Scholkmann et al., 2014), the latter designed to remove correlated (systemic) influences (e.g., head movement). Specifically, this preprocessing step incorporate a standard deviation ratio of HbO and HbR that is used in their subsequent linear combination yielding: ${}_{}^{CBSI}HbT=\frac{1}{2}\text{ }\cdot \left(HbR-\left({}^O\sigma \left(N\right)/{}^R\sigma \left(N\right)\right)\cdot HbR\right)+\frac{1}{2}\text{ }\cdot$
$\left(HbO-\left({}^R\sigma \left(N\right)/{}^O\sigma \left(N\right)\right)\cdot HbO\right).$ The improved HbT signal is regarded as a representation of anticorrelated hemodynamics attributable to local hemodynamic influences and oxygen consumption accompanying neuronal activity.

Multifractal Analyses

Multifractal Scaling Analysis

Signal summation conversion (SSC) method

The multifractal scaling functions of HbT, HbO, and HbR were calculated by the multifractal SSC method (Mukli et al., 2015) as the basis of our approach to evaluate the multifractality of our bimodal rsNIRS signals (Figure 1). For detailed description of the MF-SSC method the reader is referred to Refs. (Mukli et al., 2015; Nagy et al., 2017). Briefly, multifractal SSC uses a measure depending on the temporal scale (s) of signal window—bridge-detrended variance σ(s), see Ref. (Eke et al., 2000)—and a set of q-order statistical moments, to create corresponding moment-wise generalized variance profiles, S_σ(q, s) of Equation 1, spanning across a range of temporal scales within the chosen analytical SR.

$\begin{array}{ll}{S}_{\sigma }(q,s)={\left[\frac{1}{N_s}{\displaystyle \sum _{v=1}^{N_s}}{\mu }_v^q(s)\right]}^{1/q}\propto s^{H(q)}& (1)\end{array}$

FIGURE 1

Figure 1. Various representations of the measured NIRS signals. (A) Resting-state raw NIRS signal components. (B) Variance profiles of different signal components. (C) Variance profiles of HbT calculated at different set cross-correlation levels. (D) Fractal correlation coefficient between HbO and HbR as a function of scale. In this paper multifractal scaling function values are actually given by S_σ(q, s) [or by S_σ(s) for q = 2], where the subscript σ refers to the measure of the chosen multifractal method (SSC).

Specifically, 60 logarithmically spaced time scales were chosen between s_min = 16 and s_max = 8192 for computing σ_v(s) in each non-overlapping time window [v = 1, 2, …, N_s = int(N/s)] of cumulatively summed signal. The low temporal scales dominated by fast fluctuations with weak, non-fractal autocorrelation (Eke et al., 2006) were excluded and the fairly high scales (Nagy et al., 2017) with well-manifested scale-free processes were secured. The selected moment orders ranging from −15 to 15 were adequate¹ to capture both large- and small-variance dynamics in a fluctuating rsNIRS signal; the former being emphasized by variance profiles corresponding to positive moments, the latter by those corresponding to negative moments (Grech and Pamuła, 2012; Mukli et al., 2015). Please note that statistical estimates obtained at q = 2 obey rules of linearity, while in fact, those for q≠2 capture non-linearity in system dynamics (Gómez-Extremera et al., 2016; Bernaola-Galván et al., 2017).

Relationship between variance profiles of hemoglobin chromophores

Since HbT = HbO+HbR, it is the generalized Bienaymé formula (Equation 2) which describes the relationship between their variance profiles²: ^TS_σ(2, s), ^OS_σ(2, s), and ^RS_σ(2, s) in an exact form³ (Nagy et al., 2017).

$\begin{array}{ll}{}^{T}S{}_{\sigma }{\left(2,s\right)}^2={}^{O}S{}_{\sigma }{\left(2,s\right)}^2+{}^{R}S{}_{\sigma }{\left(2,s\right)}^2+2r_{\sigma }\left(s\right)·{}^{O}S{}_{\sigma }\left(2,s\right)·{}^{R}S{}_{\sigma }\left(2,s\right).& (2)\end{array}$

The unknown factor is the scale-wise fractal correlation coefficient, r_σ(s); all the others are scaling function values directly calculated from the measured signal.

Focus-Based Scaling-Range Adaptive Analyses

Focus-based multifractal analysis (FMF-SSC) were carried out (Mukli et al., 2015) to estimate the generalized Hurst exponent, H(q), which is essentially a set of slopes of the scaling function profiles fitted in the analytical SR, as expressed by

$\begin{array}{ll}{S}_{\sigma }(q,s)\propto s^{H(q)}.& (3)\end{array}$

H(q) describes the moment-wise or global distribution of fractal correlation (essentially that of the fractal dimension) in the signal, thereby generalizing H(2), the usual outcome of monofractal time series analysis. Taking H(q) as its input, the multifractal formalism (Frisch and Parisi, 1985; Barabási and Vicsek, 1991; Muzy et al., 1993; Eke et al., 2012) via the multiscaling exponent, τ(q) = qH(q)−1, and Legendre transformation will yield D(h), the local distribution of the fractal dimension or singularity strength, D, as a function of roughness or Hölder exponent, h.

$\begin{array}{llll}h& =& \frac{d\tau (q)}{dq},& (4)\end{array}$

$\begin{array}{ll}D\left(h\right)=\underset{q}{\text{inf}}\left(qh-\tau (q)\right).& (5)\end{array}$

Incorporating the focus—obtained as a fitted intersection of scaling function profiles at s = N–ensures monotonously decreasing Ĥ(q) and thus stable, uncorrupted D(h). Enforcing the focus of scaling function, ln(Ŝ_σ(N)), when regressing for Ĥ(q) was recognized as a prerequisite to obtain stable results of multifractal analysis (Mukli et al., 2015).

In an attempt to find the best fitting model for the observed scaling functions we adapted the concept of bimodality originally described in the frequency domain (Eke et al., 2006). This pattern can be recognized as two scale-free processes separated by moment-wise breakpoint scales (s'(q)) (Nagy et al., 2017). The less correlated (“fast”) component dominates the lower range of scales while a more correlated (“slow”) component is characteristic in the higher scales (Figure 2A). A breakpoint-based bilinear regression model (denoted as moment-wise SR adaptive) was implemented as described in Ref. (Nagy et al., 2017). Briefly, it is an iterative process by estimating breakpoint scales that minimize sum of squared error (SSE) of the residuals for each and every moment as

$\begin{array}{l}SSE(s\prime (q))={\displaystyle {\sum }_{q=-15}^{15}[{\displaystyle {\sum }_{s=s_{\min }}^{s^{\prime }(q)}(}{}^f\widehat{H}(q)·(\ln s-\ln N)}\\ +\ln {}^f\widehat{S}{}_{\sigma }(N)-\ln S_{\sigma }(q,s){)}^2+({}^f\widehat{H}(q)·(\ln s-\ln N)\\ +\ln {}^s\widehat{S}{}_{\sigma }(N)-\ln S_{\sigma }(q,s){)}^2],& (6)\end{array}$

where ln(^sŜ_σ(N)) and ln(^fŜ_σ(N)) are the moment-independent foci (Figure 2B) associated with the slow and fast components, respectively. ^sĤ(q), ^fĤ(q) are the slopes (Figure 2C) of the fitted two lines of regression⁴ (i.e., the generalized Hurst exponent functions of the two components). This iteration thus adaptively yields the best segregation of scaling ranges. $\widehat{H}$(q) and D(h) (Figure 2D) are obtained for the fast and the slow component, respectively.

FIGURE 2

Figure 2. Steps of bimodal multifractal SSC analysis. (A) Scaling function of SSC as moment-wise generalization of variance profiles. Separate components are marked in blue (fast) and red (slow). (B) The sets of power-law functions fitted separately for the two components with focus-based regression. (C) Generalized Hurst exponent functions, H(q)s of the two components as functions of moment order q. The focus point and H(q) for both components (f , fast; s, slow) were iterated for finding the scale with minimum SSE(s'(q)) as the true breakpoint at a given moment, which process finally yields ln (^sŜ_σ(N)), ln (^fŜ_σ(N)), ^sĤ(q), and ^fĤ(q) with the best fit. (D) Singularity spectra of the two components. Multifractal endpoint parameters: the highest singularity strength (D = 1) is associated with the maximal Hölder exponent (h_max), which usually correlates well with H(2), a measure of global LRC in the signal. Distribution of local scale-free behavior is captured in ΔH₁₅ = H(−15) – H(15) and in full-width at half maximum (fwhm) of D(h) respectively, reflecting degree of multifractality.

We calculated global multifractal endpoint parameters to characterize the degree of autocorrelation [maximal Hölder exponent, h_max and monofractal Hurst exponent, $\widehat{H}$(2)] and multifractality [ΔH₁₅ = H(−15)–H(15), and full-width at half maximum (fwhm) of D(h) (Wink et al., 2008; Grech and Pamuła, 2012)] in the measured cerebral hemodynamic signals as illustrated in lower panels of Figure 2.

Evaluating True Multifractality

Since multifractal tools—including FMF-SSC—readily produce seemingly realistic values for multifractal measures such as D(h) even in the absence of true multifractality (Kantelhardt et al., 2002; Clauset et al., 2009; Grech and Pamuła, 2012); it is mandatory to evaluate our empirical signals in this regard using the following framework. Accordingly, because our FMF-SSC method always produces uncorrupted D(h) irrespective whether the signal is a true multifractal or not, this property needs to be tested separately (Figure A1). Verification of true multifractality consists of three subsequent steps: (i) identifying general scale-free behavior with power spectral density (PSD) analysis, (ii) distinguishing true multifractality from multifractal noise, and (iii) determining the origin of the expressed multifractal scaling. Therefore, in these tests, true, long-range correlated multifractals are to be distinguished from processes lacking scale-free properties or not showing autocorrelation. Failing to pass in any of the aforementioned tests resulted in the exclusion of the subject in question from further analysis. Details of this framework are explained in the Appendix (Supplementary Material) and in Ref. (Racz et al., 2018).

Characterizing HbO-HbR Relationship

Scale-wise fractal cross-correlation coefficient

One approach to assess the relationship of HbO and HbR fluctuations is to calculate a measure of cross-correlation by using variance profiles. After rearrangement of Equation 2 it is possible to express r_σ(s):

$\begin{array}{ll}{r}_{\sigma }\left(s\right)=\frac{{}^{T}S{}_{\sigma }{\left(2,s\right)}^2-{}^{O}S{}_{\sigma }{\left(2,s\right)}^2-{}^{R}S{}_{\sigma }{\left(2,s\right)}^2}{2·{}^{O}S{}_{\sigma }\left(2,s\right)·{}^{R}S{}_{\sigma }\left(2,s\right)}.& (7)\end{array}$

This measure indicates whether the fluctuations at a given scale are enhanced (r_σ(s) > 0), diminished (r_σ(s) < 0) by coupling HbO and HbR signal or their relationship is insensitive to coupling between them (r_σ(s) = 0).

The fractal cross-correlation analysis yielding r_σ(s) is strikingly similar—in terms of calculation steps and order—to the detrended cross correlation analysis, the major difference of these methods concerns their measure (Podobnik and Stanley, 2008). The cited approach uses fluctuation while ours calculates bridge-detrended variance to assess correlation of coupled non-stationary time series. When comparing the above scale-wise approaches with the standard means of calculating cross-correlation (i.e., Pearson and Spearman), the first major difference is that the prerequisite of stationarity applies to the latter methods. Furthermore, the sequence of calculation steps is critical, too: because when the standard cross-correlation is calculated it is followed by a step of averaging effectively abolishing the scale-wise information. It is worth of noting that Spearman proved a more robust standard measure of correlation than the Pearson coefficient as the latter assumes not only stationarity but linearity, too (Zimeo Morais et al., 2018).

Multifractal covariance analysis

The other approach is based on the extension of FMF-SSC in order to analyze the scaling of long-range cross-correlation. This method assesses the multifractality emerging genuinely from coupled oxy- and deoxyhemoglobin fluctuations. Scale- and moment-wise bridge-detrended covariances (_Cov) were calculated of HbO and HbR signals to yield an estimate of bivariate generalized Hurst exponent function, ^ORH(q).

$\begin{array}{ll}{}_{}^{OR}{S}_{Cov}\left(q,s\right)={\left[\frac{1}{N_s}{\displaystyle \sum _{v=1}^{N_s}{}^{OR}\left|Cov\right|{}_v^q(s)}\right]}^{1/q}\propto s{}^{OR}H\left(q\right)& (8)\end{array}$

Now applying Equations 3–5 yields the corresponding singularity spectrum and multifractal endpoint parameters. Covariance truly scales only if the ^ORĤ(q) function differs from the average of ^RĤ(q) and ^OĤ(q). Therefore this comparison must be carried out after obtaining the distribution of scaling exponents and output parameters of multifractal analysis. Prior to that, moment-wise bimodal regression analysis had been performed on the q-wise (generalized) variance profiles of HbO, HbR, and HbT and on the HbO-HbR covariance profile in the same manner.

Descriptive Statistical Analyses

Normal distributions in each independent sample were checked by Shapiro-Wilk's test. Difference between group means or medians were considered significant in case of p < α_s, α_s = 0.05 (level of significance). Homogeneity of variances was confirmed by Levene's test.

Assessing the Effect of Age and Gender

Two-way univariate ANOVA were performed treating multifractal endpoint parameters and r_σ(s) as dependent variables and assuming that there was an interaction between age and gender (categorical factors). Most of the results presented in this paper are based on the output of two-way ANOVA with Tukey's post hoc test. Had the prerequisites of ANOVA been not met, group means were compared with two-sample t-test (with Welch's correction for inhomogeneous variances in case of significant Levene's test). In the absence of normal distribution, the comparison of group medians was performed by Mann-Whitney U test.

Multivariate ANOVA (MANOVA) was performed to detect age- and gender-related differences between direct descriptors of the scaling function: ^sĤ(2), ^fĤ(2), ln(^sŜ_σ(N)), ln(^fŜ_σ(N)). As changes occurring by coincidence is of concern, MANOVA tests were used taking Ĥ(2), ln(^sŜ_σ(N)), ln(^fŜ_σ(N)) and r_σ(s) (at scales corresponding to the slow component s = 4,344, 4,828, 5,367, and 5,965) as dependent variables. The p < 0.05 of Wilks's test suggests that the same subjects are responsible for each of the between-group differences. Finally, examining the degree of correlation of dependent variables (expressed as r²) enables to identify the most relevant endpoint parameters of multifractal analysis.

Explaining the Variance Profiles of HbT

In order to identify a link between the altered scale-wise fractal cross-correlation coefficient and the altered multifractal endpoint parameters, the influence of r_σ(s) as an independent variable were evaluated on ^TS_σ(2, s) in a general linear model (GLM). The Bienaymé formula expresses an explicit relationship between ^TS_σ(2, s) and the scaling function values obtained for HbO, HbR, and r_σ(s) (Equation 7). Accordingly, the variability of ^TS_σ(2, s) were explained with the aid of multiple regression tests that were performed for each temporal scales. The regressors of this model were ^OS_σ(2, s), ^RS_σ(2, s), and r_σ(s) but categorical predictors (age and gender) were not included.

Subsequently, – in addition to the regressors of the previously described test—we took into account the effect of age and gender by applying scale-wise covariance analysis (AnCova). Specifically, age and gender were treated as categorical factors and ^OS_σ(2, s), ^RS_σ(2, s), and r_σ(s) as covariates. Before AnCova, homogeneity of slopes model was evaluated for each s, significant result of this test means that AnCova is inapplicable due to an interaction between categorical predictors and covariates. In these cases the separate slopes model was used which include these interaction terms.

Software

For a more detailed description of our analytical flowchart as a guide for the FMF-SSC analysis, see (Eke et al., 2000, 2012; Mukli et al., 2015; Nagy et al., 2017). The above aspects of multifractal analyses of rsNIRS time series have been implemented in Matlab 2012 (The MathWorks, Inc., Natick, MA, U.S.A.) by custom scripts written by the authors based on the recently published “MultiFracTool” software (Mukli et al., 2015; Nagy et al., 2017; Racz et al., 2018). The toolbox can be requested from the corresponding author. Statistical analyses were performed with StatSoft Statistica 13.2.

Results

The Presence of True Bimodal Multifractality

All measured signals showed an apparent bimodal scaling function. True multifractality was confirmed in 44 subjects, two of them with an unacceptable fit of the bimodal model. For further details, see Appendix in Supplementary Material.

Impact of Age-and Gender on Multifractal Endpoint Parameters

In case of raw HbT signals, the degree of autocorrelation for the slow component, marked by ^sĤ(2), significantly increased with age unlike for the fast component (Figure 3A). Conversely, the neural component obtained with CBSI significantly decreased with age for the fast, but not for the slow component (Figure 3B). In the elderly group, the h_max of the slow component of the raw signal, ^sh_max, was found increased (Figure 3C), while it did not change with the CBSI-treated signal. These changes were the opposite—similarly to Hurst exponent—regarding ^fh_max in the young group (Figure 3D). Regarding the foci, for the raw HbT signals they were statistically the same in both age groups (Figure 3E), while for the fast component of the CBSI-enhanced HbT signal they were lower in elderly subjects (Figure 3F). Given that the (SD(HbO)/SD(HbR) ratios—the key element in CBSI model—did not differ between young and elderly groups (p = 0.543), the above significant differences should be regarded as real.

FIGURE 3

Figure 3. Results of univariate statistical analysis of multifractal endpoint parameters. Ĥ(2) of slow and fast components in the two age groups calculated from raw (A) and on CBSI-pretreated (B) HbT signals. h_max obtained from raw (C) and CBSI-treated (D) HbT signals. Focus of raw (E) and CBSI-pretreated (F) HbT signals. Recall that the raw HbT signal represents neuronal and non-neuronal events combined, while the CBSI-pretreated signal is a enhanced representation of the underlying neurodynamics. Accordingly, the fact that we found significant differences in the slow component for the raw HbT signal and in the fast component for the CBSI-pretreated signals identifies the slow component emerging dominantly from vascular events, while the fast component attributed mainly to neurodynamics.

Age-related differences remained significant when the slow and fast components were compared. Specifically, component-wise contrast – defined as ln(^sŜ_σ(N))/ln(^fŜ_σ(N))—turned out to be significantly different between the age groups (p = 0.03). A concomitant—albeit non-significant—decrease characterizing the fast component (p = 0.405) for focus contributed to an overall increased ratio of foci (p = 0.067).

Comparing the multifractal parameters of cerebral hemodynamic fluctuations of female and male subjects, the only significant difference was observed for their HbO slow component. Incorporating age groups rendered the gender-related differences non-significant.

In order to prove that the significant differences in endpoint parameters seen in Figure 3 attribute to alterations in a single subject, endpoint parameters related to slow and fast component were statistically evaluated in combination as dependent variables in multivariate analysis. When scale-free endpoint parameters of the same kind [Ĥ(2), h_max] were combined, MANOVA revealed a strong correlation (r² > 0.7). Furthermore, significant multifractal endpoint parameters of slow component of raw HbT and fast component of CBSI HbT showed joint significance in a multivariate analysis (p = 0.045, Wilk's test) suggesting the coincident change of both components. Taken it together, these findings indicate that the observed alterations in the multifractal endpoint parameters resulted from subject-wise aging. As to the key geometrical parameters of the multifractal scaling functions, p-values obtained from Wilk's test indicated an overall non-significant influence of age. Of note the two main estimates of the analysis – Ĥ(2) and ln(Ŝ_σ(N) – showed positive correlation for the fast component of the CBSI-pretreated (r² = 0.46) and the slow component of the raw HbT signal (r² = 0.58). The p-values of the statistical analyses are summarized Table 1.

TABLE 1

Table 1. Significance of gender-related differences (p-values).

Influence of Age and Gender on the HbO-HbR Relationship

Scale-Wise Fractal Cross-Correlation

The mean fractal scale-wise cross-correlation coefficient was found higher at all scales in the elderly group (Figure 4A). Importantly, this decrease in $\overline{r_{\sigma }(s)}$was more pronounced in young participants at higher values of s achieving significance (two-way ANOVA, confirmed by Tukey's post-hoc test) at a cluster of the corresponding high temporal scales (see the probability profile on Figure 4A). Based on this parameter calculated at these large scales, the oxy- and deoxyhemoglobin fluctuations are uncorrelated (random) among the elderly and were found anticorrelated in the young group. MANOVA revealed a statistical coincidence between the age-related increase in r_σ(s) at specific high scales (corresponding to 2,172, 2,414, 2,683, and 2,982 s) and the same alteration observed for each multifractal endpoint parameters [^sh_max, ^sĤ(2) of the raw HbT and ^fh_max, ^fĤ(2) of CBSI-treated HbT]. In addition, concerning r_σ(s) and ^TS_σ(2, s) as dependent variables their close relationship specifically appears at some of these aforementioned temporal scales. Since r_σ(s) is determined by the dynamics of oxygen delivery and its extraction in the brain (i.e., supply and demand), these are key results for discussing the impact of aging on cerebral oxygenation.

FIGURE 4

Figure 4. Age-related differences revealed by scale-wise and scale-free analysis of oxy- and deoxyhemoglobin relationship. (A) Fractal correlation coefficient of HbO and HbR as functions of scale in the elderly (upper) and young (middle) groups with the corresponding p-values (lower). At s_min, these were above 0 and decreased gradually toward-1 as s approached s_max in both age groups. (B) Multifractal covariance of ^fh_max (upper) and ln(^fŜ_σ(N); lower) of the fast component in the young and elderly groups.

Multifractal Covariance

In contrast to scale-wise fractal correlation analysis when covariance is normalized by σ(s), the multifractal covariance analysis allows for a moment-wise characterization of the scaling properties of HbO-HbR coupling extending also for q ≠2. This method revealed a statistically significant age-related difference concerning the fast component [see ^fh_max and ln(^fŜ_Cov(N)), Figure 4B], which is markedly correlated with ^fh_max (r² = 0.55) and ln(^fŜ_Cov(N)) (r² = 0.69) obtained for CBSI-pretreated HbT signals across age groups.

Power-law scaling of the multifractal HbO-HbR covariance function may either originate from the independent scale-free variance profiles of HbO and HbR or as from the coupled fluctuations of the two. In case of the fast component, the significant contribution from the latter is confirmed for the whole study population, given that ^ORH(q) derived from scale-wise covariances differed from (^OH(q)+^RH(q))/2 both obtained at q = 2 (p = 0.003, Wilcoxon matched pairs test). Moreover this pattern was absent in the elderly group (p = 0.116), but not in the young group (p = 0.006).

The correlation (r²) between dependent variables (i.e., multifractal endpoint parameters) captures the percentage of mutual variance of multifractal HbO-HbR covariance profiles and those obtained by SSC for the variance of CBSI-pretreated HbT signals. The percentage of mutual variance was the highest for ln(^fŜ_Cov(N)) (r² = 0.81), and there was a strong relationship between their focus ratios (r² = 0.72). However, the correlation was moderate for ^fh_max (r² = 0.56). To sum it up, results of multifractal covariance analysis seems rather consistent with an altered fast component of the CBSI-pretreated HbT signals.

Significance of Fractal Scale-Wise Cross-Correlation

For the sake of comparison, we calculated fractal scale-wise cross-correlation and averaged running Pearson and Spearman correlation coefficients at the same time scales (Figure 5).

FIGURE 5

Figure 5. Comparing means of different correlation coefficients between HbO and HbR calculated for all subjects. We explain the less correlated values observed for r_σ(s) with the effect of calculation order on trends in the signal. In case of the mean Pearson- and Spearman-coefficients, detrending and averaging step come after calculating the correlation effectively neglecting the fundamentally scale-free character of the processes. However, considering fractal scale-wise correlation coefficient, this latter step is the last preceded by averaging of bridge-detrended variances at a given scale. Please also note the characteristic transient—indicated by red arrow—only appearing in r_σ(s).

Variability of the HbT scaling function profiles at q = 2 and for all scales was explained as related to the independent variables–r_σ(s), ^OS_σ(2, s), and ^RS_σ(2, s)—based on the Bienaymé-formula given in Equation 2 (Figure 6 top). First, a series of multiple regression tests—not yet accounting for the effect of age and gender—were performed across all temporal scales that yielded positive correlation between each regressors and ^TS_σ(2, s). Importantly, the standardized regression coefficients proved to be significant for r_σ(s) in case of all temporal scales. The highest estimated effects were observed in case of ^OS_σ(2, s) while this test revealed the weakest effect of ^RS_σ(2, s) being a non-significant regressor of ^TS_σ(2, s) at high scales (note the empty blue circles on Figure 6 bottom).

FIGURE 6

Figure 6. Significance of multiple regression tests. (A) Scaling function value of HbT as a function of ^OS(q, s) (left), ^RS(q, s) (middle) and scale-wise fractal cross-correlation coefficient (right) acquired at equal scale and moment (exemplary case at q = 2 and s = 4828). (B) Regression coefficients and their significance (closed circle for cases of p < 0.05) related to scaling function values of HbO, HbR; and r_σ as functions of scale s. The estimated effect were significant for all s values in case of ^OS(q and r_σ(s), and the lower standardized regression coefficients indicated the weak (not significant for all scales) effect of ^RS(q, s). For further details, see main text.

Subsequently gender and age were incorporated in the multiple regression model as categorical predictors to evaluate their influence on the p-value of correlation between the covariates and variance profiles related to CBV fluctuations (^TS_σ(2, s)). In the GLM framework the appropriate choice was AnCova or separate slopes model depending on the prerequisites of each approach. Age in of itself turned out to be not determinant of scale-wise hemodynamic fluctuations captured by ^TS_σ(2, s). It is the scale-wise cross-correlation on the basis of which the impact of aging could be explained (Table 2). These results confirm the outcome of multiple regression analysis showing a significant and strong relationship between ^TS_σ(2, s) and ^OS_σ(2, s) and a less steep but still significant relationship between the dependent variable and the fractal scale-wise cross-correlation coefficient.

TABLE 2

Table 2. Homogeneity of Slopes Model/Separate Slopes Model/Covariance analysis results.

Discussion

In this study, we found that hemodynamics in the human brain cortex captured by NIRS technology in most of the cases exhibited a bimodal multifractal scaling emerging from a range of low and high temporal scales (slow and fast components, respectively). We suggest relating the slow component of the raw HbT signal to dominantly vascular (vasogenic, i.e., non-neural) dynamics, while we consider the fast component primarily resulting from neurovascular coupling (i.e., neurogenic component). In order to demonstrate the impact of healthy aging, CBSI pretreatment of the raw HbT signal was necessary to enhance the neurovascular contribution in the fast component. Our main result is two-fold: first, we demonstrate that the vasogenic hemodynamics (CBV fluctuations proportional to HbT concentration changes) show increased long-range autocorrelation in the elderly group compared to the young group which is in agreement with what we had found previously applying monofractal analysis (^lowPSD_w,e method⁵) within comparable scaling ranges (Eke et al., 2006). Second, we show that the fluctuations of the neural component are less correlated in the elderly group. This opposite influence of healthy aging on the slow vasogenic and the fast neurogenic fluctuations is consistent with an attenuated NVC. In support of this notion, we evidence that age-dependent alterations in HbO-HbR relationship is a manifestation of altered neurovascular coupling and also a determinant of scaling properties of CBV dynamics. Specifically, in in silico experiments we substantiated that (i) a decreased correlation in neurogenic fluctuations is attributed to decreased incoming signaling, (ii) HbO-HbR relationship became more correlated due to aging either as a result of decreased incoming signaling concomitant to lesser hemodynamic response or from increased vascular stiffening. These alterations do indicate that linear CBV dynamics is susceptible to aging. In contrast, non-linear CBV dynamics is spared by aging as demonstrated in unaltered degree of multifractality.

Multifractality of cerebral hemodynamics has been investigated extensively in case of blood oxygen level dependent (BOLD) signals using functional magnetic resonance imaging (fMRI). The influence of brain activity was demonstrated in the pioneering work of Shimizu et al. (2004). Later, findings on the topology of multifractal parameters and methodological refinements were reported (Wink et al., 2008; Ciuciu et al., 2012). On data from a publicly available imaging repository (Biswal et al., 2010), the effect of age and gender on multifractal spectrum was shown for resting-state fMRI-BOLD signals (Ni et al., 2014). While rsNIRS signals have been made subject to multifractal analysis in an earlier feasibility study (Dzung, 2010; Quang Dang Khoa and Van Toi, 2012), our study—to the best of our knowledge—is the first reporting on it elucidating some of the key underlying mechanisms using a scrutinized dataset with proven true multifractality.

Multifractal CBV Dynamics

CBV dynamics in the human brain cortex was shown to follow a complex, scale-free temporal pattern in the frequency domain that could be captured in the 1/f ^β model, where β is spectral index (Eke et al., 2006). The variance profile at q = 2 is analogous to the power spectrum, thus their apparent similarity can be readily shown (Figure 7). The estimated β obtained by ^lowPSD_w,e method is directly related to Ĥ(2), given the explicit relation of H(2) = (β+1)/2 (Eke et al., 2000). In fact, the power spectrum is equivalent to the Fourier transform of the signal's autocorrelation function according to Wiener-Khinchin theorem (He et al., 2010)⁶. Since its decay follows a power law with 2Ĥ(2)−1 as its exponent, Ĥ(2) and $\widehat{\beta }$ are interpreted as equivalent measures of LRC in the fractional Gaussian noise / fractional Brownian motion framework (Eke et al., 2002).

FIGURE 7

Figure 7. Relationship between power spectrum and scaling function profiles obtained by the FMF-SSC method. Although multifractal analyses were performed in the time domain, its design and results can be interpreted in the frequency domain as well owing to the explicit relationship between representation of dynamics in the temporal and the spectral domains both for the raw (A) and CBSI-pretreated (B) signal. The ordinate shows both frequency and temporal scale. It is the variance profile at q = 2 that corresponds well with the mean spectral estimates. The range of spectral estimates fall in between S(q, s) profiles for q = 15 and q = −15. The fractal scaling range emerges from the (non-fractal) noise (of biological origin) dominating the low temporal scales. It contains a breakpoint, which separates S(2, s) into a slow (associated with high temporal scales) and a fast (associated with intermediate temporal scales) scale-free component. This appears as a low- and a very-low frequency spectral band in the frequency domain. For further details, see text. Component-based focus is indicated by asterisk (*) as ln(Ŝ_σ(N)) for the raw (A) and ln(Ŝ_σ(N)) for the CBSI-pretreated HbT signal (B).

It should be recalled that Fourier transform builds on independency of frequency components. As multifractality can emerge from interactions between multiple time scales (Ihlen and Vereijken, 2010), it cannot be captured in the power spectrum alone. Nevertheless, its presence still could be detected in the form of phase-amplitude coupling [nested frequency, see Ref. (He et al., 2010)]. Thus capturing multifractality in the time domain can reveal the underlying multiplicative interactions between the temporal scales of the observation with ΔH₁₅ and fwhm as the measures of these cross-scale interactions (Ihlen and Vereijken, 2010).

Separation of Neurogenic and Vasogenic Multifractal Dynamics: CBSI-Pretreatment

In pursuit of the physiological origin of hemodynamic fluctuations, the analyses were performed both on raw and CBSI-pretreated (Cui et al., 2010) data. In addition, we carried out in silico experiments to substantiate the need of this preprocessing step to identify components dominated by vasogenic and neurogenic influences, respectively (see Supplementary Material). CBSI method builds on the assumption that maximally correlated fluctuations of HbO and HbR are not related to neural activity. Given our recent demonstration of CBSI pretreatment enhancing the neuronal component in the signal (Racz et al., 2017) our two-teared approach of signal processing allows for a distinction between influences of neuronal and non-neuronal nature in this study; aspects of particular interest in the aging process. Accordingly, the age-related differences of the calculated multifractal measures revealed for the fast component of the pretreated signal should reflect altered neurogenic fluctuations. Vice versa, we found that the significant age-related differences in the multifractal indices obtained for the slow component of raw NIRS signals disappeared after applying CBSI. This indicated that the slow component was non-neural, referred to as vasogenic.

Some authors pointed out specific limitations of CBSI to isolate the neural component in the signal in fNIRS studies using various stimulus response paradigms (Scholkmann et al., 2014). Although the neural activity readily and always induces anticorrelated dynamics in HbO and HbR—as postulated in CBSI –, the response to a specific task may elicit systemic confounding effects, too. In this regard, CBSI cannot be considered to be immune to global effects (Tachtsidis et al., 2008). Nevertheless, as we present results based on the analysis of resting-state NIRS observations we do not need to deal with such confounding influences in task.

In principle, we could not a priori exclude that CBSI-pretreatment would not distort the signal. The formulation of the pretreatment algorithm allows for analytical considerations about the influence of CBSI on multifractal parameters, which in cases of other pretreatment algorithms would be more difficult to make. In particular, the above step in the algorithm can be shown to affect only the scale-dependent measures such as ln(Ŝ_σ(N)) and breakpoint scales but not the scale-free parameters (Ĥ(2), h_max, ΔH₁₅, and fwhm), unlike with various filtering methods (Valencia et al., 2008). Nevertheless, these scale-dependent influences spared the impact of age on scale-free parameters because mean ^Oσ(N)/^Rσ(N) were found very similar in all measurement groups.

Origin of Multifractality in Resting-State Hemodynamic Fluctuations

Fluctuations of rsNIRS signal related to neural activity should be viewed as a sampled representation of an interim stage from intrinsic signal generation throughout the brain to the region of interest (ROI), where it is transformed into hemodynamic fluctuations. The regionally recorded rsNIRS signal—aside from systemic influences—is produced by the NVC driven by incoming signaling. Directly, it represents the hemodynamics within a population of vessels behaving like viscoelastic balloons in the ROI. As to the non-neural component of the rsNIRS signal, a likely origin of scale-free behavior is the numerous weekly coupled vascular source (diameter-dependent segmental oscillations along the arterial tree) blending into a fractally correlated pattern (Colantuoni et al., 1994).

Signal generation also raises questions about the spatial dynamics and resting-state functional connectivity. Regarding the incoming neural activity, electrocorticography records captured across various locations in the brain cortex have been shown scale-free temporal structuring (He et al., 2010; He, 2011). Moreover, the power spectral density of cortical EEG exhibits scale-free structuring not only in the temporal but—in an interrelated manner—in the spatial domain, too, (Freeman et al., 2003). Specifically, as demonstrated by these authors, fluctuations spanning from high-frequency/low-power bands to the low-frequency/high-power ranges reflect upon neural events propagating across the micro-meso-macro scales representing contributions from ion channels, across gyri all the way to those of lobes, respectively. Hence a PSD and scaling function representations of neurodynamics and coupled hemodynamics could be viewed as capturing the information flow within the system from its sources via inhomogeneous network routes eventually converging onto the signaling input of the monitored ROI (Buzsaki, 2006). The sampled representation of this process will typically show inhomogeneously distributed fluctuations, visible as intermittent periods of small and large variability; genuine properties of multifractal processes (Ihlen and Vereijken, 2010).

When the multiplicative cascading process was extended into the spatial domain, a description was obtained comparable to the one by the self-organizing branching process (Zapperi et al., 1995); a refined extension of the SOC model (Bak et al., 1987). The inference is that the intermittent in essence multifractal temporal patterns and the inhomogeneous incoming network connections are manifestations of the same phenomenon: emergence of intermittent regional activity from multiple sites of the brain converging via multiplicative interactions between spatiotemporal scales as integrated incoming signaling in the ROI.

Indeed, our findings related to fractal dynamics could potentially reflect the presence of SOC (Bak et al., 1987) in the observed physiological subsystems shaping cerebral hemodynamics in the ROI. SOC, substantiating the 1/f noise-type neurodynamics of the human brain, builds on the notion that the brain dissipates the local low-frequency perturbation elicited by external or internal stimuli without any particular spatial or temporal scales (Stam, 2005; Bullmore et al., 2009; Chialvo, 2010; Sporns, 2011). When interactions between scales occur, multifractality can readily emerge in systems showing properties of SOC (Tebaldi et al., 1999; Lima et al., 2017). Of relevance, in a recent rsNIRS study using 16 channels sampling of resting-state hemodynamics in the PFC, it has been evidenced that the presence of critical state in resting-state dynamic functional connectivity (Racz et al., 2018). In sum, the measured signals are considered as a composite of hemodynamic fluctuations of vasogenic and those elicited by incoming signaling of neurogenic origins with NVC as the link between the two.

Interpretation of Multifractal Endpoint Parameters

The vasogenic component of the rsNIRS signal in terms of its observed multifractal temporal patterns can be interpreted as a consequence of attenuated neurovascular coupling—reflected by an anticorrelated → random shift in the fractal cross-correlation of HbO and HbR. As to the neurogenic fluctuations here we explain the altered multifractal endpoint parameters resulting from interactions between multiple time scales along functional connections (Ihlen and Vereijken, 2010). Along with our focus-based multifractal formalism (Mukli et al., 2015) and a small world implementation (Watts and Strogatz, 1998) of the concept of self-organized criticality (Mandelbrot, 1974; Ihlen and Vereijken, 2010) offer a concise framework for the interpretation of the results obtained in this study the way outlined in the followings. Further details are provided below and please also see the Supplementary Material.

Linear dynamics: H(2) vs. h_max

Since the degree of global LRC is quantified by H(2), its changes can be interpreted as increased or decreased persistence, meaning correlation of a non-stationary process (Eke et al., 2000; Herman et al., 2011). Furthermore, Deligniéres et al. established a relationship between network degeneracy—meaning partial overlap in heterogeneous functional connections—and output signal correlation (Delignières and Marmelat, 2013). Accordingly, the Hurst exponent does not only reflect upon global scale-free properties emerging from a network, but its degree of degeneracy, too.

While LRC—and thus H(2)—reflects global scale-free properties, the Hölder trajectory is a local scale-free measure varying along the signal. Although multifractal analysis of physiological data usually shows tightly correlating changes of h_max and Ĥ(2), the interpretation of h_max as a measure of correlation within the signal is only approximate, since it is associated with q = 0, not q = 2.

Non-linear dynamics: Δh₁₅ and fwhm

Though ΔH₁₅ is defined on H(q) and fwhm is derived from D(h), their excellent correlation – owing to the deterministic formalism established by Equations 3–5 – offers a rationale to interpret them together. The applications of q-order statistics reveals non-linear properties in scaling of the signal (Ashkenazy et al., 2003). Thus these two multifractal endpoints are indeed equivalent measures of multiplicative interaction between temporal scales of the observed dynamics process (Ihlen and Vereijken, 2010). Importantly, ΔH₁₅ and fwhm should be regarded as indicators of non-linear dynamics (Gómez-Extremera et al., 2016; Bernaola-Galván et al., 2017).

Asymmetry of D(h)—an occasionally observed phenomenon—could be incorporated in the analysis of multifractality in terms of W = W+/W– where fwhm is equal to the sum of W+ and W– (Wink et al., 2008), corresponding to the width of left and right-half of the singularity spectrum. We calculated the W and found the shape of our singularity spectra symmetric and not affected by age and gender. We stress that testing for true multifractality based on ΔH₁₅ statistics (for details see Appendix in Supplementary Material) is an essential prerequisite when it comes to evaluating changes in these endpoint parameters in regard of multifractal CBV dynamics.

A scale-dependent measure of hemodynamic power: focus

The focus of the scaling function is a key element of our regression scheme for obtaining H(q) thus securing a robust estimate of D(h) free of inversion (Mukli et al., 2015). Importantly, it is also a robust scale-dependent statistics estimated at signal length as a point of convergence for the scaling function profiles. Since our analysis is based on the SSC method using bridge-detrended variance, Ŝ_σ(N) is essentially the variance associated with the whole signal. Given that the coefficient of variance for our rsNIRS signals were the same, Ŝ_σ(N) is also the measure of the signal mean, which is consistent with our SOC-simulations (Figure S3). In the frequency domain, it is analogous with the power of the DC-component of the signal (see the behavior of spectra and scaling functions on Figure 7 as f → 0).

Since its value is influenced by numerous other variables, conclusions regarding hemodynamic alterations can be drawn if focus is interpreted together with Ĥ(2). Given their analogous frequency domain parameter the hemodynamic power corresponding to a spectral band can be estimated. The separated components of the rsNIRS signal have distinct scaling ranges, and the area under the scaling function corresponding to such SRs is an approximation of summed logarithmic variance of the given component. This area can be explicitly calculated as it is proportional with ln(Ŝ_σ(N)) and width of SR, and it also increases with decreased Ĥ(2). Given the straightforward relationship between summed variance and total hemodynamic power in a given temporal/frequency range, the obtained results for the area should be viewed as power of hemodynamic fluctuations associated with the isolated components.

The Inference of Bimodality

The majority of the measured rsNIRS signal showed bimodal scaling that was statistically confirmed by comparing the errors of fits for the bimodal and unimodal models. We used a scaling-range adaptive method to assess scaling exponents and multifractal endpoint parameters of the two components, which approach has already been used in our previous study (Nagy et al., 2017) and in other studies as well (Ge and Leung, 2012; Kuznetsov et al., 2013). The apparent structural heterogeneity in the scaling functions of our dataset (convex and/or concave transient range) prompted us to choose the robust moment-wise SR-adaptive method instead of an alternative (decomposition of scaling function) that was specifically designed to assess additive properties of bimodal scaling.

The bimodal analysis is an adaptive tool that separates two fractal SRs within an overall range of scales between the lowest temporal scale of 8 seconds (s_min = 16) and 4096 seconds (s_max = 8192). Although the Fourier transform does not assume an exact relationship between temporal scales and frequencies, we can still assign the non-fractal “noise” component below s_min to a spectral range from 0.125 to 1 Hz (i.e., Nyquist frequency). It is known that this band is dominated by fluctuations of systemic origin such as cardiac pulsation and respiration (Tian et al., 2009; Sassaroli et al., 2012), which is corroboratory regarding the exclusion of time scales shorter than 8 s.

Earlier we introduced the slow and fast components as dominantly neurogenic and vasogenic, in the followings the reader is taken through arguments substantiating this division. The fast and slow components were identified in the time domain within their respective analytical scaling ranges in the time domain. In a frequency domain representation, they correspond to the low- and very low-frequency oscillations (LFO and VLFO, respectively; Figure 7). The LFO/VLFO classification was used by authors evaluating fundamental spectral aspects of NIRS dynamics (Obrig et al., 2000; Schroeter et al., 2004; Li et al., 2013; Vermeij et al., 2014). LFO is generally regarded of neural origin and thus is a commonly investigated in NIRS studies of cognition (Chance et al., 1993) and functional connectivity (Sasai et al., 2011). Also, it substantiated the concept of fMRI-based functional connectivity ever since Biswal et al. (1995) analyzed cross-correlations in paired signals by clipping their power spectra to zero above 0.1 Hz in order to identify temporal coincidences in local activities thus excluding the global influences. It seems reasonable to separate the VLFO component since its dynamics is dominated by non-neural (Schroeter et al., 2004; Li et al., 2013; Vermeij et al., 2014), particularly endothelium-related mechanisms (Li et al., 2013; Chen et al., 2014). VLFO would manifest as the slow, vasogenic component in our study within the breakpoints and s_max, while the range down to s_min corresponds to LFO (Figure 7). Nevertheless, the debate is still ongoing over further contributors to CBV dynamics such as vasomotion (Elwell et al., 1999) and Mayer-waves (Sassaroli et al., 2012); effectively being excluded from our analyses by setting s_min to 8 s corresponding to 0.125 Hz. In fact, vasomotion may show up, but within a narrow range of scales and at low fluctuation amplitudes thus having a weak effect if any on our scaling analysis. Please note that the assessment of bimodality and scale-free properties of both LFO and VLFO were only possible because the analyses were performed with a high-enough s_max as recommended by Nagy et al. (2017). The consequent statistical instability is compensated by our focus-based regression model.

Healthy Aging Is Associated With Altered Complexity of Cerebral Hemodynamics

In the present study, decreased Ĥ(2), h_max, and ln(Ŝ_σ(N)) were seen in the fast—neurogenic—component of the CBSI-pretreated NIRS signals of the elderly participants. In addition, we found an increased Ĥ(2) and h_max in the elderly group for the slow—vasogenic—component of the raw rsNIRS signal (Figure 3) with spared ln(Ŝ_σ(N)) and multifractality (i.e., no difference in ΔH₁₅ and fwhm). These changes in the multifractal endpoint parameters are consistent with attenuated neurovascular coupling concomitant to declining incoming signaling and impaired vascular responses.

Altered Neurogenic Component Due to Declining Neurodynamics

Multifractal measures of the fast component (~ LFO) of the CBSI-pretreated NIRS signal revealed a difference between young and elderly subjects. Specifically, Ĥ(2), h_max, and ln(^fŜ_σ(N))—a key parameter characterizing the overall decline in the neurogenic component—was found decreased among the elderly participants. In principle, a decrease of Ĥ(2) and h_max could have resulted from the contribution of biological noise. The decreased ln(^fŜ_σ(N)) can be directly interpreted as a decreased power in LF oscillations which was concluded by other studies, too (Schroeter et al., 2004; Li et al., 2013; Vermeij et al., 2014). Although the decline in power appears at all frequencies, a disproportionate decrease could still manifest from biological noise present across the lower frequencies first seen in our previous study (Eke et al., 2006). A recent numerical study clearly showed that varying signal/noise ratio by adding white noise yielded underestimation of H(q) (Ludescher et al., 2011). Hence our decreased Ĥ(q) and h_max of the fast component can at least in part be explained by the relative impact of biological noise progressively dominating the higher temporal scales. Although our choice of s_min exceeds this range of temporal scales, this factor must be taken into account in our interpretation of the observed alterations of the neurogenic component associated with the LFO range in our rsNIRS signals. In spite of the spurious estimates of Ĥ(q) attributable to the increased relative impact of biological noise, the conclusion of an overall decline in the neurogenic component can be still held. Nevertheless, please note that CBSI-pretreatment of raw HbT signal effectively removes the impact of non-anticorrelated dynamics in the signal. Therefore the impact of biological noise is unlikely and the observed changes should be regarded as real.

This allows us to interpret the decreased focus as a sign of decreased incoming signaling in line with the dominantly neurogenic origin of the LFO and the aforementioned in silico observations (Supplementary Material). As multifractality can be viewed as resulting from cross-scale spatiotemporal interactions (Monto, 2012), unchanged ΔH₁₅ and fwhm suggest that healthy aging spared these interactions within the incoming networks.

As to the multiple mechanisms involved in neurogenic aspects of the aging process, based on our results we argue that even healthy aging leads to progressive attenuation in incoming signaling that to some unknown extent could be masked by impaired neurovascular coupling. The aging process is known leading to gray matter atrophy associated with dropping neuronal count (McGeer et al., 1984) and lower gray matter density (Sowell et al., 2003) along with impaired synaptic activity. The latter is known to be prevailing in postmenopausal women due to lower levels of estrogen (Gibbs and Aggarwal, 1998; Khan et al., 2013). In addition, an age-related decrease in the hemodynamic response elicited by cognitive task has been observed in the human prefrontal cortex by NIRS (Schroeter et al., 2004). Several studies have evidenced that significant changes in resting state functional connectivity take place in the aging brain at large and small spatial scales alike with inference to temporal dynamics (Ferreira and Busatto, 2013). For example, the structural changes occurring in the aging brain imply changes in its functional connectivity readily manifesting in altered parameters of complexity parameters (Sun et al., 2012). Specifically, global deleterious effects of older age has been reported on functional networks (Achard and Bullmore, 2007) mounting to “topological marginalization” like that of the prefrontal cortex due to segregation in the global network. It was suggested that even healthy aging disrupts the underlying networks by severing the long-range connections, especially in higher-function areas like the prefrontal cortex (Ferreira and Busatto, 2013). In a sample of subjects (19–80 years) a linear effect of age associated with impaired resting-state functional connectivity has been demonstrated (Mevel et al., 2013).

Indeed, when eliminating the multiplicative neural interactions across spatial scales in an interaction-dominant model of human cognition (Ihlen and Vereijken, 2010), more aggregated, focal activities remain. Putatively, this model behavior (global → focal shift in intrinsic activities of the brain) is capable of explaining the loss of complexity in our rsNIRS signal as the sign of decreased incoming signaling. Nevertheless, this takes place in the PFC during healthy aging to such an extent that would still not interfere with the cross-scale spatiotemporal interactions in the observed dynamics captured in unaltered parameters of multifractality.

A further theory suggests the role of dedifferentiation in the aging of physiological subsystems like the brain function (Sleimen-Malkoun et al., 2014). Experiment using a motor task paradigm indeed demonstrated that in addition to stronger activation in dedicated regions to a motor task, older adults generally exhibit activation of additional areas of the brain not or only marginally involved in young participants (Sleimen-Malkoun et al., 2014). Approaches relying on fMRI-based connectivity studies accounting for the difference in task and connectivity paradigms demonstrated that higher levels of activity coexists with disrupted connectivity (Sala-Llonch et al., 2015).

Altered Vasogenic Component Due to Impaired Vascular Responses

The observations regarding the vasogenic component are compatible with those of our previous monofractal study using raw rsNIRS signals reporting on an increased spectral index β in the elderly suggesting the impact of age-related vascular sclerosis on CBV dynamics (Eke et al., 2006). The increased Ĥ(2) and focus indicate a more correlated vascular dynamics with decreased hemodynamic power in the elderly group. The mean singularity spectrum for the vasogenic component was found shifted to the right (reflected by h_max) with maintained shape and width (Figure 8). Consistent with this view, we hypothesize that healthy aging leads to this increased correlation pattern in the vasogenic (VLF) component as a result of vessel stiffening and a decline in the endothelium-mediated (metabolic) regulation of cerebrovascular smooth muscle tone. Our evidence to this hypothesis is indirect, though: the result of CBSI-pretreatment and the temporal scales characterizing VLF hemodynamics. Nevertheless, this interpretation of our results is in agreement with several studies on this component of cerebral hemodynamics. In an rsNIRS study investigating the effect of age, Li et al. demonstrated (Li et al., 2013) decreased average amplitudes of spontaneous oscillations in the elderly. The latter authors assumed that the oscillations in the 0.005–0.02 Hz range originated from the endothelium. Aging has also been shown decreasing the responsiveness of distal segments of the arterial tree due to endothelial dysfunction (Toda, 2012) and increasing wall stiffness (Schroeter et al., 2004; Zhu et al., 2011; Wardlaw et al., 2013; Vermeij et al., 2014). The functional hyperemia studied by selectively interrupting endothelial signaling in the somatosensory cortex of rats confirmed its key role in mediating the very slow—maintained—hemodynamic response brought about by NVC (Chen et al., 2014). The vascular decline in the elderly thus may attenuate the local hemodynamic response, too.

FIGURE 8

Figure 8. Input and output parameters of the focus-based multifractal formalism in the young and elderly groups for the slow (vasogenic) and fast (neurogenic) components of CBV fluctuations. Average scaling functions (A,B) and average singularity spectra (C,D) obtained for the raw and CBSI-pretreated HbT signal. Age coding: young—red (A,C), blue (C,D); elderly—green (A,C), orange (B,D). Preserved focus and increased correlation characterize the age-related changes of the vasogenic component (A). Conversely, these parameters of the neurogenic component indicate a decline apparent in the average scaling functions (B). The alteration of h_max is a good indicator of the right-ward shift of D(h) for the vasogenic (C) and a left-ward shift for the neurogenic component (D). It also reveals altered correlation properties if the shape and distribution of singular behavior is not much affected, which happens to be the case in this study as the distribution of Hölder exponents was found symmetrical around h_max irrespective of age and gender.

As to the attenuation in local vascular dynamics, indeed, the local hemodynamic response elicited by incoming neural activity is known to be driven locally by fast glutamate-mediated signaling, and more globally by amine- and ACh-mediated neural systems (Attwell and Iadecola, 2002). Both have been shown declining with age (McGeer et al., 1984; Gibbs and Aggarwal, 1998; Attwell and Iadecola, 2002). The above scenario may indicate that an impaired neurovascular coupling (Fabiani et al., 2014; Tarantini et al., 2017) could to some degree mask the effects of deterioration of connectivity on neurogenic signal complexity.

Implications of HbO-HbR Relationship

Age-related physiological dysregulation is essentially a gradual and typically irreversible loss of regulatory control originating from structural instabilities in regulatory systems (Cohen, 2016). We shall discuss this phenomenon concerning parameters reflecting neurovascular coupling with coupled HbO-HbR dynamics in its focus. We used multifaceted approach for its characterization: (i) scale-wise fractal cross-correlation, (ii) multifractal covariance and (iii) statistical analysis of Bienaymé-formulation of the HbO-HbR relationship. At this end, we interpret our results as the impact of age on the output of an integrated system of neurodynamics, coupled HbO-HbR dynamics and hemodynamics as genuinely interrelated aspects of neurovascular coupling.

Age-Related Increase in Scale-Wise Fractal Cross-Correlation

In elderly participants, the higher r_σ(s) capturing HbO-HbR relationship indicate the relatively larger contribution of correlating systemic hemodynamics to the recorded NIRS signal (which) is also consistent with declining neurodynamics. Furthermore, scale-wise correlation coefficients were found significantly elevated for rather high temporal scales directly influencing this domain of the variance profiles associated with the vasogenic component (Figure 4A). In a certain range of scales characterizing VLFO, the increase in r_σ(s) values were significantly associated with the ^sh_max and ^sĤ(2) of the raw HbT. Given the demonstrated association between endothelial-mediated responses and VLFO (Stefanovska et al., 1999), the pattern found in our study should be regarded as evidence supporting the endothelial contribution to the age-related increase in ^sh_max and ^sĤ(2).

Analyzing the coupled fluctuations of oxy- and deoxyhemoglobin has been used to assess cerebral oxygenation changes and the underlying processes (Reinhard et al., 2006; Wylie et al., 2009; Pierro et al., 2012). There is scarce evidence to determine the effects of aging on their relationship and the rationale behind the utilized parameters is also debatable. Please note the difference in the scale-wise pattern of r_σ(s) only showing a clear gradual decrease below ≈50–100 s (see the red arrow on Figure 5). It is reasonable to assume that such pattern reflects the relative contribution of mechanisms eliciting correlated or anticorrelated chromophore dynamics. Supposing the origin of the anticorrelated dynamics within the local balloons, and regarding the correlated dynamics mainly of systemic origin, this pattern may be informative of a relative impact of local and systemic hemodynamics on our measured rsNIRS signals. Specifically, local determinants of oxygen supply and extraction dominate the correlation across wide range of frequencies below 0.01 Hz (Stefanovska et al., 1999) emphasizing the potential contribution of spatiotemporally sustained response mediated by astrocyte-endothelial signaling (Iadecola and Nedergaard, 2007).

The Significance of Non-linear Relationship Revealed by Multifractal Covariance Analysis

While scale-wise cross-correlation analysis does not reveal scale-free properties, multifractal covariance analysis is specifically designed to characterize LRCs and multifractality in the coupled HbO-HbR dynamics. The covarying fluctuations of HbO and HbR originate from their respective individual fluctuations and the directly coupled dynamics of the exchange between the two compartments. Following the approach of Kristoufek (2011) the deviation of ^ORĤ(2) from (^OĤ(2)+^RĤ(2))/2 indicates oxygen exchange within the hemoglobin pool. We found a clearly significant decrease of ln(^fŜ_Cov(N)) in the elderly (Figure 4B) indicating weakening of the coupled fluctuation of HbO and HbR. Indeed, the inference is an upset balance between oxygen demand and supply. Given the strong association between the foci of the fast component of the covariance scaling function and of CBSI-pretreated signal scaling function, this change supposedly can be attributed to the age-related decline in neurodynamics. Although at a weaker significance and lower correlation with the corresponding multifractal endpoint parameter (obtained for multifractal covariance analysis and multifractal analysis of CBSI-pretreated signals, respectively), the decreased ^fh_max of HbO-HbR covariance in the elderly participants besides its similar pattern with ^fh_max of neurogenic fluctuations supports this notion. Clearly, more evidence is needed coming from measurements and synthesized datasets to elaborate and consolidate physiological models for this interpretation.

Although our approach is essentially similar to the multifractal detrended- (Zhou, 2008), height cross-correlation (Kristoufek, 2011) analyses and cross-wavelet analysis (Ciuciu et al., 2014; Jiang et al., 2016), this is the first paper describing and utilizing the bivariate adaptation of multifractal SSC analysis. Owing to the similarities between our implementation and other time domain algorithms a reliable parameter estimation was expected. A performance characterization of multifractal covariance analysis on synthetic signals with known degree of correlation, while is certainly desirable, but is beyond the scope of the present study.

The Importance of Coupled HbO-HbR Fluctuations Driving Cerebral Hemodynamics

Physiologically, oxygen consumption related to neural activity is not blood flow-limited, but rather neural activity controls the CMRO₂ and cerebral hemodynamics (Raichle and Mintun, 2006). Consequently, it is plausible to explain variations of CBV fluctuations (proportional to the measured HbT signal) as a function of HbO-HbR relationship providing insight in variations in cerebral oxygenation. The explicit deterministic, quantitative relationship between key variables of the study established by Bienaymé-formula (Equation 2) offers an explicit way assigned by GLM to evaluate their role with special emphasis on r_σ(s).

Scale-wise cross-correlation is a robustly significant determinant of HbT fluctuations at all temporal scales within the framework of GLM. Omitting either age and gender from the statistical analyses (i.e., multiple regression, Figure 6) or incorporating these categorical variables in it (AnCova/separate slopes) yields essentially the same outcome on the role of r_σ(s). Considering the impact of age on r_σ(s) and the age-corrected dependence of r_σ(s) CBV fluctuations (see Table 2) it can be concluded that age exerts its effect on ^TS_σ(2, s) in part via r_σ(s) shifting toward uncorrelated chromophore dynamics. The impact of weakened HbO-HbR coupling in the elderly—supposedly due to declining neurodynamics—is clearly seen in the age-related changes of multifractal endpoint parameters characterizing cerebral hemodynamics.

Limitations and Future Perspectives

Acquiring additional modalities would have allowed for a more straightforward physiological interpretation of the signal. Specifically, transcranial Doppler measurement (capturing blood flow velocity in the middle cerebral artery) or continuous blood pressure records offer measures on systemic influences that could have been regressed out. To compensate, we applied CBSI to remove these influences from the signal. Further development and testing of robust pre-processing methods are in place to further enhance the interpretation of multifractal measures of physiological processes.

The intimate relationship between multifractal hemodynamic fluctuations and functional connectivity has been demonstrated and characterized on intrinsic fMRI networks (Ciuciu et al., 2014) and by revealing scale-free network dynamics in the prefrontal cortex captured by NIRS (Racz et al., 2018). In this comparison, our single-region measurement measurement appears like a limitation (Novi et al., 2016). However, in fact, we present an analytical framework which is capable of integrating aspects on incoming signaling with those of regional hemodynamics elicited by neurovascular coupling.

The tools developed for characterizing the coupled HbO-HbR dynamics have potentials in other applications, too. The fractal scale-wise correlation analysis captures linear aspects while the multifractal covariance analysis adds a non-linear dimension to the assessment of coupling between non-stationary time series. These methods open new ways in uni- and multimodal applications to investigate functional connectivity or neurovascular coupling, respectively.

Conclusions

Mono- and multifractal approaches have greatly enriched our insight of biological complexity in particular that of the brain (Bullmore and Sporns, 2009; Herman et al., 2009; Ihlen and Vereijken, 2010; Nagy et al., 2017; Racz et al., 2018). To the best of our knowledge, this study is the first using consolidated datasets with tested and proven correlation-type multifractality for an in-depth characterization of resting-state NIRS fluctuations. Here we interpret the multifractality of single-region cerebral hemodynamics as resulting from neurogenic oscillations via cross-scale interactions blending into a scale-free intermittent arrhythmic pattern via neurovascular coupling. However, the intrinsic and endothelium-evoked heterogeneous oscillations of the vascular smooth muscle tone give rise to multifractality of vasogenic fluctuations. The multifractal endpoint parameters obtained from raw and pre-processed signal attest to the impact of healthy aging on cerebral hemodynamic fluctuations in the human prefrontal cortex. Specifically, we report that the total power of very-low frequency—vasogenic—oscillations of CBV decreased due to a preserved value of focus and an increase of LRC, (^sĤ(2)); the latter is concomitant to the right-shifting singularity spectra [D(h) along with its ^sh_max]. As to the hemodynamic fluctuations elicited by neural activity—related to low-frequency oscillations—we show a general decline indicated by decreased Ĥ(2), h_max, and focus of the neurogenic component. On the contrary, parameters reflecting degree of multifractality are the same in the group of young and elderly subjects which demonstrates non-pathological aging spares non-linear hemodynamics. In case of the elderly participants, the anticorrelation of HbO and HbR fluctuations were barely present at high temporal scales, while an attenuated cross-correlation was revealed by multifractal covariance analysis. We show that the impact of age on the parameters of neuro- and vasogenic components must have resulted from the age-related alterations in HbO-HbR coupling. In our study, the HbO-HbR relationship appears as a key element directly influenced by the neuronal activity and directly coupled to CBV dynamics via neurovascular coupling which seems like sensitive to aging. We suggest that decreased incoming signaling and the prevalence of an altered pattern of moment-to-moment HbO-HbR coupling could contribute to the mismatch between oxygen demand and supply. Together with vascular dysfunction could well be considered as factors behind the observed age-dependent alterations of cerebral hemodynamics.

Author Contributions

PM adapted the multifractal analytical framework to resting-state NIRS signals, performed the analysis on measured data and wrote the manuscript. ZN developed methods and helped writing the manuscript. FR implemented a testing framework and assessed true multifractality of the measured data. PH carried out the measurements. AE helped developing and writing the manuscript and provided conceptual guidance in the study.

Conflict of Interest Statement

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

Acknowledgments

The authors acknowledge the support from the Hungarian Scientific Research Found OTKA grant T 034122 awarded to AE.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01072/full#supplementary-material

Abbreviations

BOLD, blood oxygen level dependent; CBSI, correlation based signal improvement; ΔH₁₅, the difference between the H(−15) and H(15) values; ^f, measure describing the fast signal component dominating over the low-frequency region; FMF, focus-based multifractal formalism (an approach using focus-based regression scheme); fwhm, full-width of the singularity spectrum, D(h), at half of its maximum; h_max, maximal Hölder exponent at which singularity strength (D) is equal to 1; LRC, long-range correlation; N_s, number of non-overlapping segments; PSD, power spectral density; rsNIRS, resting-state near-infrared spectroscopy; r_σ(s), scale-wise fractal correlation coefficient; ^s, measure describing the slow signal component dominant over the lower scales; s', scaling boundary (possible breakpoint); ^XS_σ(q, s), scaling function value at a given moment order (q) and temporal scale (s), calculated from signal X, with σ as measure; ln(^XŜ(N)), the focus of the scaling function for signal X; SR, scaling range; SSC, signal summation conversion (method); SSE, sum of squared error; v, the order of non-overlapping segments v = 1, …, N_s; (V)LFO, (very) low-frequency oscillation.

Footnotes

1. ^q takes all integer values in this interval.

2. ^Different signal types (T – HbT; O – HbO; R – HbR; OR – relationship of oxy- and deoxyhemoglobin) are indicated by a left superscript to the variable of multifractal analysis.

3. ^Note that the Bienaymé formula is not specific to fractal time series.

4. ^Parameters obtianed directly from multifractal regression analysis are denoted as estimates with a hat: H(q) and ln(Ŝ_σ(N)).

5. ^low – right half (high frequencies) of the spectrum is excluded, w – windowing, e – endmatching.

6. ^Since Wiener-Khinchin theorem only applies to wide-sense stationary processes, the autocorrelation function of the increment process is considered according to the fGn6fBm framework.

References

Achard, S., and Bullmore, E. (2007). Efficiency and cost of economical brain functional networks. PLoS Comput Biol. 3:e17. doi: 10.1371/journal.pcbi.0030017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashkenazy, Y., Havlin, S., Ivanov, P. C., Peng, C. K., Schulte-Frohlinde, V., and Stanley, H. E. (2003). Magnitude and sign scaling in power-law correlated time series. Phys. A Statist. Mech. Appl. 323, 19–41. doi: 10.1016/S0378-4371(03)00008-6

CrossRef Full Text | Google Scholar

Attwell, D., and Iadecola, C. (2002). The neural basis of functional brain imaging signals. Trends Neurosci. 25, 621–625. doi: 10.1016/S0166-2236(02)02264-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: an explanation of the 1/f noise. Phys. Rev. Lett. 59, 381–384. doi: 10.1103/PhysRevLett.59.381

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabási, A. L., and Vicsek, T. (1991). Multifractality of self-affine fractals. Phys. Rev. A 44, 2730–2733. doi: 10.1103/PhysRevA.44.2730

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabási, A. L., Szépfalusy, P., and Vicsek, T. (1991). Multifractal spectra of multi-affine functions. Physica A 178, 17–28. doi: 10.1016/0378-4371(91)90072-K

CrossRef Full Text | Google Scholar

Bassingthwaighte, J. B., Liebovitch, L. S., and West, B. J. (1994). Fractal Physiology. New York, NY: Oxford University Press.

PubMed Abstract | Google Scholar

Beckers, F., Verheyden, B., and Aubert, A. E. (2006). Aging and nonlinear heart rate control in a healthy population. Am. J. Physiol. Heart Circ. Physiol. 290, H2560–H2570. doi: 10.1152/ajpheart.00903.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernaola-Galván, P. A., Gómez-Extremera, M., Romance, A. R., and Carpena, P. (2017). Correlations in magnitude series to assess nonlinearities: Application to multifractal models and heartbeat fluctuations. Phys. Rev. E 96:032218. doi: 10.1103/PhysRevE.96.032218

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswal, B. B., Mennes, M., Zuo, X. N., Gohel, S., Kelly, C., Smith, S. M., et al. (2010). Toward discovery science of human brain function. Proc. Natl. Acad. Sci. U.S.A. 107, 4734–4739. doi: 10.1073/pnas.0911855107

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar mri. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullmore, E., Barnes, A., Bassett, D. S., Fornito, A., Kitzbichler, M., Meunier, D., et al. (2009). Generic aspects of complexity in brain imaging data and other biological systems. NeuroImage 47, 1125–1134. doi: 10.1016/j.neuroimage.2009.05.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain. New York, NY: Oxford University Press.

Google Scholar

Caccia, D. C., Percival, D., Cannon, M. J., Raymond, G., and Bassingthwaighte, J. B. (1997). Analyzing exact fractal time series: evaluating dispersional analysis and rescaled range methods. Phys. A Statist. Mech. Appl. 246, 609–632. doi: 10.1016/S0378-4371(97)00363-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Chance, B. (1994). Current state of methodology on hemoglobin oximetry in tissues. Adv. Exp. Med. Biol. 345, 23–32.

PubMed Abstract | Google Scholar

Chance, B., Nioka, S., and Zhao, Z. (2007). A wearable brain imager. IEEE Eng. Med. Biol. Mag. 26, 30–37. doi: 10.1109/MEMB.2007.384093

PubMed Abstract | CrossRef Full Text | Google Scholar

Chance, B., Zhuang, Z., UnAh, C., Alter, C., and Lipton, L. (1993). Cognition-activated low-frequency modulation of light absorption in human brain. Proc. Natl. Acad. Sci. U.S.A. 90, 3770–3774. doi: 10.1073/pnas.90.8.3770

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B. R., Kozberg, M. G., Bouchard, M. B., Shaik, M. A., and Hillman, E. M. (2014). A critical role for the vascular endothelium in functional neurovascular coupling in the brain. J. Am. Heart Assoc. 3:e000787. doi: 10.1161/JAHA.114.000787

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/Nphys1803

CrossRef Full Text | Google Scholar

Ciuciu, P., Abry, P., and He, B. J. (2014). Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks. Neuroimage 95, 248–263. doi: 10.1016/j.neuroimage.2014.03.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciuciu, P., Varoquaux, G., Abry, P., Sadaghiani, S., and Kleinschmidt, A. (2012). Scale-free and multifractal time dynamics of fmri signals during rest and task. Front. Physiol. 3:186. doi: 10.3389/fphys.2012.00186

PubMed Abstract | CrossRef Full Text | Google Scholar

Clauset, A., Shalizi, C. R., and Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Rev. 51, 661–703. doi: 10.1137/070710111

CrossRef Full Text | Google Scholar

Cohen, A. A. (2016). Complex systems dynamics in aging: new evidence, continuing questions. Biogerontology 17, 205–220. doi: 10.1007/s10522-015-9584-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Colantuoni, A., Bertuglia, S., and Intaglietta, M. (1994). Microvascular vasomotion: origin of laser Doppler flux motion. Int. J. Microcirc. Clin. Exp. 14, 151–158.

PubMed Abstract | Google Scholar

Cope, M., Delpy, D. T., Reynolds, E. O., Wray, S., Wyatt, J., and van der Zee, P. (1988). Methods of quantitating cerebral near infrared spectroscopy data. Adv. Exp. Med. Biol. 222, 183–189. doi: 10.1007/978-1-4615-9510-6_21

PubMed Abstract | CrossRef Full Text | Google Scholar

Csermely, P. (2006). Weak Links: The universal key to the stability of networks and complex systems. Berlin; Heidelberg: Springer.

Google Scholar

Cui, X., Bray, S., and Reiss, A. L. (2010). Functional near infrared spectroscopy (NIRS) signal improvement based on negative correlation between oxygenated and deoxygenated hemoglobin dynamics. Neuroimage 49, 3039–3046. doi: 10.1016/j.neuroimage.2009.11.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, R. B., and Harte, D. S. (1987). Tests for hurst effect. Biometrika 74, 95–101. doi: 10.1093/biomet/74.1.95

CrossRef Full Text | Google Scholar

Delignières, D., and Marmelat, V. (2013). Degeneracy and long-range correlations. Chaos 23:043109. doi: 10.1063/1.4825250

PubMed Abstract | CrossRef Full Text | Google Scholar

Devor, A., Dunn, A. K., Andermann, M. L., Ulbert, I., Boas, D. A., and Dale, A. M. (2003). Coupling of total hemoglobin concentration, oxygenation, and neural activity in rat somatosensory cortex. Neuron 39, 353–359. doi: 10.1016/S0896-6273(03)00403-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Drake, C. T., and Iadecola, C. (2007). The role of neuronal signaling in controlling cerebral blood flow. Brain Lang. 102, 141–152. doi: 10.1016/j.bandl.2006.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, S. (2010). Eeg pattern of normal and epileptic rats: monofractal or multifractal? Fractals 18, 425–431. doi: 10.1142/S0218348x10005081

CrossRef Full Text | Google Scholar

Dzung, N. T. (2010). “Multifractality in NIRS data of brain activity,” in IFMBE Proceedings (Ho Chi Minh City), 80–83.

Google Scholar

Eke, A., and Hermán, P. (1999). Fractal analysis of spontaneous fluctuations in human cerebral hemoglobin content and its oxygenation level recorded by NIRS. Adv. Exp. Med. Biol. 471, 49–55. doi: 10.1007/978-1-4615-4717-4_7

PubMed Abstract | CrossRef Full Text | Google Scholar

Eke, A., Hermán, P., and Hajnal, M. (2006). Fractal and noisy CBV dynamics in humans: influence of age and gender. J. Cereb. Blood Flow Metab. 26, 891–898. doi: 10.1038/sj.jcbfm.9600243

PubMed Abstract | CrossRef Full Text | Google Scholar

Eke, A., Hermán, P., Bassingthwaighte, J. B., Raymond, G. M., Percival, D. M., Cannon, M., et al. (2000). Physiological time series: distinguishing fractal noises from motions. Pflügers Archiv. Eur. J. Physiol. 439, 403–415. doi: 10.1007/s004249900135

PubMed Abstract | CrossRef Full Text | Google Scholar

Eke, A., Herman, P., Kocsis, L., and Kozak, L. R. (2002). Fractal characterization of complexity in temporal physiological signals. Physiol. Meas. 23, R1–R38. doi: 10.1088/0967-3334/23/1/201

PubMed Abstract | CrossRef Full Text | Google Scholar

Eke, A., Herman, P., Sanganahalli, B. G., Hyder, F., Mukli, P., and Nagy, Z. (2012). Pitfalls in fractal time series analysis: Fmri BOLD as an exemplary case. Front. Physiol. 3:417. doi: 10.3389/fphys.2012.00417

PubMed Abstract | CrossRef Full Text | Google Scholar

Elwell, C. E., Springett, R., Hillman, E., and Delpy, D. T. (1999). Oscillations in cerebral haemodynamics. Implications for functional activation studies. Adv. Exp. Med. Biol. 471, 57–65.

PubMed Abstract | Google Scholar

Fabiani, M., Gordon, B. A., Maclin, E. L., Pearson, M. A., Brumback-Peltz, C. R., Low, K. A., et al. (2014). Neurovascular coupling in normal aging: a combined optical, ERP and fMRI study. Neuroimage 85 (Pt 1), 592–607. doi: 10.1016/j.neuroimage.2013.04.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, L. K., and Busatto, G. F. (2013). Resting-state functional connectivity in normal brain aging. Neurosci. Biobehav. Rev. 37, 384–400. doi: 10.1016/j.neubiorev.2013.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Firbank, M., Okada, E., and Delpy, D. T. (1998). A theoretical study of the signal contribution of regions of the adult head to near-infrared spectroscopy studies of visual evoked responses. Neuroimage 8, 69–78. doi: 10.1006/nimg.1998.0348

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, M. D., and Raichle, M. E. (2007). Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700–711. doi: 10.1038/nrn2201

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, W. J., Holmes, M. D., Burke, B. C., and Vanhatalo, S. (2003). Spatial spectra of scalp EEG and EMG from awake humans. Clin. Neurophysiol. 114, 1053–1068. doi: 10.1016/s1388-2457(03)00045-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, U., and Parisi, G. (1985). “Turbulence and predictability in geophysical fluid dynamics and climate dynamics,” in Fully Developed Turbulence and Intermittency Appendix: On the Singularity Structure of Fully Developed Structure, eds M. Ghil, R. Benzi, and G. Parisi (Amsterdam: North-Holland), 823.

Ge, E., and Leung, Y. (2012). Detection of crossover time scales in multifractal detrended fluctuation analysis. J. Geogr. Syst. 15, 115–147. doi: 10.1007/s10109-012-0169-9

CrossRef Full Text | Google Scholar

Gibbs, R. B., and Aggarwal, P. (1998). Estrogen and basal forebrain cholinergic neurons: implications for brain aging and Alzheimer's disease-related cognitive decline. Horm. Behav. 34, 98–111. doi: 10.1006/hbeh.1998.1451

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L., Amaral, L. A., Hausdorff, J. M., Ivanov, P. C. h., Peng, C. K., and Stanley, H. E. (2002). Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. U.S.A. 99 (Suppl. 1), 2466–2472. doi: 10.1073/pnas.012579499

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez-Extremera, M., Carpena, P., Ivanov, P. c. h., and Bernaola-Galván, P. A. (2016). Magnitude and sign of long-range correlated time series: Decomposition and surrogate signal generation. Phys. Rev. E 93:042201. doi: 10.1103/PhysRevE.93.042201

PubMed Abstract | CrossRef Full Text | Google Scholar

Grech, D., and Pamuła, G. (2012). Multifractal background noise of monofractal signals. Acta Phys. Pol. 121(2 B), B34–B39. doi: 10.12693/APhysPolA.121.B-34

CrossRef Full Text | Google Scholar

He, B. J. (2011). Scale-free properties of the functional magnetic resonance imaging signal during rest and task. J. Neurosci. 31, 13786–13795. doi: 10.1523/JNEUROSCI.2111-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B. J. (2014). Scale-free brain activity: past, present, and future. Trends Cogn. Sci. 18, 480–487. doi: 10.1016/j.tics.2014.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B. J., Zempel, J. M., Snyder, A. Z., and Raichle, M. E. (2010). The temporal structures and functional significance of scale-free brain activity. Neuron 66, 353–369. doi: 10.1016/j.neuron.2010.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Herman, P., Kocsis, L., and Eke, A. (2009). Fractal characterization of complexity in dynamic signals: application to cerebral hemodynamics. Methods Mol. Biol. 489, 23–40. doi: 10.1007/978-1-59745-543-5_2

PubMed Abstract | CrossRef Full Text | Google Scholar

Herman, P., Sanganahalli, B. G., Hyder, F., and Eke, A. (2011). Fractal analysis of spontaneous fluctuations of the BOLD signal in rat brain. NeuroImage 58, 1060–1069. doi: 10.1016/j.neuroimage.2011.06.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Iadecola, C., and Nedergaard, M. (2007). Glial regulation of the cerebral microvasculature. Nat. Neurosci. 10, 1369–1376. doi: 10.1038/nn2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ihlen, E. A., and Vereijken, B. (2010). Interaction-dominant dynamics in human cognition: beyond 1/f(alpha) fluctuation. J. Exp. Psychol. Gen. 139, 436–463. doi: 10.1037/a0019098

PubMed Abstract | CrossRef Full Text | Google Scholar

Ihlen, E. A., and Vereijken, B. (2013). Multifractal formalisms of human behavior. Hum. Mov. Sci. 32, 633–651. doi: 10.1016/j.humov.2013.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C. h., Ma, Q. D., Bartsch, R. P., Hausdorff, J. M., Amaral, L. A., Schulte-Frohlinde, V., et al. (2009). Levels of complexity in scale-invariant neural signals. Phys. Rev. E 79 (4 Pt 1):041920. doi: 10.1103/PhysRevE.79.041920

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Nunes Amaral, L. A., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Stanley, H. E., et al. (2001). From 1/f noise to multifractal cascades in heartbeat dynamics. Chaos 11, 641–652. doi: 10.1063/1.1395631

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Z. Q., Gao, X. L., Zhou, W.-X., and Stanley, H. E. (2016). Multifractal cross wavelet analysis. Fractals 25:1750054. doi: 10.1142/S0218348X17500542

CrossRef Full Text | Google Scholar

Jöbsis, F. F. (1977). Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters. Science 198, 1264–1267. doi: 10.1126/science.929199

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., and Stanley, H. E. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A 316, 87–114. doi: 10.1016/S0378-4371(02)01383-3

CrossRef Full Text | Google Scholar

Khan, M. M., Dhandapani, K. M., Zhang, Q. G., and Brann, D. W. (2013). Estrogen regulation of spine density and excitatory synapses in rat prefrontal and somatosensory cerebral cortex. Steroids 78, 614–623. doi: 10.1016/j.steroids.2012.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Khoa, T. Q., and Nakagawa, M. (2008). Recognizing brain activities by functional near-infrared spectroscope signal analysis. Nonlinear Biomed. Phys. 2:3. doi: 10.1186/1753-4631-2-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kocsis, L., Herman, P., and Eke, A. (2006a). Mathematical model for the estimation of hemodynamic and oxygenation variables by tissue spectroscopy. J. Theor. Biol. 241, 262–275. doi: 10.1016/j.jtbi.2005.11.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Kocsis, L., Herman, P., and Eke, A. (2006b). The modified Beer-Lambert law revisited. Phys. Med. Biol. 51, N91–N98. doi: 10.1088/0031-9155/51/5/N02

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristoufek, L. (2011). Multifractal height cross-correlation analysis: a new method for analyzing long-range cross-correlations. Epl 95:6. doi: 10.1209/0295-5075/95/68001

CrossRef Full Text | Google Scholar

Kuznetsov, N., Bonnette, S., Gao, J., and Riley, M. A. (2013). Adaptive fractal analysis reveals limits to fractal scaling in center of pressure trajectories. Ann. Biomed. Eng. 41, 1646–1660. doi: 10.1007/s10439-012-0646-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Zhang, M., Xin, Q., Luo, S., Cui, R., Zhou, W., et al. (2013). Age-related changes in spontaneous oscillations assessed by wavelet transform of cerebral oxygenation and arterial blood pressure signals. J. Cereb. Blood Flow Metab. 33, 692–699. doi: 10.1038/jcbfm.2013.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Lima, G. Z. d. S., Corso, G., Correa, M. A., Sommer, R. L., Ivanov, P. C., and Bohn, F. (2017). Universal temporal characteristics and vanishing of multifractality in Barkhausen avalanches. Phys. Rev. E 96:022159. doi: 10.1103/PhysRevE.96.022159

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370–1377. doi: 10.1523/JNEUROSCI.21-04-01370.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Lipsitz, L. A. (2003). “Aging as a process of complexity loss,” in Complex Systems Science in Biomedicine, eds T. S. Deisboeck, D. Kresh, and J. K. Kepler (Boston, MA: Springer), 641–654.

Google Scholar

Lo, C. C., Amaral, L. A. N., Havlin, S., Ivanov, P. C., Penzel, T., Peter, J. H., et al. (2002). Dynamics of sleep-wake transitions during sleep. Europhys. Lett. 57, 625–631. doi: 10.1209/epl/i2002-00508-7

CrossRef Full Text | Google Scholar

Ludescher, J., Bogachev, M. I., Kantelhardt, J. W., Schumann, A. Y., and Bunde, A. (2011). On spurious and corrupted multifractality: the effects of additive noise, short-term memory and periodic trends. Physica A 390, 2480–2490. doi: 10.1016/j.physa.2011.03.008

CrossRef Full Text | Google Scholar

Mandelbrot, B. B. (1974). Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier. J. Fluid Mech. 62, 331–358. doi: 10.1017/S0022112074000711

CrossRef Full Text | Google Scholar

Maxim, V., Sendur, L., Fadili, J., Suckling, J., Gould, R., Howard, R., et al. (2005). Fractional Gaussian noise, functional MRI and Alzheimer's disease. NeuroImage 25, 141–158. doi: 10.1016/j.neuroimage.2004.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

McGeer, P. L., McGeer, E. G., Suzuki, J., Dolman, C. E., and Nagai, T. (1984). Aging, Alzheimer's disease, and the cholinergic system of the basal forebrain. Neurology 34, 741–745. doi: 10.1212/WNL.34.6.741

PubMed Abstract | CrossRef Full Text | Google Scholar

Mevel, K., Landeau, B., Fouquet, M., La Joie, R., Villain, N., Mezengé, F., et al. (2013). Age effect on the default mode network, inner thoughts, and cognitive abilities. Neurobiol. Aging 34, 1292–1301. doi: 10.1016/j.neurobiolaging.2012.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Monto, S. (2012). Nested synchrony-a novel cross-scale interaction among neuronal oscillations. Front. Physiol. 3:384. doi: 10.3389/fphys.2012.00384

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukli, P., Nagy, Z., and Eke, A. (2015). Multifractal formalism by enforcing the universal behavior of scaling functions. Physica A 417, 150–167. doi: 10.1016/j.physa.2014.09.002

CrossRef Full Text | Google Scholar

Muzy, J. F., Bacry, E., and Arneodo, A. (1993). Multifractal formalism for fractal signals: the structure-function approach versus the wavelet-transform modulus-maxima method. Phys. Rev. E Stat. Nonlin. Soft. Matter. Phys. 47, 875–884. doi: 10.1103/PhysRevE.47.875

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagy, Z., Mukli, P., Herman, P., and Eke, A. (2017). Decomposing multifractal crossovers. Front. Physiol. 8:00533. doi: 10.3389/fphys.2017.00533

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, H. J., Huang, X. L., Ning, X. B., Huo, C. Y., Liu, T. B., and Ben, D. (2014). Multifractal analysis of resting state fMRI series in default mode network: age and gender effects. Chinese Sci. Bull. 59, 3107–3113. doi: 10.1007/s11434-014-0355-x

CrossRef Full Text | Google Scholar

Novi, S. L., Rodrigues, R. B., and Mesquita, R. C. (2016). Resting state connectivity patterns with near-infrared spectroscopy data of the whole head. Biomed. Optics Express 7, 2524–2537. doi: 10.1364/BOE.7.002524

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunes Amaral, L. A., Ivanov, P. C., Aoyagi, N., Hidaka, I., Tomono, S., Goldberger, A. L., et al. (2001). Behavioral-independent features of complex heartbeat dynamics. Phys. Rev. Lett. 86(26 Pt 1), 6026–6029. doi: 10.1103/PhysRevLett.86.6026

PubMed Abstract | CrossRef Full Text | Google Scholar

Obrig, H., Neufang, M., Wenzel, R., Kohl, M., Steinbrink, J., Einhäupl, K., et al. (2000). Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage 12, 623–639. doi: 10.1006/nimg.2000.0657

PubMed Abstract | CrossRef Full Text | Google Scholar

Pierro, M. L., Sassaroli, A., Bergethon, P. R., Ehrenberg, B. L., and Fantini, S. (2012). Phase-amplitude investigation of spontaneous low-frequency oscillations of cerebral hemodynamics with near-infrared spectroscopy: a sleep study in human subjects. Neuroimage 63, 1571–1584. doi: 10.1016/j.neuroimage.2012.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Podobnik, B., and Stanley, H. E. (2008). Detrended cross-correlation analysis: a new method for analyzing two nonstationary time series. Phys. Rev. Lett. 100:084102. doi: 10.1103/PhysRevLett.100.084102

PubMed Abstract | CrossRef Full Text | Google Scholar

Quang Dang Khoa, T., and Van Toi, V. (2012). Multifractals Properties on the Near Infrared Spectroscopy of Human Brain Hemodynamic. Math. Probl. Eng. 2012, 1–12. doi: 10.1155/2012/670761

CrossRef Full Text | Google Scholar

Racz, F. S., Mukli, P., Nagy, Z., and Eke, A. (2017). Increased prefrontal cortex connectivity during cognitive challenge assessed by fNIRS imaging. Biomed. Optics Express 8, 3842–3855. doi: 10.1364/BOE.8.003842

PubMed Abstract | CrossRef Full Text | Google Scholar

Racz, F. S., Mukli, P., Nagy, Z., and Eke, A. (2018). Multifractal dynamics of resting-state functional connectivity in the prefrontal cortex. Physiol. Meas. 39:024003. doi: 10.1088/1361-6579/aaa916

PubMed Abstract | CrossRef Full Text | Google Scholar

Raichle, M. E., and Mintun, M. A. (2006). Brain work and brain imaging. Annu. Rev. Neurosci. 29, 449–476. doi: 10.1146/annurev.neuro.29.051605.112819

PubMed Abstract | CrossRef Full Text | Google Scholar

Rattan, S. I. S. (2014). Aging is not a disease: implications for intervention. Aging Disease 5, 196–202. doi: 10.14336/AD.2014.0500196

CrossRef Full Text | Google Scholar

Reinhard, M., Wehrle-Wieland, E., Grabiak, D., Roth, M., Guschlbauer, B., Timmer, J., et al. (2006). Oscillatory cerebral hemodynamics–the macro- vs. microvascular level. J. Neurol. Sci. 250, 103–109. doi: 10.1016/j.jns.2006.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, S. G., Venugopal, V., Fienberg, K., Arneodo, A., and Foufoula-Georgiou, E. (2009). Evidence for inherent nonlinearity in temporal rainfall. Adv. Water Resour. 32, 41–48. doi: 10.1016/j.advwatres.2008.09.007

CrossRef Full Text | Google Scholar

Sala-Llonch, R., Bartrés-Faz, D., and Junqué, C. (2015). Reorganization of brain networks in aging: a review of functional connectivity studies. Front. Psychol. 6:663. doi: 10.3389/fpsyg.2015.00663

PubMed Abstract | CrossRef Full Text | Google Scholar

Sasai, S., Homae, F., Watanabe, H., and Taga, G. (2011). Frequency-specific functional connectivity in the brain during resting state revealed by NIRS. Neuroimage 56, 252–257. doi: 10.1016/j.neuroimage.2010.12.075

PubMed Abstract | CrossRef Full Text | Google Scholar

Sassaroli, A., Pierro, M., Bergethon, P. R., and Fantini, S. (2012). Low-frequency spontaneous oscillations of cerebral hemodynamics investigated with near-infrared spectroscopy: a review. IEEE J. Sel. Topics Quant. Elect. 18, 1478–1492. doi: 10.1109/jstqe.2012.2183581

CrossRef Full Text | Google Scholar

Saupe, D. (1988). “Algorithms for random fractals,” in The Science of Fractal Images, eds H. O. Peitgen and D. Saupe. (New York, NY: Springer-Verlag), 71–136.

Google Scholar

Scholkmann, F., Kleiser, S., Metz, A. J., Zimmermann, R., Mata Pavia, J., Wolf, U., et al. (2014). A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology. NeuroImage 85, 6–27. doi: 10.1016/j.neuroimage.2013.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroeter, M. L., Schmiedel, O., and Von Cramon, D. Y. (2004). Spontaneous low-frequency oscillations decline in the aging brain. J. Cereb. Blood Flow Metab. 24, 1183–1191. doi: 10.1097/01.WCB.0000135231.90164.40

PubMed Abstract | CrossRef Full Text | Google Scholar

Schumann, A. Y., and Kantelhardt, J. W. (2011). Multifractal moving average analysis and test of multifractal model with tuned correlations. Phys. A 390, 2637–2654. doi: 10.1016/j.physa.2011.03.002

CrossRef Full Text | Google Scholar

Shimizu, Y., Barth, M., Windischberger, C., Moser, E., and Thurner, S. (2004). Wavelet-based multifractal analysis of fMRI time series. Neuroimage 22, 1195–1202. doi: 10.1016/j.neuroimage.2004.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, L. E., Silva, C. A., Salgado, H. C., and Fazan, R. (2017). The role of sympathetic and vagal cardiac control on complexity of heart rate dynamics. Am. J. Physiol. Heart Circ. Physiol. 312, H469–H477. doi: 10.1152/ajpheart.00507.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Sleimen-Malkoun, R., Temprado, J. J., and Hong, S. L. (2014). Aging induced loss of complexity and dedifferentiation: consequences for coordination dynamics within and between brain, muscular and behavioral levels. Front. Aging Neurosci. 6:140. doi: 10.3389/fnagi.2014.00140

PubMed Abstract | CrossRef Full Text | Google Scholar

Sowell, E. R., Peterson, B. S., Thompson, P. M., Welcome, S. E., Henkenius, A. L., and Toga, A. W. (2003). Mapping cortical change across the human life span. Nat. Neurosci. 6, 309–315. doi: 10.1038/nn1008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sporns, O. (2011). The human connectome: a complex network. Ann. N. Y. Acad. Sci. 1224, 109–25. doi: 10.1111/j.1749-6632.2010.05888.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stam, C. J. (2005). Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. Clin. Neurophysiol. 116, 2266–2301. doi: 10.1016/j.clinph.2005.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanley, H. E., Amaral, L. A., Goldberger, A. L., Havlin, S., Ivanov, P. C. h., and Peng, C. K. (1999). Statistical physics and physiology: monofractal and multifractal approaches. Phys. A 270, 309–324. doi: 10.1016/S0378-4371(99)00230-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefanovska, A., Bracic, M., and Kvernmo, H. D. (1999). Wavelet analysis of oscillations in the peripheral blood circulation measured by laser Doppler technique. IEEE Trans. Biomed. Eng. 46, 1230–1239.

PubMed Abstract | Google Scholar

Sun, J., Tong, S., and Yang, G. Y. (2012). Reorganization of brain networks in aging and age-related diseases. Aging Dis. 3, 181–193.

PubMed Abstract | Google Scholar

Tachtsidis, I., Leung, T. S., Devoto, L., Delpy, D. T., and Elwell, C. E. (2008). Measurement of frontal lobe functional activation and related systemic effects: a near-infrared spectroscopy investigation. Adv. Exp. Med. Biol. 614, 397–403. doi: 10.1007/978-0-387-74911-2_44

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarantini, S., Tran, C. H. T., Gordon, G. R., Ungvari, Z., and Csiszar, A. (2017). Impaired neurovascular coupling in aging and Alzheimer's disease: contribution of astrocyte dysfunction and endothelial impairment to cognitive decline. Exp. Gerontol. 94, 52–58. doi: 10.1016/j.exger.2016.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Tebaldi, C., De Menech, M., and Stella, A. L. (1999). Multifractal scaling in the Bak-Tang-Wiesenfeld sandpile and edge events. Phys. Rev. Lett. 83, 3952–3955. doi: 10.1103/PhysRevLett.83.3952

CrossRef Full Text | Google Scholar

Thurner, S., Windischberger, C., Moser, E., Walla, P., and Barth, M. (2003). Scaling laws and persistence in human brain activity. Phys. A 326, 511–521. doi: 10.1016/S0378-4371(03)00279-6

CrossRef Full Text | Google Scholar

Tian, F., Chance, B., and Liu, H. (2009). Investigation of the prefrontal cortex in response to duration-variable anagram tasks using functional near-infrared spectroscopy. J. Biomed. Opt. 14:054016. doi: 10.1117/1.3241984

PubMed Abstract | CrossRef Full Text | Google Scholar

Toda, N. (2012). Age-related changes in endothelial function and blood flow regulation. Pharmacol Ther. 133, 159–176. doi: 10.1016/j.pharmthera.2011.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Tulppo, M. P., Kiviniemi, A. M., Hautala, A. J., Kallio, M., Seppänen, T., Mäkikallio, T. H., et al. (2005). Physiological background of the loss of fractal heart rate dynamics. Circulation 112, 314–319. doi: 10.1161/CIRCULATIONAHA.104.523712

PubMed Abstract | CrossRef Full Text | Google Scholar

Valencia, M., Artieda, J., Alegre, M., and Maza, D. (2008). Influence of filters in the detrended fluctuation analysis of digital electroencephalographic data. J. Neurosci. Methods 170, 310–316. doi: 10.1016/j.jneumeth.2008.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Vandeput, S., Verheyden, B., Aubert, A. E., and Van Huffel, S. (2012). Nonlinear heart rate dynamics: circadian profile and influence of age and gender. Med. Eng. Phys. 34, 108–117. doi: 10.1016/j.medengphy.2011.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Vermeij, A., Meel-van den Abeelen, A. S., Kessels, R. P., van Beek, A. H., and Claassen, J. A. (2014). Very-low-frequency oscillations of cerebral hemodynamics and blood pressure are affected by aging and cognitive load. Neuroimage 85 (Pt 1), 608–615. doi: 10.1016/j.neuroimage.2013.04.107

PubMed Abstract | CrossRef Full Text | Google Scholar

Wardlaw, J. M., Smith, C., and Dichgans, M. (2013). Mechanisms of sporadic cerebral small vessel disease: insights from neuroimaging. Lancet Neurol. 12, 483–497. doi: 10.1016/s1474-4422(13)70060-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of “small-world” networks. Nature 393, 440–442. doi: 10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

West, B. J. (1991). Fractal Physiology and Chaos in Medicine. World Scientific Pub. Co. Inc.

Google Scholar

Wink, A. M., Bullmore, E., Barnes, A., Bernard, F., and Suckling, J. (2008). Monofractal and multifractal dynamics of low frequency endogenous brain oscillations in functional MRI. Hum. Brain Mapp. 29, 791–801. doi: 10.1002/hbm.20593

PubMed Abstract | CrossRef Full Text | Google Scholar

Wylie, G. R., Graber, H. L., Voelbel, G. T., Kohl, A. D., DeLuca, J., Pei, Y., et al. (2009). Using co-variations in the Hb signal to detect visual activation: a near infrared spectroscopic imaging study. Neuroimage 47, 473–481. doi: 10.1016/j.neuroimage.2009.04.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada, T., Umeyama, S., and Matsuda, K. (2012). Separation of fNIRS Signals into Functional and Systemic Components Based on Differences in Hemodynamic Modalities. PLoS ONE 7:0050271. doi: 10.1371/journal.pone.0050271

PubMed Abstract | CrossRef Full Text | Google Scholar

Zapperi, S., Bækgaard Lauritsen, K., and Stanley, H. E. (1995). Self-organized branching processes: mean-field theory for avalanches. Phys. Rev. Lett. 75, 4071–4074.

PubMed Abstract | Google Scholar

Zhao, L., Li, W., Fenu, A., Podobnik, B., Wang, Y., and Stanley, H. E. (2017). The q-dependent detrended cross-correlation analysis of stock market. arXiv preprint arXiv:1705.01406. doi: 10.1088/1742-5468/aa9db0

CrossRef Full Text | Google Scholar

Zhou, W. X. (2008). Multifractal detrended cross-correlation analysis for two nonstationary signals. Phys. Rev. E Stat. Nonlin. Soft. Matter Phys. 77(6 Pt 2), 066211. doi: 10.1103/PhysRevE.77.066211

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y. S., Tseng, B. Y., Shibata, S., Levine, B. D., and Zhang, R. (2011). Increases in cerebrovascular impedance in older adults. J. Appl. Physiol. 111, 376–381. doi: 10.1152/japplphysiol.01418.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimeo Morais, G. A., Scholkmann, F., Balardin, J. B., Furucho, R. A., de Paula, R. C. V., Biazoli, C. E. Jr., et al. (2018). Non-neuronal evoked and spontaneous hemodynamic changes in the anterior temporal region of the human head may lead to misinterpretations of functional near-infrared spectroscopy signals. Neurophotonics 5:011002. doi: 10.1117/1.NPh.5.1.011002

PubMed Abstract | CrossRef Full Text | Google Scholar