Multiscale detrended cross-correlation coefficient: estimating coupling in non-stationary neurophysiological signals

The brain consists of a vastly interconnected network of regions, the connectome. By estimating the statistical interdependence of neurophysiological time series, we can measure the functional connectivity (FC) of this connectome. Pearson’s correlation (r_P) is a common metric of coupling in FC studies. Yet r_P does not account properly for the non-stationarity of the signals recorded in neuroimaging. In this study, we introduced a novel estimator of coupled dynamics termed multiscale detrended cross-correlation coefficient (MDC₃). Firstly, we showed that MDC₃ had higher accuracy compared to r_P and lagged covariance using simulated time series with known coupling, as well as simulated functional magnetic resonance imaging (fMRI) signals with known underlying structural connectivity. Next, we computed functional brain networks based on empirical magnetoencephalography (MEG) and fMRI. We found that by using MDC₃ we could construct networks of healthy populations with significantly different properties compared to r_P networks. Based on our results, we believe that MDC₃ is a valid alternative to r_P that should be incorporated in future FC studies.

Introduction

Neuroscientific research has undergone a profound transformation in the last 100 years. Berger’s invention of electroencephalography (EEG) (Berger, 1929) made it possible to record and evaluate neural activity in a non-invasive manner. Initially, studies relied on univariate (i.e., single time series) analysis of the brain dynamics. This started to change toward the end of the 20th century with the first functional connectivity (FC) studies (Friston et al., 1993; Biswal et al., 1995). This new field does not rely only on anatomical connections, it rather studies functional connections that can be created between directly or indirectly coupled neuronal populations. In more mathematical terms, the brain regions are considered nodes on a graph, interconnected by edges (Rubinov and Sporns, 2010). These edges are defined by the statistical relationship of the neuronal time series under investigation.

Several different FC estimators have been introduced with Pearson’s correlation (r_P) being one of the first applied in FC studies (Friston et al., 1993; Biswal et al., 1995). Some drawbacks of this method (e.g., unreliable assessment of non-linear relationships) and the growing interest in exploring other aspects of FC, lead to the introduction of newer methodologies such as phase locking value (PLV) (Lachaux et al., 1999; Bruña et al., 2018), phase lag index (PLI) (Stam et al., 2007), synchronization likelihood (SL) (Stam and Van Dijk, 2002) and mutual information (MI) (Steuer et al., 2002; van den Heuvel and Fornito, 2014). The use of different FC estimators can greatly influence the topology of the networks (Lindquist, 2020; Mukli et al., 2021; Stylianou et al., 2021a). Such differences can be especially problematic when non-healthy populations are being investigated, − e.g., in Alzheimer’s disease patients (Jalili, 2016) – complicating the reproducibility and meta-analysis of studies. It is then important that an informed choice should be made for selecting an FC estimator. Nevertheless, r_P is still widely used (Fornito et al., 2016) due to its simplicity and interpretability. An important advantage of r_P is the capacity to identify positive and negative correlations, which is not always the case with other estimators.

Signals can be divided into two categories: (i) stationary and (ii) non-stationary. A time series X_t – where t indicates the discrete time – is completely stationary when the joint probability distributions of {X_t1, X_t2, X_t3 …, X_tn} and {X_t1 + k, X_t2 + k, X_t3 + k …, X_tn + k} are identical for any set of time points t1, t2, t3…, tn and any integer k. While this definition is easily understood, it is rather unrealistic. Hence, a less strict definition for weak stationarity has been used to classify physiological signals. According to this, the mean and variance of a time series remain constant. In line with that, the covariance of two weakly stationary signals will also be constant throughout the propagation of time. On the other hand, non-stationary signals have varying mean and variance. Additionally, the covariance between two non-stationary signals will be time-dependent (Priestley, 1988). Figure 1 shows an exemplary case of these weakly-stationary and non-stationary signals. From now on, any reference to stationary signals corresponds to weakly-stationary signals. Most biosignals are non-stationary (Semmlow, 2018). As a result, calculating the r_P – a standardized covariance – of two biosignals can be misleading. A solution to this issue was given with the introduction of the detrended cross-correlation coefficient (DCCC) (Zebende, 2011). DCCC makes use of the averaged variance and covariance of smaller sections of the signals (see Section “Multiscale detrended cross-correlation coefficient” below). In this study, we propose an extension of DCCC termed multiscale detrended cross-correlation coefficient (MDC₃). Contrary to DCCC, the output of MDC₃ does not depend on the scale (window length) resulting in easier interpretation of the results. To show this, we compared MDC₃ to r_P (and its directed equivalent) using simulated time series with: (i) known coupling and (ii) known causal interactions [i.e., effective connectivity (EC)]. We also demonstrated the differences between the two estimators in magnetoencephalography (MEG) and functional magnetic resonance imaging (fMRI) recordings.

Figure 1

Figure 1. Example of weakly-stationary and non-stationary signals generated using auto-regressive fractionally integrated moving-average (ARFIMA) processes (see Simulated time series). The mean and variance of weakly-stationary signals remain constant throughout time, while they vary in non-stationary signals.

Methods

Multiscale detrended cross-correlation coefficient

Before introducing MDC₃ we briefly describe DCCC (Zebende, 2011), upon which MDC₃ is based. DCCC was introduced as a more accurate coupling estimator between non-stationary time series. DCCC is calculated for several scales (s) (or window lengths) as follows. For every scale (window length), the two signals X and Y are divided into N non-overlapping windows¹ of length s. In every window the linear trend is removed, leaving the detrended signals ${\widehat{X}}_i$ and ${\widehat{Y}}_i$ , where i is the index of the window. Detrending is performed in order to counteract (at least partially) any spurious coupling emerging due to autocorrelation effects (Horvatic et al., 2011). Then, the covariance between the two signals and the variances of the two signals are estimated for every window. Finally, the ratio of average covariance and the square root of the product of average variances is calculated. Equation 1 provides the mathematical formulation of these steps.

$\begin{array}{ll}DCC{C}_s=\frac{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{cov}_{{\widehat{X}}_i{\widehat{Y}}_i}}{\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{var}_{{\widehat{X}}_i}\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{var}_{{\widehat{Y}}_i}}}& (1)\end{array}$

DCCC is reminiscent of r_P since both estimators range between −1 and 1 with negative values corresponding to anticorrelation and positive values corresponding to correlation (Podobnik et al., 2011). In 2014 Kristoufek showed that DCCC was more accurate than r_P (Kristoufek, 2014) in synthetic non-stationary signals of known coupling. These results warrant the use of DCCC in FC studies, since neuronal time series are non-stationary (Semmlow, 2018). Unfortunately, the use of a multitude of scales (window lengths) makes it hard to interpret. Figure 2 shows a case where different scales (window lengths) result in different coupling estimation, sometimes even with a different sign. Are the two signals correlated or anticorrelated and to what extent? It is not possible to draw a clear conclusion. We believe that MDC₃ could offer a mathematically sound solution to this problem.

Figure 2

