Spontaneous brain activity measured by functional magnetic resonance imaging (fMRI) has provided the evidence that the human brain is intrinsically organized into large-scale functional networks. Functional magnetic resonance imaging can provide evidence for functional specialization and integration via its ability to simultaneously localize neural activity in the entire brain. Friston et al. (1993) defined functional connectivity (FC) as the temporal correlation between spatially remote neurophysiological events. FC is one way to characterize such interactions. Understanding functional connectivity within the brain is crucial to understanding neural function, which supports highly distributed neural circuits.

Temporally varying information provides insight into the fundamental properties of brain networks. However, most fMRI studies have assumed that the essential connectivity between time series from distinct brain regions is constant. Recently there have been increasing attempts to quantify the dynamic changes in FC, such as Lindquist et al. (2014).

Though it is of increasing importance, interpreting temporal fluctuations in FC is difficult due to low signal-to-noise ratio, physiological artifacts, and variation in BOLD signal mean and variance over time (Hutchison et al., 2013). For these reasons, it is often difficult to determine whether observed fluctuations in FC should be attributed to neuronal activity or whether they are due to random noise, and thus research is still ongoing in this area.

Most current approaches to examining functional connectivity (FC) implicitly assume that relationships are constant throughout the length of the recordings. However, investigations of intrinsic brain organization based on resting-state fMRI provides the presence and potential of temporal variability. Temporal trends in the occurrence of different FC states motivate theories regarding their functional roles and relationships with vigilance/arousal.

A variety of other approaches to identify FC states are also possible, including using topological descriptions of brain connectivity as features, e.g., modularity or community membership (Bassett et al., 2011; Jones et al., 2012; Kinnison et al., 2012). Lindquist et al. (2007) detected dynamic changes in BOLD response using EWMA related approach. Formal models for detecting change points in connectivity are introduced by Hampson et al. (2002), Robinson et al. (2010). Cribben et al. (2012) and Aston and Kirch (2012) proposed other promising approaches to improve methods for identifying FC states and state transitions. Huang et al. (2019) proposed windowless dynamic correlation using a heat kernel. Recently, phase synchronization (PS) methods were proposed as a measuring method of the level of synchrony between time series from two different regions of interest (ROI)s in the brain (Glerean et al., 2012; Pedersen et al., 2018). Honari et al. (2021) suggested assessing the relationship between signals from different brain regions and measuring their phase synchronization (PS) across time. Kim et al. (2021) proposed the change-point estimation based on random matrix theory for the epilepsy fMRI. They utilized a PS metric and a sliding window for instantaneous phase synchronization (IPS).

Dynamic functional connectivity (DFC) uses the correlations with the fMRI time series according to the time-points. This temporally varying information may be used to provide insight into the fundamental properties of brain networks. Here the correlation measurement method is an issue in studying functional connectivity. It is often difficult to determine whether observed fluctuations in FC should be due to neuronal activity or random noise. Thus significant research is still in need in the area. In particular, more research should be needed to add the appropriate analysis strategy.

This paper focuses on pairwise correlations between time courses from two brain regions. We consider the most commonly used sliding-window approach (Chang and Glover, 2010; Handwerket et al., 2012), EWMA, as proposed in Lindquist et al. (2007), the recently suggested volatility model-based DCC approach in Lindquist et al. (2014), and VAR model-based dynamic causality approach. In particular, we try the method with the difference in time series and a moving model with recursive residuals. Overall, we suggest that the study of time-varying aspects of FC can unveil flexibility in the operational coordination between different neural systems. The exploitation of these dynamics and the corresponding FC estimation method may improve our understanding of behavioral shifts and adaptive processes.

This paper is organized as follows. Section 2 provides pairwise correlation estimation methods. Section 3 gives simulation results, and the final discussions are in Section 4.

In this section, we set up the problem for pairwise correlation and introduce its estimation methods for dynamic connectivity.