Figure 2. Detrended cross-correlation coefficient (DCCC) values for a 4 s-long pair of MEG signals at different scales (window lengths).

The estimation of MDC₃ starts by calculating DCCC for different scales (window lengths). To avoid any arbitrary choice of scales (window lengths), we define frequencies (f) for which we would like to study the coupling of the time series. These frequencies can be converted to scales (window lengths) using the sampling rate (SR) of the signals (s = SR/f). First the DCCC for every frequency is calculated. Then, the two signals are detrended – in this case as a whole – and their cross-spectral density is estimated. We finally calculate the weighted average of DCCC, based on the relative power of each frequency in the cross-spectral density (Equation 2). The distribution of DCCC – similarly to r_P‘s distribution – can be skewed, so DCCC values are normalized using Fisher’s z transform (Alexander, 1990; Corey et al., 1998) before the calculation of the weighted average.

$\begin{array}{ll}MD{C}_3=tanh\left\{{\displaystyle \sum }_{s=\alpha }^{\omega }{w}_s\left[{tanh}^{-1}\left(DCC{C}_s\right)\right]\right\}& (2)\end{array}$

Where $w_s$ is the weight of every scale, ${tanh}^{-1}$ is the inverse hyperbolic tangent², $tanh$ is the hyperbolic tangent³, α is the minimal scale and ω is the maximal scale.

In its current form MDC₃ cannot construct directed graphs, i.e., the FC matrix obtained is symmetric. We also developed the directed MDC₃ (dMDC₃). For the estimation of dMDC₃ instead of calculating the covariance between ${\widehat{X}}_i$ and ${\widehat{Y}}_i$ , we calculate the lagged covariance (LG). ${\widehat{X}}_i$ is shifted from -L-1 to L-1 datapoints, where L is the length of the two signals in datapoints. Such shifts are term lags, e.g., a − 50 lag means that ${\widehat{X}}_i$ was shifted 50 datapoints earlier than ${\widehat{Y}}_i$ . We then estimate the covariance of the two signals for every lag. Negative lags correspond to the cases where ${\widehat{X}}_i$ is leading and positive lags correspond to the cases where ${\widehat{X}}_i$ is lagging the connection. Here the terms leading/lagging indicate the causal effect or effective connectivity (EC) between the two signals. ${\widehat{X}}_i$ leading ${\widehat{Y}}_i$ means that changes in ${\widehat{X}}_i$ will influence ${\widehat{Y}}_i$ . On the contrary, ${\widehat{X}}_i$ lagging ${\widehat{Y}}_i$ means that changes in ${\widehat{Y}}_i$ will influence ${\widehat{X}}_i$ . For both the leading and lagging cases we estimated the maximal covariance, in absolute terms. As a result, for every connection we had two covariance values, one for when the ${\widehat{X}}_i$ is leading and one for when it is lagging. Details about MDC₃ can be found in Figure 3 and the pseudo-code in Table 1. MATLAB, Python, and R versions of MDC₃ are available at: https://github.com/BrainModes/mdc3.

Figure 3

Figure 3. Demonstration of multiscale detrended cross-correlation coefficient (MDC₃) using a 4 s-long pair of MEG signals with a sampling rate of 1,000 Hz. (A) The two signals (green and purple) are divided into smaller non-overlapping windows of length s, in this example s = 500. (B) Each window is detrended. (C) The variances (upper panel) and covariance (lower panel) are calculated for every window. (D) The detrended cross-correlation coefficient (DCCC) is estimated for several scales (window lengths). The black bar is the DCCC when s = 500. (E) The cross-spectral density of the two time series is calculated. The red asterisks correspond to the frequencies used for the estimation of DCCC, while the blue disk corresponds to 2 Hz (i.e., s = 500). MDC₃ is calculated by taking the weighted average of DCCC, where the weight of each frequency is defined by the relative proportion of its power to the total cross-spectral power.

Table 1

Table 1. Multiscale detrended cross-correlation coefficient (MDC₃) pseudo-code.

Simulated time series

ARFIMA processes

In order to validate the efficacy of MDC₃ we simulated pairs of auto-regressive fractionally integrated moving-average (ARFIMA) processes with known cross-correlation, as in Kristoufek (2014). These series are created as follows:

$\begin{array}{ll}A={\displaystyle \sum }_{n=0}^{100}{\alpha }_n\left(d\right){\epsilon }_{A,t-n}& (3)\end{array}$

$\begin{array}{ll}B={\displaystyle \sum }_{n=0}^{100}{\alpha }_n\left(d\right){\epsilon }_{B,t-n}& (4)\end{array}$

${\epsilon }_A$ is sampled from a standard normal distribution. In order to inject cross-correlation ρ between the two time series, we set ${\epsilon }_B=\rho {\epsilon }_A+\epsilon \sqrt{1-{\rho }^2}$ , with $\epsilon$ being sampled from a standard normal distribution (see Appendix for proof). ${\alpha }_n\left(d\right)=\frac{\Gamma \left(n+d\right)}{\Gamma \left(n+1\right)\Gamma \left(d\right)}$ , where Γ is the gamma function. The parameter $d$ defines the non-stationarity of the simulated signal; $d<0.5$ corresponds to stationary time series, $d\ge 0.5$ corresponds to non-stationary time series. Higher values of $d$ indicate a higher level of non-stationarity.

We wanted to study the coupling for both stationary and non-stationary time series. So we employed the same parameters as Kristoufek (2014): (i) $d=\left[0.1,1.4\right]$ with increments of 0.1 and (ii) $\rho =\left[-0.9,0.9\right]$ with increments of 0.1. To demonstrate the benefits of MDC₃ in real-life neuronal time series, our simulations consisted of two types. The first type aimed to emulate EEG/MEG signals with three different lengths: 1000, 5,000 and 10,000 data points. We assumed that their sampling rate was 250 Hz, corresponding to 4, 20 and 40 s of recordings. MDC₃ was calculated for frequencies between 0.5 and 31 Hz with increments of 0.5. In the second type, we wanted to study how lower sampling rates, seen in fMRI, will affect our methodology. The created signals consisted of 100, 200 and 500 data points. In this case we assumed that the sampling rate was 1 Hz, meaning that the simulated time series corresponded to 100, 200 and 500 s. MDC₃ was calculated for frequencies between 0.01 and 0.12 Hz with increments of 0.01. In both types, the maximum frequencies were selected so there were at least 8 data points in every window. We decided to detrend the time series using a second-degree polynomial, since preliminary analysis showed better results compared to linear detrending. We ran 1,000 simulations for each model.

We wanted to see how closely the two estimators (MDC₃ and r_P) are to the real coupling. For every $d$ , ρ and signal length we calculated the root mean squared error (RMSE) of MDC₃ and r_P. Then, simulations of the same $d$ and signal length were grouped together. As a result, we ended up with 14 pairs (one for each value of d) of 19-points (one for each value of ρ) distributions, for every signal length (see Figure 4 for a graphical representation of the distributions). We compared every pair of distributions using a paired t-test or Wilcoxon signed rank test, depending on the normality of the underlying distributions (evaluated using Lilliefors test). Finally, Benjamini-Hochberg (BH) correction (for each signal length, i.e. 14 p values) (Benjamini and Hochberg, 1995) was used to counteract the effect of multiple comparisons. Throughout the manuscript a comparison was considered statistically significant when BH-adjusted p < 0.05.

Figure 4

Figure 4. Root mean squared error (RMSE) of multiscale detrended cross-correlation coefficient (MDC₃) and Pearson’s correlation for different levels of non-stationairity (d) and signal length (panels A–F). We simulated auto-regressive fractionally integrated moving-average (ARFIMA) processes with varying d, signal length and coupling strength (ρ). ρ was used to estimate the RMSE of MDC₃ and Pearson’s correlation. Pairs of distributions whose difference was statistically significant (Benjamini-Hochberg adjusted p < 0.05) are fully colored.

Simulated fMRI

While ARFIMA processes can create signals with known coupling, they do not represent realistic neuronal time series. We simulated the fMRI of 100 “subjects” using The Virtual Brain (TVB) (Sanz Leon et al., 2013; Schirner et al., 2022). Based on the structural connectivity (SC) matrix of each subject (see next paragraph), we simulated the fMRI signal of 68 brain regions – according to the Desikan-Killiany atlas (Desikan et al., 2006) – using the Reduced Wong Wang (Deco et al., 2013) neural mass model:

$\begin{array}{ll}{x}_k=w{J}_N{S}_k+I_o+J_{N}G{\displaystyle \sum }_j{C}_{kj}{S}_j& (5)\end{array}$

$\begin{array}{ll}H\left(x_k\right)=\frac{a{x}_k-b}{1-exp\left[-d\left(a{x}_k-b\right)\right]}& (6)\end{array}$

$\begin{array}{ll}{\dot{S}}_k=-\frac{S_k}{{\tau }_s}+\left(1-S_k\right)H\left(x_k\right)\gamma & (7)\end{array}$

$H\left(x_k\right)$ and $S_k$ correspond to the firing rate and synaptic gating variable of the population at the k^th cerebral region, respectively. $G$ is a global scaling factor and $C_{kj}$ is the structural connection strength between the k^th and j^th regions. The description and default values of the rest of parameters can be found in Table 12 of Sanz-Leon et al. (2015).

The simulated SC matrices were based the real SC matrix retrieved from https://zenodo.org/record/4263723#.Y7-8Q-zMLMI (found in “QL_20120814_Connectivity.zip”). The real SC matrix was divided into 4 quadrants. The values within each quadrant were randomly shuffled. Additionally, 30% of the connections of each quadrant were changed. Their new value was randomly selected from a normal distribution of mean and standard deviation based on the SC values of each quadrant. This shuffling and random allocation of values was also done in the accompanying tract lengths matrix created after loading “QL_20120814_Connectivity.zip” on TVB. These steps ensured that the simulated brains were different enough from the template, but they were still biologically plausible. We then proceeded with simulating 21 min of fMRI time series using the Reduced Wong Wang model. The selection of appropriate parameters in brain simulations is crucial. A common practice is to perform a grid search with different combinations of parameters and compare it to properties of empirical brain activity. We varied G, w and J, while using the default values of the rest of the parameters. G was in the [0.1, 29.9] range with increments of 0.1. w was in the [0,1] range with increments of 0.1. Finally, J was in the [0.2609, 0.4609] range with increments of 0.05. We estimated the FC matrix of each simulated fMRI dataset using r_P. We also estimated the FC of the empirical fMRI signal using r_P.⁴ We then compared the similarities of empirical and simulated FC using Spearman’s correlation. The most realistic simulation (Spearman’s correlation 0.34) was produced when G = 0.2, w = 0.1 and J = 0.42 while the rest of the parameters were kept in their default values.

Having available the SC matrices for every “subject” allowed us to use dynamic causal modeling (DCM) (Friston et al., 2003) to calculate the EC. Investigation of whole-brain networks with traditional DCM is a time-consuming process, which can be accelerated with regression dynamic causal modeling (rDCM) (Frässle et al., 2017, 2018, 2021b) [available at the Translational Algorithms for Psychiatry-Advancing Science (TAPAS) toolbox (Frässle et al., 2021a)]. rDCM offers a simplified version of DCM without severe loss in accuracy [for further details please see Frässle et al.]. While FC is simple to understand and estimate, it is merely a statistical relationship between signals. On the other hand, rDCM’s constraints allow for a depiction of brain connectivity based on a more realistic network model of the brain. Hence, the EC captured by rDCM was chosen as the ground truth of our comparison. In rDCM a realistic SC connectivity matrix is used as a template. Applying a forward model to the underlying SC can simulate fMRI signals. A parameter of this forward model is an EC matrix, which can be fine-tuned to produce realistic fMRI time series. Since rDCM can capture the direction of the connection, we employed dMDC₃ which we compared with the LG (see Multiscale detrended cross-correlation coefficient). In order to study the effect of signal length we analyzed the first 5, 10, 15 and 20 min of the simulated fMRI. This resulted in 12 matrices (4 signal lengths x 3 metrics) (Table 2) for every simulated brain. Since the rDCM and LG are not constrained between −1 and 1 as dMDC₃, we calculated the Z-scores of every rDCM, dMDC₃ and LG matrix, which we then used for the comparisons. Using rDCM as our ground truth, we calculated the RMSE of dMDC₃ and LG for each simulation. This resulted in 8 (2 EC estimators x 4 signal lengths) 100-point (100 simulated brains) distributions. We compared every pair of distributions using a paired t-test or Wilcoxon signed rank test, depending on the normality of the underlying distributions (evaluated using Lilliefors test). The 4 p values were adjusted using BH correction. dMDC₃ was calculated for the frequencies between 0.011 to 0.17 Hz with increments of 0.01. 0.17 Hz was selected as the highest cutoff so each window during the estimation of dMDC₃ had 8 datapoints. Second-degree polynomials were fitted for the detrending in dMDC3.

Table 2

Table 2. Connectivity matrices used in the analysis of simulated fMRI signals.

Empirical time series

MEG dataset

The MEG dataset consisted of eyes closed resting-state recordings of 20 elderly healthy participants (12 females, aged 71.5 ± 4.03 years), acquired using a 306-channel (102 magnetometers and 204 planar gradiometers) Vectorview MEG system (Elekta AB, Stockholm, Sweden) placed inside a magnetically shielded room (VacuumSchmelze GmbH, Hanau, Germany) located at the Laboratory of Cognitive and Computational Neuroscience (Madrid, Spain). MEG data were acquired with a sampling rate of 1,000 Hz and an online [0.1–330] Hz anti-alias band-pass filter. All participants provided informed consent. To allow subject-specific source reconstruction, individual T1-weighted MRI scans were also available for each participant. MRI images were recorded at the Hospital Universitario Clínico San Carlos (Madrid, Spain) using a 1.5 T General Electric MRI scanner with a high-resolution antenna and a homogenization PURE filter (fast spoiled gradient echo sequence, with parameters: repetition time/echo time/inversion time = 11.2/4.2/450 ms; flip angle = 12°; slice thickness = 1 mm; 256 × 256 matrix; and field of view = 256 mm).