The bivariate time series y_1t and y_2t, measured over two separate ROIs in the brain, at equally spaced time-points t = 1, . . ., T. Let y_t = (y_1t, y_2t)′ be a vector containing the values of both time series at time-point t. Assume that

where μ_t = (μ_1t, μ_2t)′ is the conditional mean of y_t using all information in the time series observed up to time t, denoted E_t(y_t). The noise term e_t has mean zero and its conditional covariance matrix at time t as

The diagonal terms in (2.2) represents the conditional variance of y_i,t using all the information in the time course up to time t. σi,t2 is typically referred to as the volatility of the time series in the finance literature. The off-diagonal term is σ_12,t = ρ_tσ_1,tσ_2,t where

represents the conditional correlation coefficient. This definition shows that the conditional correlation at time t relies on information that is obtained up to time t – 1. The quantity lies in [−1, 1].

Without loss of generality we assume that μ_t = 0 and therefore y_t = e_t. The conditional covariance matrix in (2.2) can be written as

where D_t is a diagonal matrix such as D_t = diag{σ_1,t, σ_2,t} and R_t is the correlation matrix as

Our primary interest is in developing estimation of the conditional covariance whose components are σ_1,t, σ_2,t and ρ_t. In the following section, we discuss several methods for our comparison.

2.1. Sliding-window approach

The simplest approach to estimating the elements of the covariance matrix is to use the sliding-window technique. A time window of fixed length w is selected, and data points within the window are used to calculate the correlation coefficients. The window is moved across time, and a new correlation coefficient is computed for each time-point.

Chang and Glover (2010) define the sliding-window correlation at time t at the resting-state as follows:

We consider the sliding-window correlation based only past values as a more suitable estimate of the conditional correlation as

The sliding-window technique allows for a simple approach for exploring changes in connectivity. However, it has some obvious shortcomings. First, it gives equal weight to all observations less than w time points in the past and 0 importance to all others. Highly influential outlying data points will cause a sudden change in the dynamic correlation and may be mistakenly a critical aspect of brain connectivity. To circumvent this problem, Allen et al. (2014) suggested using a tapered sliding window.

2.2. EWMA approach

The EWMA (exponentially weighted moving average) (Hunter, 1986) approach applies declining weights to the past observations and the most weight on recent observations in the time series based on a parameter λ as follows:

where ∑_t is the conditional covariance matrix.

Decomposing the covariance matrix (2.8), the conditional variances and covariances are expressed in the component as follows:

and

With recursive computaion, a small value of λ gives high weight to the recent time-points, while a large value produces the estimates more gradually to new information. The parameter λ takes values between 0 and 1, This value determines how many data points are included in the correlation calculation and serves equivalently to the window size used in the sliding-window approach.

If one assumes that y_t are bivariate normal, the optimal λ can be chosen through maximum likelihood estimation as:

2.3. DCC method

Lindquist et al. (2014) introduced DCC (dynamic conditional correlation) as an approach to estimating the conditional variances and correlations that has become increasingly popular in the finance literature. DCC incorporates the GARCH process (Bollerslev, 1986; Engle, 1982) used to model volatility in univariate time series, providing conditional variance in the following. Let y_t = e_t again then with GARCH expression we have

Define

where D_t = diag{σ_1,t, σ_2,t}. Let

where Q̄ is the unconditional covariance matrix. Decompose the covariance matrix as

where Λt=diag{Q}t-1/2Qtdiag{Q}t-1/2. Here θ₁ and θ₂ are non-negative scalars satisfying 0 < θ₁+θ₂< 1. The estimation procedure of θ₁ and θ₂ is in Appendix A in Lindquist et al. (2014).

Wiener (1956) first discussed the causality between the variables in the observed multivariate time series. Granger (1969) studied Wiener’s idea and introduced a few concepts related to causality, mainly in the framework of bivariate AR modeling. Recently a similarity between the causality study in economics and neuroscience has been recognized.