The MEG recordings were preprocessed offline using a tempo-spatial filtering algorithm (Taulu and Hari, 2009) (Maxfilter Software v2.2, correlation limit of 0.9 and correlation window of 10 s) to eliminate magnetic noises and compensate for head movements during the recording. The continuous MEG data were imported into MATLAB (R2017b, Mathworks, Inc.) using the Fieldtrip Toolbox (Oostenveld et al., 2011), where an independent component-based algorithm was used to remove the effects of ocular and cardiac signals from the data, together with external noises. Source reconstruction was performed using minimum norm estimates (Hämäläinen and Ilmoniemi, 1994) with the software Brainstorm (Tadel et al., 2011). In order to model the orientation of macrocolumns of pyramidal neurons the dipole orientations were considered to be normal to the cortical surface of the participant [(see Tadel et al., 2019)]. Neural time series were finally collapsed to the regions of interest (ROI) of the Desikan-Killiany atlas (Desikan et al., 2006) by using the mean operator across all vertex-level constrained time series within that ROI. The data were band-pass filtered between 0.5 and 45 Hz using FIR filtering.

For every participant we analyzed multiple (ranging from 45 to 61) 4 s segments. We estimated the FC of each segment using MDC₃ and r_P. Then, we calculated the node strength of the brain regions by summing up the strength of every incoming and outgoing connection for every cortical area. Finally, we averaged the node strengths for all segments, so every participant had one set of node strength values. Again, we employed a series of paired t-tests or Wilcoxon signed rank tests – depending on the normality of the distributions (Lilliefors test) – to compare the node strengths of the MDC₃ and r_P created networks. The p-values of each comparison group were adjusted using BH correction. MDC₃ was calculated for the frequencies between 0.5 and 45 Hz. Second-degree polynomials were fitted for the detrending in MDC₃.

fMRI dataset

Finally, we analyzed 767 healthy, young adults (426 females) from the Human Connectome Project (HCP) (Smith et al., 2013). The fMRI time series were already preprocessed according to the HCP standards (Glasser et al., 2013). Details about the participants can be found in the attached CSV file in the Supplementary material (fMRI Subjects Information).

For the FC estimation we used only the first eyes open resting-state period of 14.4 min. The dataset had a left-to-right and right-to-left echo-planar imaging (EPI) encoding. We calculated the FC using MDC₃ and r_P for both EPI. We then averaged the FC matrices of the two EPI using Fisher’s z transform, as suggested by Smith et al. (2013). This resulted in one MDC₃ and one r_P FC matrix per subject. We compared the strength of each connection through a series of Wilcoxon signed rank tests that were later corrected using BH. MDC₃ was calculated for the frequencies between 0.011 and 0.17 Hz with increments of 0.01. 0.17 Hz was selected as the highest cutoff, so each window had 8 datapoints. Second degree polynomials were fitted for the detrending in MDC3.

Results

Simulated time series

As shown in Figure 4 MDC₃ is a more accurate estimator of coupling in the simulated ARFIMA signals in almost every case. Only some small difference can be observed for stationary signals (d < 0.5); but as we transition to non-stationary time series (d ≥ 0.5), the RMSE of r_P is significantly higher. All BH-adjusted p-values can be found in the Supplementary material (Statistics).

The same results can be seen in realistic fMRI simulations. As Figure 5 shows, the RMSE was significantly smaller when dMDC₃ was used as an FC estimator in all signal lengths. We also see that as the signal length increases, the RMSE of LG increases while the RMSE of MDC₃ decreases. All BH-adjusted p-values can be found in the Supplementary material (Statistics).

Figure 5

Figure 5. Root mean squared error (RMSE) of directed multiscale detrended cross-correlation coefficient (dMDC₃) and lagged covariance (LG), for four different signal lengths (5 min, 10 min, 15 min and 20 min). We simulated realistic fMRI signals using The Virtual Brain. The effective connectivity of the simulated brains – calculated using regression dynamic causal modeling (rDCM)– was used to estimate the RMSE of MDC₃ and LG.

Neurophysiological time series

Figure 6 shows the difference of the node strengths between the MDC₃ and r_P networks as estimated using MEG tracings. Significant differences can be seen in 7 channels (10%), where the r_P network had mainly higher node strengths seen by the blue color. All BH-adjusted p-values can be found in the Supplementary material (Statistics).

Figure 6

Figure 6. Difference between the node strengths calculated during eyes closed resting-state magnetoencephalography: lateral view (up); medial view (down). The colors represent the difference (MDC₃-r_P) in the node strengths while the numbers indicate the brain regions whose node strength was significantly different between the two estimators (BH-adjusted p < 0.05). The numbers correspond to the regions of interest as defined in the Desikan-Killiany atlas (Desikan et al., 2006), list provided in Supplementary Table S2 (Additional Analysis).

For the last real-life dataset, we analyzed fMRI recordings from HCP. As Figure 7 shows, the two networks had different connectivity strengths. In some instances, r_P found higher coupling than MDC₃ and in some other cases lower. These observations were validated statistically, since 97% (69,599 out of 71,631) of the comparisons were significantly different. All BH-adjusted p-values can be found in the Supplementary material (Statistics).

Figure 7

Figure 7. Averaged functional connectivity matrices using multiscale detrended-cross correlation coefficient (MDC₃), Pearson’s correlation (r_P), and the difference between them (MDC₃-r_P) using eyes open resting-state functional magnetic resonance imaging.

Discussion

In this study we introduced the statistical metric MDC₃ – a weighted average of DCCC – for estimating coupling in a system. Our simulations with signals of known coupling showed that MDC₃ is a more accurate estimator of the model’s coupling parameters than r_P or LG. The exemplary FC analysis of MEG and fMRI data also showed that the use of MDC₃ could lead to significant differences in the connectivity matrices compared to r_P.

We simulated 1,000 pairs of time series of different coupling strengths, signal lengths and degrees of non-stationarity. For each pair we calculated MDC₃ and r_P. As explained in the Introduction, and shown in Figure 1, the variance and covariance of stationary signals remain constant, meaning that MDC₃ and r_P will be similar. This is not the case for non-stationary series whose variance and covariance heavily depend on time. Our simulations confirm that, since the RMSE of MDC₃ was significantly smaller in every case, except for fairly stationary signals (Figure 4). The discrepancy between the two estimators increased greatly with higher levels of non-stationarity. Similar findings have been reported for DCCC in Kristoufek (2014). Several studies (Kristoufek, 2014; Zhang et al., 2020; Guedes et al., 2021; Racz et al., 2024) using DCCC integrate (cumulatively sum) the signals before its estimation. When ARFIMA processes were used as ground truth, the integrating version of MDC₃ performed worse than r_P (see Supplementary Figure S1 in Additional Analysis). A possible explanation is that by integrating the signal, more non-stationary characteristics are introduced, making MDC₃ unreliable. We suggest that future studies take this into account before integrating their signals. We also simulated a series of fMRI signals using TVB. We could not simulate realistic neuronal time series with known coupling, so we decided to use the rDCM matrices of the simulations as ground truth. The results showed that dMDC₃ is closer to rDCM compared to LG (Figure 5). We observed that as the length of the signals increased the accuracy of dMDC₃ slightly increased, while the accuracy of to LG slightly decreased. At this point someone might question the benefits of MDC₃ over rDCM, since we considered rDCM the ground truth. In our opinion, the biggest limitation of rDCM is the need of a SC matrix, which is often not available in several neuroimaging studies. On the contrary, MDC₃ (as well as the rest of FC estimators) can capture functional interdependence without prior knowledge of SC. Smith et al. (2011) decided to validate FC estimators using the underlying SC as ground truth. While we considered this approach, we decided to use EC instead. The choice was based on the two following reasons. Firstly, SC cannot entirely predict FC (Honey et al., 2009). Secondly, the lack of negative values in SC would not allow for accurate study of negatively correlated brain regions. For the sake of completeness, we also compared MDC₃ and r_P of the simulated fMRI signals using SC as ground truth. This time, r_P was found to be a better estimator, albeit with a narrow margin (see Supplementary Figure S2 in Additional Analysis). An interesting byproduct of this analysis was that r_P and MDC₃ were similar to SC, while EC and dMDC₃ were similar to the tract length matrices used for the construction of the simulations. While this finding is interesting, it is beyond the scope of this study and should be revised in future studies. The matrices of each simulation can be found in the Supplementary material (TVB Matrices). Finally, we repeated our MDC₃ and r_P comparisons this time using the simulations from Smith et al. (2011). In most cases MDC₃ was more accurate, especially when EC was used as ground truth. The complete results of the additional analysis can be found in the Supplementary material (Additional Analysis).

Of course, statistical significance in simulations without real-life benefits would not warrant the use of MDC₃. To demonstrate its advantages, we used MEG and fMRI datasets. As shown in Figure 6, using MDC₃ and r_P as FC estimators resulted in significantly different brain networks. In some cases, the node strengths of the r_P networks were higher, while in others they were lower. After analyzing the FC matrices of the fMRI dataset, we saw that almost all connections were significantly different between the two matrices (Figure 7). Once again, some connections were stronger and some weaker when r_P was used. A homogenous overestimation or underestimation would not have been a serious drawback since FC studies usually rely on relative comparisons and not on the exact values themselves. But it seems that in some regions r_P would give lower values and in others higher, presenting a rather false picture of the brain network. At a first glance, someone might be dismissive of this, since it is well known that different estimators can lead to different FC matrices (Jalili, 2016; Mukli et al., 2021; Stylianou et al., 2021a). This would have been the case if we had not seen the higher reliability of MDC₃ in simulations (Figure 4; Figure 5). We then suggest that MDC₃ should be preferred over r_P. Even if MDC₃ is computationally more expensive, today’s computational capabilities make the time difference negligible.

Finally, it should be noted that MDC₃ is still a linear FC estimator. Non-linear estimators like PLV, MI, PLI, and SL still capture dynamics that MDC₃ cannot. In spite of that, we believe that MDC₃ is a valuable addition to the FC field due to its ability to capture the sign of correlation (i.e., correlation vs. anticorrelation); something that the aforementioned non-linear estimators cannot do. A common practice in FC studies is the exclusion of anticorrelations (Rubinov and Sporns, 2010). Since the human brain operates with several negative feedback loops, we believe that it is necessary to study anticorrelation in order to obtain more accurate brain architectures, as suggested by previous studies (Chen et al., 2011; Zhan et al., 2017). We decided to explore this further in the Supplementary material (Additional Analysis) using the MEG dataset. Briefly, we compared the FC matrices as estimated with MDC₃, r_P and PLV using two different source reconstruction pipelines, i.e., with constrained and unconstrained dipoles. The first method makes it possible to obtain a more realistic phase (and sign) of the reconstructed time series. This benefit can be overshadowed by the inability of most FC estimators to capture the sign of coupling, including PLV. As a result, such metrics could mistakenly identify correlation for anticorrelation and vice versa. As expected, both MDC₃ and r_P detected more differences between the reconstructions with constrained and unconstrained dipoles than PLV (Supplementary Figure S3). It then seems that MDC₃ and r_P offer an advantage when the sign of the correlation is important. In the main body of the manuscript, we compared the MDC₃ and r_P using constrained dipoles as source reconstruction for the MEG dataset. In the Supplementary material (Additional Analysis) we repeated this analysis using unconstrained dipoles. The results (Supplementary Figure S4 right panel) showed that even with unconstrained dipoles differences were found between the two metrics. MDC₃ estimated exclusively higher values for all node strengths that were found statistically significant in the unconstrained version. On the contrary, when we constrained the dipoles, we saw both higher and lower values estimated by MDC₃ (Figure 6; Supplementary Figure S4 left panel). It is clear then that irrespective of the source reconstruction method, substantial differences were observed between the two FC estimators. Of course, there are several other preprocessing steps where different pipelines can be implemented. In some of these pipelines MDC₃ and r_P could capture the same dynamics. Nonetheless, considering that MDC₃ had lower RMSE in our simulations (Figure 4; Figure 5), we believe that MDC₃ should be preferred over r_P.

DCCC and its extension MDC₃ are closely related to the scale-free analysis of signals. The numerator and denominator of Equation 1 are integral parts of the detrended fluctuation analysis (Peng et al., 1994) and detrended cross-correlation analysis (Podobnik and Stanley, 2008) analysis, respectively. DCCC has been incorporated in surrogate testing of fractal (scale-free) coupling already (Podobnik et al., 2011; Blythe et al., 2016; Stylianou et al., 2021a, 2021b, 2022). We wanted to explore how the scale-free character of the signals affects the performance of the MDC₃ and r_P [see Supplementary material (Additional Analysis)]. To achieve this, we estimated the Hurst exponent (i.e., degree of autocorrelation or H) of the simulated ARFIMA processes and correlated it with the RMSE of the two FC estimators. We found that H positively correlated with RMSE (Supplementary Figure S5; Supplementary Table S1). As explained in the Supplementary material (Additional Analysis) high autocorrelation corresponds to high non-stationarity. Hence, these results agree with Figure 4 which showed that RMSE increases with higher non-stationarity (and by extension higher H). Despite the almost perfect correlation between RMSE and H for both estimators, RMSE increased in a faster rate when r_P was used. This validates our original findings and shows that MDC₃ is a better option when the scale-free nature of the signals is under consideration. DCCC has also been employed in multifractal FC (Kwapień et al., 2015); where different exponents capture different sizes of fluctuations. Theoretically, a multifractal MDC₃ could be created as well. This is beyond the scope of the current study because we focused on improving the interpretability of DCCC. The calculation of MDC₃ using different scaling exponents would add another layer of complexity to the interpretation of the outputs. Recently, a real-time algorithm for the estimation of DCCC was presented (Kaposzta et al., 2022, 2023), which can be extended for MDC₃ as well. This means that MDC₃ can be used in brain-computer interfaces or clinical monitoring of patients, where constant tracking of network dynamics is needed.