2.4. VAR model

Goebel and Roebroeck (2003), and Roebroeck and Formisano (2005) have proposed the use of VAR (vector autoregressive model) and shown their utility in the analysis of fMRI experiments.

where u_t = (u_1,t, u_2,t)′ is a random error vector with mean zero and covariance matrix ∑ given by

This can be expressed with the autoregressive order p as

where

and u_t is an error vector of random variables with zero mean and covariance matrix ∑(t) given by

This Granger causality is shown helpful for inferring functional brain connectivity. However, VAR modeling is an adequate approach for stationary time series i.e. the autoregressive coefficients and error matrix covariance are time-invariant.

We let VARw as the VAR method with the sliding window to localize the correlation calculation. Within the observed data in the window, we apply the VAR method locally.

2.5. Difference AR coefficient model

Consider the AR(1) structure for each time series since there usually exists autocorrelation at lag1 at least,

Estimate each AR coefficient as ϕ̂₁ and ϕ̂₂ and then compute the followings

After the autoregressive and linear trend is removed, the pairwise correlation with the window w is calculated such as

We also try the conditional correlation method with the sliding window as denoted DCCw to localize the correlation calculation. AR (autogressive) of order 1 and temporal linear trends are estimated for each observed series in the window and then their residuals are used in the pairwise correlation calculation. This allows observed points to enter and exit from the window as it moves across time removing the linear trend.

2.6. Linear moving model with recursive residuals

Starting with Brown et al. (1975), residuals became one of the most important tools in change-point analysis for testing the constancy of regression relationship. The recursive (rec) residual for l = 1, 2 is defined for each time series such as

where Xt-1′=[x1,…,xt-1], Yl,t-1′=[y1,1,…,yl,t-1] and β^t-1=(Xt-1′Xt-1)-1Xt-1′Yt-1 is the ordinary least-squares estimator. Here x_i = t_i is the same time-points for each times series. Assuming that ε_i’s follow the normal distribution, Brown et al. (1975) showed that under H₀, w_l,k+1, . . .,w_l,T are independent N(0, σ²).

These recursive residuals are used for correlation calculation as

Since the recursive residuals are obtained with removing the linear trend, the correlation with these recursive residuals may not be mixed with the linear trend.

In this section, we set up the problem situations for pairwise correlation and compare several of its estimation methods for dynamic connectivity. The two time courses are designed to be generated from the bivariate normal distributions. In each case, the covariance matrix of the distribution was set as a possible structure. The following cases help us understand the different underlying change patterns. We choose T = 300 for the total time-points. The simulated data are mean corrected for calculation. For change cases, we set up one change-point at T₁ = T/2, For two change cases, we consider change-points at T₁ = T/3 and T₂ = (2/3)T.

where the covariance term p(t) are allowed to vary across time for t = 1, . . ., T, which are allowed to control the dynamic relationship between the two time courses y_1,t and y_2,t. The number of time-points T = 300 is considered in the simulation. The following simulation setup is considered.

• Case 1. Uncorrelated and no change in correlations

  p(t)=0,         t=1,…,T.

• Case 2. Slowly varying smooth change in correlation.

  p(t)=0.5sin(t/Δ),         t=1,…,T, Δ=1024/2k, k=3 (chosen).