Conclusion

We presented a new estimator of coupling between time series termed multiscale detrended cross-correlation coefficient. Using simulated data, we showed a higher accuracy over r_P and LG. The differences between the estimators were made apparent in MEG and fMRI datasets of healthy populations. Here we explored the benefits of MDC₃ only in neuronal time series. We believe that our new method has the potential to be used in several other disciplines where linear coupling of non-stationary signals is investigated. Of course, appropriate validation pipelines specific to each field are recommended before any prior use.

Data availability statement

The data that support the findings of this study are available from the corresponding authors upon reasonable request. Requests to access these datasets should be directed to Petra Ritter, cGV0cmEucml0dGVyQGJpaC1jaGFyaXRlLmRl.

Ethics statement

The studies involving humans were approved by Ethikkommission Ethikausschuss am Campus Benjamin Franklin, Charite, Berlin. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

OS: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. GS: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Validation, Visualization, Writing – review & editing. MH: Writing – review & editing. IS-M: Data curation, Writing – review & editing. DL-S: Data curation, Writing – review & editing. MS: Data curation, Investigation, Methodology, Writing – review & editing. PR: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. P.R. acknowledges support by Digital Europe TEF-Health 101100700, EU H2020 Virtual Brain Cloud 826421, Human Brain Project SGA2 785907; Human Brain Project SGA3 945539, ERC Consolidator 683049; German Research Foundation SFB 1436 (project ID 425899996); SFB 1315 (project ID 327654276); SFB 936 (project ID 178316478); SFB-TRR 295 (project ID 424778381); SPP Computational Connectomics RI 2073/6–1, RI 2073/10–2, RI 2073/9–1; PHRASE Horizon EIC grant 101058240; Berlin Institute of Health & Foundation Charité, Johanna Quandt Excellence Initiative; ERAPerMed Pattern-Cog, the Virtual Research Environment at the Charité Berlin – a node of EBRAINS Health Data Cloud.

Acknowledgments

A preprint of the study can be found at: https://www.biorxiv.org/content/10.1101/2024.04.16.589689v1.abstract.

Conflict of interest

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

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

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

Supplementary material

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

Footnotes

1. ^Preliminary analysis with 50% overlapping windows did not show significant benefits compared to non-overlapping windows. For the sake of computational speed, non-overlapping windows were chosen.

2. ^ $\mathit{tan}{h}^{-1}\left(x\right)=\frac{1}{2}\mathit{log}\left(\frac{1+x}{1-x}\right)$

3. ^ $\mathit{tanh}\left(x\right)=\frac{e^{2x}-1}{e^{2x}+1}$

4. ^Also retrieved from https://zenodo.org/record/4263723#.Y7-8Q-zMLMI

References

Alexander, R. A. (1990). A note on averaging correlations. Bull. Psychon. Soc. 28, 335–336. doi: 10.3758/BF03334037

Crossref Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B Methodol. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

Crossref Full Text | Google Scholar

Berger, H. (1929). Uber das Elektrenkephalogramm des Menschen (On the human elec- troencephalogram). Arch. F Psychiatr. U Nervenkrankh. 87, 527–570. doi: 10.1007/BF01797193

Crossref Full Text | Google Scholar

Biswal, B., Zerrin Yetkin, F., 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

Blythe, D. A. J., Nikulin, V. V., and Müller, K.-R. (2016). Robust statistical detection of Power-law cross-correlation. Sci. Rep. 6:27089. doi: 10.1038/srep27089

PubMed Abstract | Crossref Full Text | Google Scholar

Bruña, R., Maestú, F., and Pereda, E. (2018). Phase locking value revisited: teaching new tricks to an old dog. J. Neural Eng. 15:056011. doi: 10.1088/1741-2552/aacfe4

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, G. G., Chen, G. G., Xie, C., and Li, S.-J. J. (2011). Negative functional connectivity and its dependence on the shortest path length of positive network in the resting-state human brain. Brain Connect. 1, 195–206. doi: 10.1089/brain.2011.0025

PubMed Abstract | Crossref Full Text | Google Scholar

Corey, D. M., Dunlap, W. P., and Burke, M. J. (1998). Averaging correlations: expected values and bias in combined Pearson rs and fisher’s z transformations. J. Gen. Psychol. 125, 245–261. doi: 10.1080/00221309809595548

Crossref Full Text | Google Scholar

Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann, P., and Corbetta, M. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J. Neurosci. 33, 11239–11252. doi: 10.1523/JNEUROSCI.1091-13.2013

PubMed Abstract | Crossref Full Text | Google Scholar

Desikan, R. S., Ségonne, F., Fischl, B., Quinn, B. T., Dickerson, B. C., Blacker, D., et al. (2006). An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage 31, 968–980. doi: 10.1016/j.neuroimage.2006.01.021

PubMed Abstract | Crossref Full Text | Google Scholar

Fornito, A., Zalesky, A., and Bullmore, E. T. (2016). “Chapter 1 - an introduction to brain networks” in Fundamentals of brain network analysis (San Diego: Academic Press).

Google Scholar

Frässle, S., Aponte, E. A., Bollmann, S., Brodersen, K. H., Do, C. T., Harrison, O. K., et al. (2021a). TAPAS: an open-source software package for translational Neuromodeling and computational psychiatry. Front. Psych. 12:680811. doi: 10.3389/fpsyt.2021.680811

PubMed Abstract | Crossref Full Text | Google Scholar

Frässle, S., Harrison, S. J., Heinzle, J., Clementz, B. A., Tamminga, C. A., Sweeney, J. A., et al. (2021b). Regression dynamic causal modeling for resting-state fMRI. Hum. Brain Mapp. 42, 2159–2180. doi: 10.1002/hbm.25357

PubMed Abstract | Crossref Full Text | Google Scholar

Frässle, S., Lomakina, E. I., Kasper, L., Manjaly, Z. M., Leff, A., Pruessmann, K. P., et al. (2018). A generative model of whole-brain effective connectivity. NeuroImage 179, 505–529. doi: 10.1016/j.neuroimage.2018.05.058

Crossref Full Text | Google Scholar

Frässle, S., Lomakina, E. I., Razi, A., Friston, K. J., Buhmann, J. M., and Stephan, K. E. (2017). Regression DCM for fMRI. NeuroImage 155, 406–421. doi: 10.1016/j.neuroimage.2017.02.090

Crossref Full Text | Google Scholar

Friston, K. J., Frith, C. D., Liddle, P. F., and Frackowiak, R. S. J. (1993). Functional connectivity: the principal-component analysis of large (PET) data sets. J. Cereb. Blood Flow Metab. 13, 5–14. doi: 10.1038/jcbfm.1993.4

PubMed Abstract | Crossref Full Text | Google Scholar

Friston, K. J., Harrison, L., and Penny, W. (2003). Dynamic causal modelling. NeuroImage 19, 1273–1302. doi: 10.1016/S1053-8119(03)00202-7

Crossref Full Text | Google Scholar

Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., et al. (2013). The minimal preprocessing pipelines for the human connectome project. NeuroImage 80, 105–124. doi: 10.1016/j.neuroimage.2013.04.127

PubMed Abstract | Crossref Full Text | Google Scholar

Guedes, E. F., da Silva Filho, A. M., and Zebende, G. F. (2021). Detrended multiple cross-correlation coefficient with sliding windows approach. Phys. Stat. Mech. Its Appl. 574:125990. doi: 10.1016/j.physa.2021.125990

Crossref Full Text | Google Scholar

Hämäläinen, M. S., and Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. doi: 10.1007/BF02512476

PubMed Abstract | Crossref Full Text | Google Scholar

Honey, C. J., Honey, C. J., Sporns, O., Sporns, O., Cammoun, L., Cammoun, L., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. Proc. Natl. Acad. Sci. USA 106, 2035–2040. doi: 10.1073/pnas.0811168106

PubMed Abstract | Crossref Full Text | Google Scholar

Horvatic, D., Stanley, H. E., and Podobnik, B. (2011). Detrended cross-correlation analysis for non-stationary time series with periodic trends. EPL Europhys. Lett. 94:18007. doi: 10.1209/0295-5075/94/18007

Crossref Full Text | Google Scholar

Jalili, M. (2016). Functional brain networks: does the choice of dependency estimator and Binarization method matter? Sci. Rep. 6:29780. doi: 10.1038/srep29780

PubMed Abstract | Crossref Full Text | Google Scholar

Kaposzta, Z., Czoch, A., Mukli, P., Stylianou, O., Liu, D. H., Eke, A., et al. (2023). Fingerprints of decreased cognitive performance on fractal connectivity dynamics in healthy aging. GeroScience 46, 713–736. doi: 10.1007/s11357-023-01022-x

PubMed Abstract | Crossref Full Text | Google Scholar

Kaposzta, Z., Czoch, A., Stylianou, O., Kim, K., Mukli, P., Eke, A., et al. (2022). Real-time algorithm for Detrended cross-correlation analysis of long-range coupled processes. Front. Physiol. 13:817268. doi: 10.3389/fphys.2022.817268

PubMed Abstract | Crossref Full Text | Google Scholar

Kristoufek, L. (2014). Measuring correlations between non-stationary series with DCCA coefficient. Phys. Stat. Mech. Its Appl. 402, 291–298. doi: 10.1016/j.physa.2014.01.058

Crossref Full Text | Google Scholar

Kwapień, J., Oświęcimka, P., and Drożdż, S. (2015). Detrended fluctuation analysis made flexible to detect range of cross-correlated fluctuations. Phys. Rev. E 92:052815. doi: 10.1103/PhysRevE.92.052815

Crossref Full Text | Google Scholar

Lachaux, J.-P., Rodriguez, E., Martinerie, J., and Varela, F. J. (1999). Measuring phase synchrony in brain signals. Hum. Brain Mapp. 8, 194–208. doi: 10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C

Crossref Full Text | Google Scholar

Lindquist, M. (2020). Neuroimaging results altered by varying analysis pipelines. Nature 582, 36–37. doi: 10.1038/d41586-020-01282-z

PubMed Abstract | Crossref Full Text | Google Scholar

Mukli, P., Nagy, Z., Racz, F. S., Portoro, I., Hartmann, A., Stylianou, O., et al. (2021). Two-tiered response of cardiorespiratory-cerebrovascular network to orthostatic challenge. Front. Physiol. 12:622569. doi: 10.3389/fphys.2021.622569

PubMed Abstract | Crossref Full Text | Google Scholar

Oostenveld, R., Fries, P., Maris, E., and Schoffelen, J.-M. (2011). FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011, 1–9. doi: 10.1155/2011/156869

PubMed Abstract | Crossref Full Text | Google Scholar

Peng, C. K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., and Goldberger, A. L. (1994). Mosaic organization of DNA nucleotides. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 49, 1685–1689. doi: 10.1103/physreve.49.1685

Crossref Full Text | Google Scholar

Podobnik, B., Jiang, Z.-Q., Zhou, W.-X., and Stanley, H. E. (2011). Statistical tests for power-law cross-correlated processes. Phys. Rev. E 84:066118. doi: 10.1103/PhysRevE.84.066118

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

Priestley, M. B. (1988). Non-linear and non-stationary time series analysis. Available at: https://ui.adsabs.harvard.edu/abs/1988nlns.book.....P (Accessed June 21, 2023).

Google Scholar

Racz, F. S., Kumar, S., Kaposzta, Z., Alawieh, H., Liu, D. H., Liu, R., et al. (2024). Combining detrended cross-correlation analysis with Riemannian geometry-based classification for improved brain-computer interface performance. Front. Neurosci. 18:1271831. doi: 10.3389/fnins.2024.1271831

PubMed Abstract | Crossref Full Text | Google Scholar

Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. NeuroImage 52, 1059–1069. doi: 10.1016/j.neuroimage.2009.10.003

PubMed Abstract | Crossref Full Text | Google Scholar

Sanz Leon, P., Knock, S., Woodman, M., Domide, L., Mersmann, J., McIntosh, A., et al. (2013). The virtual brain: a simulator of primate brain network dynamics. Front. Neuroinform. 7:10. doi: 10.3389/fninf.2013.00010

Crossref Full Text | Google Scholar

Sanz-Leon, P., Knock, S. A., Spiegler, A., and Jirsa, V. K. (2015). Mathematical framework for large-scale brain network modeling in the virtual brain. NeuroImage 111, 385–430. doi: 10.1016/j.neuroimage.2015.01.002

PubMed Abstract | Crossref Full Text | Google Scholar

Schirner, M., Domide, L., Perdikis, D., Triebkorn, P., Stefanovski, L., Pai, R., et al. (2022). Brain simulation as a cloud service: the virtual brain on EBRAINS. NeuroImage 251:118973. doi: 10.1016/j.neuroimage.2022.118973

PubMed Abstract | Crossref Full Text | Google Scholar

Semmlow, J. (2018). “Chapter 10 - stochastic, nonstationary, and nonlinear systems and signals” in Circuits, signals and Systems for Bioengineers (third edition). ed. J. Semmlow (Academic Press), 449–489. Available at: https://www.sciencedirect.com/science/article/abs/pii/B9780128093955000102

Google Scholar

Smith, S. M., Beckmann, C. F., Andersson, J., Auerbach, E. J., Bijsterbosch, J., Douaud, G., et al. (2013). Resting-state fMRI in the human connectome project. NeuroImage 80, 144–168. doi: 10.1016/j.neuroimage.2013.05.039

PubMed Abstract | Crossref Full Text | Google Scholar

Smith, S. M., Miller, K. L., Salimi-Khorshidi, G., Webster, M., Beckmann, C. F., Nichols, T. E., et al. (2011). Network modelling methods for FMRI. NeuroImage 54, 875–891. doi: 10.1016/j.neuroimage.2010.08.063