• Case 3. Uncorrelated but variances change.

  (3.2)Σt={Σ1=(1002),1≤t≤T1=T/2,Σ2=(2003),T1<t≤T.

• Case 4. p(t) is linear according to time.

  (3.3)p(t)=bt,         1≤t≤T,   b=0.3.

• Case 5. p(t) has one abrupt change-point in correlation (Case 1 + one change-point at T₁).

  (3.4)p(t)={0,1≤t≤T1,0.4,T1<t≤T.

• Case 6. p(t) is linear according to time with one change-point (Case 2 + one change-point).

  (3.5)p(t)={bt,1≤t≤T1=T/2,a+bt,T1<t≤T,

  where a = 1/Δ, Δ = 500.

• Case 7. One change-point in correlation with variance change (Case 3 + one change-point).

  (3.6)p(t)={0,1≤t≤T1,0.4,T1<t≤T,

  and

  (3.7)Σt={Σ1=(1002),1≤t≤T1,Σ2=(2003),T1<t≤T.

• Case 8. p(t) is linear according to time (Case 4 + one change-point).

  (3.8)p(t)={bt,1≤t≤T1,a+bt,T1<t≤T,

  where a = 0.1, b = 0.2.

• Case 9. Two change-points (Case 1 + two change-points at T₁ and T₂).

  (3.9)p(t)={0,1≤t≤T1,0.3,T1<t≤T2,0,T2<t≤T.

• Case 10. Two change-points in periodic correlation (Case 2 + two change-points).

  (3.10)p(t)={0.1+0.3 sin(t/Δ),1≤t≤T1,sin(t/Δ),T1<t≤T2,0.1+0.3 sin(t/Δ),T2<t≤T

  with Δ = 1024/2^k, k = 3 (in the simulation).

The following abbreviation for methods are used in the comparison results:

• Sl: slide window estimation with the bandwidth h = 0.1 that is, w = 30.

• EWMA: exponentially weighted moving average.

• DCC: dynamic conditional correlation.

• DCCw: driving conditional correlation based on AR(1) and linear trend removed in the window w.

• VAR: vector autoregressive model.

• VARw: vector autoregressive model with the moving window w.

• Dc: autoregrssive difference corrected model, advantage with AR(1) model.

• LMc: correlation with linear moving recursive residuals.

The mean square error (MSE) is computed by calculating the mean of the squared difference between the estimated correlation and the true correlation at time-points. Table 1 provides the mean, median, and SD (standard deviation) of MSE of the corresponding method in each case. We get the results using Sl, EWMA, DCC, DCCw, VAR, VARw, Dc, and LMc models in 1,000 repetitions with T = 300. We use the window size w = 30 as fixed for sliding window and moving average methods. For EWMA λ = 0.3 is chosen to give more weight to the recent terms. There is some potential to slightly improve the performance based on the optimal widow size. In each case, two figures show the time-varying correlation estimates according to time with the boxplot of MSE’s. The left panel of Figures shows an example of the true correlation plotted together with the estimated correlations obtained by the considered methods for a single simulation repetition. The smallest MSE means are highlighted in each case in Table 1. The best method varies according to case situation. The right panel shows boxplots of the MSEs for all 1,000 repetitions in the simulation. In Figures 1–4 VAR and Dc have smaller mean of MSEs compared to others. VARw and DCCw show similar behavior since they both use AR coefficients. DCC works better with smooth correlation changes in Figure 2. In Figures 5–8 for one abrupt change-point in correlation function, the best performance of the methods depends on the underlying change pattern. When there are two abrupt change-points in Figures 9–10, DCC, VAR, and Dc show similar performance. Since it incorporates multivariate relations, DCC and VAR show good performance generally. The sliding window methods depend on the window size that accounts for the underlying dynamic correlation patterns. We can choose the optimal window size to minimize some criteria. However, there is an inherent difficulty in determining an appropriate window length in a real-world situation. In our simulation, the performance of Dc would improve upon DCC and VAR for the autoregressive situation. Even though DCC is proposed and preferred by Lindquist et al. (2014), the estimation should be improved when there are change-points in correlations. We should consider change-point analysis previously beforehand.

Our simulations show that the performance of each method varies according to the underlying functional form. Removing the underlying trend and considering multivariate relatedness, then estimating the bivariate correlation is a principle. We expect that there is a way to circumvent these problems. The model’s complexity and data variability affect the decision of the correlation estimation method. As a result, we suggest the practical estimation of dynamic changes and correlations. In contrast to commonly used sliding windows techniques, DCC and VAR methods are captivating options.

We consider resting-state fMRI data from 20 patients with major depressive disorder (MDD) participating in a depression treatment biomarker study. Mayberg et al. (1999) proposed that the bidirectional nature of this limbic-cortical reciprocity provides additional evidence of potential mechanisms mediating cognitive, pharmacological, and surgical treatments of mood disorders such as depression. Kemmer et al. (2017) estimated the FC networks from 20 healthy subjects and 20 patients with major depressive disorder (MDD). The MDD patients were matched with the healthy control subjects by age and gender, with an average age of 45.8 years and 50% male. The mean Hamilton depression rating scale (HAM-D) score was 19 (SD: 3.4) for the MDD subjects. This score value indicates severe depression (Hamilton, 1960). The average length of their current depression episode was 82 weeks. Functional images were collected over 150 time-points on 14 ROIs. For each subject, 150 fMRI volumes were obtained in 7.5 minutes during a resting-state, in which subjects are left to think for themselves while visually focusing on a cross displayed on a monitor in the scanner. The following data preprocessing steps are acquired, including motion correction to offset subject movement in the scanner, slice timing correction to adjust for the fact that each 3D scan represents a single time point as a series of 2D slices over time, normalization to adopt a common reference space across subjects and spatial smoothing. They focus on the interactions between some of those reported ROIs (based on the Brodmann atlas) in the prefrontal cortex (PFC) including medial frontal cortex (mF10), orbital frontal cortex (OF11), and lateral prefrontal cortex (latF9); subgenual cingulate cortex (Cg25) and anterior cingulate (Cg24) shown in Figure 11.

For each subject, we have (142)=91 pairwise correlations. For example, we apply the correlation methods for fMRI data with MDD patients. There has been a growing interest in estimating potential dynamic correlations between two brain ROIs, as it is thought to provide important information about the properties of brain networks. Figure 12 shows the fMRI data on 14 ROIs for one MDD patient (MDD001) and one healthy control (CON021). For the pairwise correlations considered in this paper, we have Figure 14 which shows difference according to the method and some different patterns due to groups. Figures 15 and 16 show the bivariate correlations between ROI1 and other ROIs by the considered methods. In this application with the window size w = 30 out of T = 150. These examplary figures give a hint that the dynamic correlations should be carefully dealt in fMRI analysis. We observed a fair amount of variability in the behavior of these correlations across time, regions and subject.

This study presents a preliminary analysis of dynamic correlations between two brain regions, providing vital information about the properties of brain networks. We use several pairwise correlation measures to analyze time-varying interactions between brain regions. This work illustrates that we should consider change situations for correlation computation. Commonly used sliding-window techniques may not be beneficial for tracking dynamic correlations.

DCC and VAR models are extensively used in the finance literature for modeling time-varying variances and correlations and are generally considered preferable to sliding-window type approaches (Bauwens et al., 2006). Both have adaptability even though there are change-points. While the sliding window estimator is non-parametric, the DCC and EWMA models are parametric. These parametric approaches are powerful if the model is reasonably accurate. Moving estimators, including differences, can also be considered to avoid parameter estimation in the model. The more flexible non-parametric methods have advantages when the underlying function is changeable.

The wealth of information provided by time-frequency connectivity analysis presents additional challenges for studying multiple subjects and spatial locations. One way of handling this information is to summarize the dynamic information along several potentially-relevant dimensions.

Modeling time-varying variances and correlations is preferable to sliding-window type approaches. When considering the N–variate case, a total of N(N + 1)/2 variances and covariances need to be estimated, some of which may be time-varying. In addition, the resulting covariance matrix must be positive definite at each time-point. When there are underlying changes, consider the estimation of change-points (Kim et al., 2021), for example. Finally, the practical estimation of dynamic correlations depends on accurately analyzing the underlying change function.

This research was supported by 2023 Duksung Women’s University Research Fund.