Crossref Full Text | Google Scholar

Stam, C. J., Nolte, G., and Daffertshofer, A. (2007). Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum. Brain Mapp. 28, 1178–1193. doi: 10.1002/hbm.20346

PubMed Abstract | Crossref Full Text | Google Scholar

Stam, C. J., and Van Dijk, B. W. (2002). Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets. Phys. Nonlinear Phenom. 163, 236–251. doi: 10.1016/S0167-2789(01)00386-4

Crossref Full Text | Google Scholar

Steuer, R., Kurths, J., Daub, C. O., Weise, J., and Selbig, J. (2002). The mutual information: detecting and evaluating dependencies between variables. Bioinformatics 18, S231–S240. doi: 10.1093/bioinformatics/18.suppl_2.S231

PubMed Abstract | Crossref Full Text | Google Scholar

Stylianou, O., Kaposzta, Z., Czoch, A., Stefanovski, L., Yabluchanskiy, A., Racz, F. S., et al. (2022). Scale-free functional brain networks exhibit increased connectivity, are more integrated and less segregated in patients with Parkinson’s disease following dopaminergic treatment. Fractal Fract. 6:737. doi: 10.3390/fractalfract6120737

PubMed Abstract | Crossref Full Text | Google Scholar

Stylianou, O., Racz, F. S., Eke, A., and Mukli, P. (2021a). Scale-free coupled dynamics in brain networks captured by bivariate focus-based multifractal analysis. Front. Physiol. 11, 1–14. doi: 10.3389/fphys.2020.615961

PubMed Abstract | Crossref Full Text | Google Scholar

Stylianou, O., Racz, F. S., Kim, K., Kaposzta, Z., Czoch, A., Yabluchanskiy, A., et al. (2021b). Multifractal functional connectivity analysis of electroencephalogram reveals reorganization of brain networks in a visual pattern recognition paradigm. Front. Hum. Neurosci. 15:740225. doi: 10.3389/fnhum.2021.740225

PubMed Abstract | Crossref Full Text | Google Scholar

Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D., and Leahy, R. M. (2011). Brainstorm: a user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci. 2011, 1–13. doi: 10.1155/2011/879716

PubMed Abstract | Crossref Full Text | Google Scholar

Tadel, F., Bock, E., Niso, G., Mosher, J. C., Cousineau, M., Pantazis, D., et al. (2019). MEG/EEG group analysis with brainstorm. Front. Neurosci. 13:76. doi: 10.3389/fnins.2019.00076

PubMed Abstract | Crossref Full Text | Google Scholar

Taulu, S., and Hari, R. (2009). Removal of magnetoencephalographic artifacts with temporal signal-space separation: Demonstration with single-trial auditory-evoked responses. Hum. Brain Mapp. 30, 1524–1534. doi: 10.1002/hbm.20627

Crossref Full Text | Google Scholar

van den Heuvel, M. P., and Fornito, A. (2014). Brain networks in schizophrenia. Neuropsychol. Rev. 24, 32–48. doi: 10.1007/s11065-014-9248-7

Crossref Full Text | Google Scholar

Zebende, G. F. (2011). DCCA cross-correlation coefficient: quantifying level of cross-correlation. Phys. Stat. Mech. Its Appl. 390, 614–618. doi: 10.1016/j.physa.2010.10.022

Crossref Full Text | Google Scholar

Zhan, L., Jenkins, L. M., Wolfson, O. E., GadElkarim, J. J., Nocito, K., Thompson, P. M., et al. (2017). The significance of negative correlations in brain connectivity. J. Comp. Neurol. 525, 3251–3265. doi: 10.1002/cne.24274

Crossref Full Text | Google Scholar

Zhang, N., Lin, A., and Yang, P. (2020). Detrended moving average partial cross-correlation analysis on financial time series. Phys. Stat. Mech. Its Appl. 542:122960. doi: 10.1016/j.physa.2019.122960

Crossref Full Text | Google Scholar

Appendix

Auto-regressive fractionally integrated moving-average processes

Assume two distributions ${\epsilon }_A$ and ${\epsilon }_{{\rm B}}$ . ${\epsilon }_A$ is a standard normal distribution, meaning E[ ${\epsilon }_A$ ] = 0 and var( ${\epsilon }_A$ ) = 1. ${\epsilon }_{{\rm B}}=\rho {\epsilon }_A+\epsilon \sqrt{1-{\rho }^2}$ , where ε is also standard normal [i.e. E[ε] = 0 and var(ε)=1]. The variance of ${\epsilon }_{{\rm B}}$ can be calculated as follows:

var( ${\epsilon }_B)=var\left(\rho {\epsilon }_A+\epsilon \sqrt{1-{\rho }^2}\right)=var\left(\rho {\epsilon }_A\right)+var\left(\epsilon \sqrt{1-{\rho }^2}\right)={\rho }^{2}var\left({\epsilon }_A\right)+\left(1-{\rho }^2\right)var\left(\epsilon \right)={\rho }^2+1-{\rho }^2=1$

Then the real coupling between the two distributions can be calculated as:

ρ( ${\epsilon }_A,{\epsilon }_B)$ = $\begin{array}{l}\frac{cov\left({\epsilon }_A,{\epsilon }_B\right)}{\sqrt{var\left({\epsilon }_A\right)var\left({\epsilon }_B\right)}}=cov\left({\epsilon }_A,{\epsilon }_B\right)=E\left[{\epsilon }_A{\epsilon }_B\right]-E\left[{\epsilon }_A\right]E\left[{\epsilon }_B\right]=E\left[{\epsilon }_A{\epsilon }_B\right]=E\left[\rho {{\epsilon }_A}^2+{\epsilon }_A\epsilon \sqrt{1-{\rho }^2}\right]=E\left[\rho {{\epsilon }_A}^2\right]+E\left[{\epsilon }_A\epsilon \sqrt{1-{\rho }^2}\right]\\ =\rho E\left[{{\epsilon }_A}^2\right]+\sqrt{1-{\rho }^2}E\left[{\epsilon }_A\epsilon \right]=\rho E\left[{{\epsilon }_A}^2\right]+\sqrt{1-{\rho }^2}E\left[{\epsilon }_A\right]E\left[\epsilon \right]=\rho E\left[{{\epsilon }_A}^2\right]=\rho \left(var\left({\epsilon }_A\right)+E{\left[{\epsilon }_A\right]}^2\right)=\rho \end{array}$

The two ARFIMA series ( $A={\displaystyle \sum }_{n=0}^{100}{\alpha }_n\left(d\right){\epsilon }_{A,t-n},B={\displaystyle \sum }_{n=0}^{100}{\alpha }_n\left(d\right){\epsilon }_{B,t-n}$ ) are the cumulative sums of ${\epsilon }_A$ and ${\epsilon }_{{\rm B}}$ multiplied by a step-specific weight [( ${\alpha }_n(d$ )]. The only source of stochasticity of A and B are ${\epsilon }_A$ and ${\epsilon }_{{\rm B}}$ , meaning that the true coupling between A and B is ρ.