Model-based estimation of dynamic functional connectivity in resting-state functional magnetic resonance imaging

Recently, we have witnessed an increase in scientific interest in understanding the dynamic nature of brain networks by evaluating dynamic functional connectivity (FC) using functional magnetic resonance imaging (fMRI). In this work, we introduce two multivariate volatility models, standardized dynamic conditional correlation, and standardized exponentially weighted moving average, both of which are built upon the framework of dynamic conditional correlation and exponentially weighted moving average models, respectively. In these two models, we use standardized residuals with the goal of determining whether the use of standardized residuals reduces the mean square rate error. Moreover, in traditional simulation studies, time series were considered with zero conditional expectation and static conditional variance which do not capture the nature of the real data. This is because of hemodynamic response function in the brain and dynamic functional connectivity of each brain region with itself during the experiment time, respectively. That is why, next, some new simulation studies are introduced which are more similar to blood-oxygen-level-dependent time series of brain regions. Then, methods’ proficiency is analyzed using previous and new simulation studies. Results show that, in both former and latter simulations, the new methods work better. Finally, the best model is utilized to model FC in an Iranian resting-state fMRI data.


Introduction
One of the fields of interest in functional magnetic resonance imaging (fMRI) is finding the functional connectivity (FC) between brain areas. In this field, the dependency of blood-oxygen-level-dependent (BOLD) response time series from various brain regions is studied.
Heretofore, these time series had been assumed to be stationary over time. Therefore, FC was considered to be constant across time (static FC) [1][2][3]. In recent years, scientists have observed indications of non-stationarity in these time series both in resting-state and in task-based fMRI experiments [4][5][6][7][8][9][10][11]. This evidence led to the introduction of a new concept in evaluating brain regions connectivity named dynamic FC which has encouraged the introduction of a large number of methods to detect it. Sliding window, for instance, has been a widely used method to assess dynamic FC [12][13][14][15][16]. Dynamic connectivity regression (DCR) is another method recommended by [4,5] to discover FC change points between brain areas. [17] proposed dynamic connectivity detection (DCD) algorithm based on DCR algorithm. This new method is faster and has no problem finding changes in connectivity for more than 100 brain regions. In this regard, network change point detection (NCPD) method was introduced by [18]. Whole-brain dynamics can be investigated through this method. As other studies in this field, [19] suggested dynamic Bayesian variable partition model (DBVPM) to model dynamic FC via a unified Bayesian framework. Product hidden Markov models (PHMMs) were recommended by [20] to assess temporal behavior of dynamic FC in resting-state fMRI. [21] pointed out the drawbacks of the sliding window method drawbacks and recommended two model-based approaches named exponentially weighted moving average (EWMA) and dynamic conditional correlation (DCC) models and demonstrated that DCC performed better than sliding window by calculating mean square error (MSE) in simulation studies. These two well-known models, which simulate dynamic correlation in finance literature, have recently drawn a lot of attention from scientists in the field of neuroscience, like real-time fMRI. For instance, [22] proposed to utilize EWMA models as well as Smooth Incremental Graphical Lasso Estimation (SINGLE) algorithm to estimate dynamic FC in real time.
Functional connectivity is evaluated by estimating the covariance matrix of two BOLD response time series. Therefore, non-stationarity in BOLD response time series leads to time-variant covariance matrix (dynamic FC). When the conditional covariance is time-variant, we have conditional heteroscedasticity which is the main focus of multivariate volatility modeling. The effective estimations derived from volatility models in neuroimaging context [21,22] motivated us to investigate and modify them. In this regard, the EWMA and DCC models introduced in [21] are further developed in this article. To be more precise, we present two new multivariate volatility models built upon the framework of DCC and EWMA models and named standardized dynamic conditional correlation (SDCC) and standardized exponentially weighted moving average (SEWMA). These two new models apply standardized residuals instead of residuals. Our goal is to examine whether exploiting standardized residuals can reduce the mean square error. This idea arose from this fact that the residuals have fat tails. Furthermore, some new simulation studies are proposed based on the simulation designed by [21]. These simulations are designed in such a way that they more closely mirror the blood-oxygen-level-dependent (BOLD) time series of brain areas. Although former simulations were well defined, they had some shortcomings. The first one is that, in all simulations, the conditional expectation was equal to zero. This condition is not realistic due to the existence of hemodynamic response function (hrf) in the brain. Furthermore, each brain region has a dynamic association with itself during the time. Thus, considering static variances in those simulations is not acceptable. The performance of these two new methods is surveyed via previous and recent simulation studies. After calculating the MSE in these simulation studies, it is proved that the new methods work better than the older ones in both previous and new simulations. Finally, the best model fits to an Iranian restingstate fMRI data.
This article is organized as follows: Two new models are introduced in Sect. ''Methods.'' First, we describe basic concepts required in estimating FC. Next, SEWMA and SDCC models are nominated. Section ''Simulation studies'' explains simulation studies and demonstrates their results. Experimental data are illustrated in Sect. ''Experimental data.'' In this part, we first state information of the Iranian resting-state fMRI data and brain regions. After that, the best model fits to the real data in order to estimate the FC. Finally, the presented results are discussed in the last section.

Methods
In this part, we first explain BOLD response time series and hrf in preliminaries. Then, the problem of evaluating dynamic FC is illustrated in problem setup.

Preliminaries
The hrf is the hypothetical BOLD response to an idealized impulse of neural activation [23]. In fact, where B t is the BOLD response time series, N t is the neural activation, h t is the hrf, and indicates convolution. This fact is clearly demonstrated by Fig. 1. It should be mentioned that resting-state fMRI is studied in this article. Thus, there is not any information about neural activity onsets.

Problem setup
Let y t ¼ ðy 1;t ; y 2;t Þ T , t ¼ 0; . . .; T, is a vector where y i;t , i ¼ 1; 2, is a time series extracted from one region of interest (ROI) in the brain. Suppose that with and where F i;tÀ1 indicates the information set available up to time t -1. Moreover, e t is the noise term with zero mean and conditional covariance matrix R t as follows: where and Residual is also defined as below: wherel t is the fitted conditional mean of y t . In fact, residual is the difference between the observed value and the predicted value. The conditional covariance matrix R t can be written as follows: where D t ¼ diagfr 1;t ; r 2;t g; ð11Þ To evaluate the connectivity of two separate brain regions, we should estimate R t and, therefore, r 1;t ; r 2;t , and q t . [24] recommended two multivariate volatility models, EWMA and DCC models, to estimate R t in finance literature, and [21] advocated these two models to estimate dynamic FC in resting-state fMRI data. In this article, we have attempted to improve both of them as below. Unlike [21], we put l t 6 ¼ 0 in Eq. (2) since in actual fMRI data, the conditional mean (l t Þ is not zero. This fact is due to the hrf in the brain.

Standardized exponentially weighted moving average (SEWMA)
This method, like EWMA, places emphasis on the basic fact that past observations are less relevant than the recent ones. The SEWMA model is where a Ã t ¼ ða Ã 1;t ; a Ã 2;t Þ is a vector of standardized residual, and a Ã i;t ; i ¼ 1; 2; is defined in the following way where r 2 i;t is the conditional variance of y i;t and is defined in Eq. (6). Furthermore, 0\k\1 indicates the decaying rate. If y t is bivariate Gaussian, then we can estimate k via maximum likelihood estimation. In this survey, we use k ¼ 0:94 as a default value. This value was recommended in RiskMetrics [25], and not only it is common in financial literature [26], but also it is proved to be appropriate for fMRI data [21].

Standardized dynamic conditional correlation model (SDCC)
This model is based on DCC model introduced by [24] which similarly has two steps. Suppose we have y t = l t þ e t with the mentioned assumptions, we have a Ã i;t which is defined in Eq. (14). Now, in the first step, r 2 i;t , i ¼ 1; 2 is obtained as below: And in the second step, R t is modeled as follows: where Q is the unconditional covariance matrix of a Ã t and h 1 ; h 2 ð Þ are nonnegative real numbers so that 0\h 1 þ h 2 \1. w 1 ; n 1 ; b 1 ; w 2 ; n 2 ; b 2 ; h 1 ; h 2 are model parameters estimated via a two-step method.

Simulation studies
In order to assess different models performance, we not only designed some new simulation studies, but we also made use of the former simulations recommended by [21]. It should be mentioned that the new simulation studies are designed with respect to the former simulations designed by [21]. In this part, we want to generate random data y t ¼ ðy 1;t ; y 2;t Þ T , t ¼ 0; . . .; T. As mentioned in Eq. (2), we have y t ¼ l t þ e t . In the previous simulation studies, the conditional expectation was equal to zero. So, Eq. (2) is reduced to where e t is the noise term with zero mean and conditional covariance matrix R t as mentioned in Eq. (5). Therefore, [21] generated y t using a mean-zero bivariate normal distribution with covariance matrix where r 12;t is the covariance term and varies across time. In fact, variances are constant (time-invariant) and only covariance is time-variant. In simulation studies, our goal is to control the relationship between y 1;t ; y 2;t by q t as below: 1. r 12;t ¼ 0; for t ¼ 0; . . .; 600; 2. r 12;t ¼ sin t 2 k À Á with k ¼ 6; 7; 8; and 9 and t ¼ 0; . . .; 600; 3. r 12;t ¼ exp À xÀ250 2r 2 À Á with r ¼ 15 Ã k; k ¼ 1; 2; 3; and 4; and t ¼ 0; . . .; 600: Each simulation study was repeated 1000 times, and the MSE of four models (EWMA, SEWMA, DCC, and SDCC) were measured.
In order to compare the results, we compute Mean Diff1 and Mean Diff2 as below: where MSE A i ð Þ is the mean square error of model A in iteration i. Positive values for Mean Diff1 and Mean Diff2 mean that SEWMA works better than EWMA and SDCC outperforms DCC in that simulation, respectively. Table 1 abstracts the results of Mean Diff1 and Mean Diff2 for each simulation.
As it is seen in Table 1, the values of Mean Diff1 and Mean Diff2 are positive. This shows that in the former simulations, the proposed models work better.
Since these values are small, it should be tested whether the models' results are significantly different. In this regard, we first used a one-sample Kolmogorov-Smirnov test for MSE EWMA , MSE SEWMA , MSE DCC , and MSE SDCC data sets. In 5% significance level, the test rejects the null hypothesis that each data set follows a standard normal distribution. Therefore, a Mann-Whitney-Wilcoxon test was utilized to compare the means of the groups as it is demonstrated in Table 2.     Table 1 are positive. Therefore, it can be concluded that, in the former simulation studies, SEWMA and SDCC outperform EWMA and DCC models, respectively.
In this work, one of our goals is to improve the former simulation studies. As mentioned before, we should not put l t ¼ 0, because the changes in fMRI signals are trigged by instantaneous neural activity which is known as the hrf. So, according to Eq. (2) (y t = l t þ e t ), we should opt for a realistic value for l t . The best choice is to utilize the conditional expectation of real BOLD response time series which we gain from fMRI data. As we assess pair-wise FC, the BOLD response time series of two brain regions are needed. To do this, we first randomly selected one person's fMRI data and extracted the BOLD response time series of two brain areas. One of the two selected regions was always posterior cingulate cortex (PCC), and the other one was randomly selected from the remaining five regions. In experimental data, we explain brain regions in detail. Next, conditional expectations, l 1;t and l 2;t , were conditional mean of those BOLD response time series. In addition, e t should be generated by a mean-zero bivariate normal distribution with the following covariance matrix where r 2 i;t ; i ¼ 1; 2, is conditional variance of those BOLD response time series which were extracted. Unlike the previous simulation studies, we have chosen time-variant variances since each brain region has dynamic FC with itself. Moreover, r 12;t was determined in three different forms: Then, for each of these three forms, the conditional correlation q t is computed as Finally, y t is achieved by substituting the values of l t and e t in Eq. (2) (y t = l t þ e t Þ. Each of these three simulations was repeated 1000 times. In each repetition, EWMA, DCC, SEWMA, and SDCC models were fitted to the simulated data and the MSEs were computed. Figure 2 demonstrates the result for simulation 1. In this simulation, we suppose that two brain regions do not have any connectivity during the experiment. For this simulation, subject 10 and ROI 1 were both randomly selected. Part (a) shows the result of fitting EWMA, DCC, SEWMA, and SDCC models in one of the iterations, and part (b) shows the box plot for MSE of these three models. According to the part (b), SEWMA and SDCC perform better than EWMA and DCC, respectively. And judging by the results, it is evident that SDCC is the best model among them. Figure 3 displays the results for the second simulation. In this simulation, we suppose that connectivity of two brain regions is in the form of sine. The information (number) about selected subjects and regions is written on top of the figures. Part (a) shows the results of fitting EWMA, DCC, SEWMA, and SDCC models in one of the iterations, and part (b) shows the box plot for MSE of these four models. According to the part (b), SEWMA and SDCC perform better than EWMA and DCC, respectively. The results for the third simulation are exhibited in Fig. 4. In this simulation study, we suppose that connectivity of two brain regions is in the Gaussian form. The information (number) about selected subjects and regions is written on top of each figure. Part (a) shows the result of fitting EWMA, DCC, SEWMA, and SDCC models in one of the iterations, and part (b) shows the box plot for MSE of these four models. According to the part (b), SEWMA and SDCC outperform EWMA and DCC, respectively. Having the lowest MSEs, SDCC is clearly the best model.

Experimental data
These study protocol and consent form were verified by the Ethics Committee of Tehran University of Medical Sciences. Before scanning, the imaging procedure was thoroughly explained for all subjects and their written informed consents were acquired. This study contained 20 healthy male volunteers (23-46 years old). Each resting-state scan took 7 min. A 2D EPI sequence with partial parallel imaging acceleration was utilized to gain 3 9 3 mm in-plane resolution in 3-mm axial slices with a 1-mm slice gap. The slice order was interleaved with TR/TE = 2500/30 ms, and flip angle of 50 and SENSE acceleration factor of 2 were applied.
The preprocessing step and extracting brain regions time series were done via SPM 12, DPARSF 4.3, and MATLAB codes. The preprocessing steps included these steps: slice time correction, realignment, segmentation, normalization, smoothing (6-mm kernel), and reducing physiological noise by CompCor [27].
Brain regions of interest (ROIs) were selected based on [13]. In the default mode network (DMN), the posterior cingulate cortex (3-mm-radius sphere centered at x = -6, y = -58, z = 28) was selected as the principal region and five other regions of interest (ROIs), demonstrated by [13] as brain regions with the most enormous variation in dynamic correlation with PCC, were also chosen. These brain regions are as follows: ROI1 Right inferior parietal cortex (3-mm-radius sphere centered at x = 38, y = -58, z = 44).
SDCC model was selected to assess dynamic FC in this study. Figure 5 represents BOLD response time series of PCC and the other ROIs, and the estimation of their dynamic FC. The information (number) about selected subjects and regions is written on top of each figure. Moreover, 95% confidence intervals and static FCs can be seen in this figure. Results show that in these three subjects, the dynamic FCs do not change a lot, and there are few points in static FCs which are outside of the confidence intervals.

Discussion
In this paper, we try to estimate dynamic FC, one of the research priorities in the world, by using some models with lower error. Our criterion for choosing the best model is to calculate the models MSE in simulation studies which are close to BOLD response time series in brain regions.
EWMA and DCC models were introduced by [24] for the first time in finance literature, and [21] made use of these two models to assess dynamic FC in the brain. In this work, we present SDCC and SEWMA models based on DCC and EWMA, respectively, with the aim of decreasing MSE. The basic difference between the former and latter models is applying the standardized residual, a Ã i;t , for i ¼ 1; 2, in the new models. This idea is due to the fact that the residuals have fat tails. Although missing rapid changes is the common shortcoming of SDCC and DCC models, they have lots of merits which cause us to prefer them over the other models. Dynamic connectivity regression (DCR) is an appropriate method for the situations with rapid changes [4,5,21]. SDCC models inherit the advantages of DCC models. For instance, random noise does not disturb them to recognize the real correlation, and models' parameters are estimated via a quasi-maximum likelihood method. Moreover, SDCC models have lower MSE than DCC models.
Furthermore, we strive to apply more realistic simulation studies. In this survey, we have l t 6 ¼ 0. This choice is due to the hrf in the brain. Hemodynamic response is because of any kind of stimulus that brain receives even in resting state. We also suppose that each time series variance is time-variant. This assumption is owing to the dynamic association of each brain region with itself during the experiment time.
To check these models, we compared these models MSE in both previous and recent simulation studies. The results showed that both new models performed better than the older ones in both former and latter simulation studies. Although, in the previous simulations, the difference between the former and proposed models was less than these values in the new simulations, it was shown that these little differences are significant in the 5% significance level. The former simulations consider the ideal generated data with l t ¼ 0 and mean-zero bivariate normal distribution for e t . So, standardized residuals and residuals have fewer differences in those generated data than in ones generated in new simulation studies.
SDCC model outperformed the other methods in these studies. Hence, we exploited this model to find dynamic FC in the Iranian resting-state fMRI data. For these 20 subjects, we calculated dynamic FC between PCC and five other regions which according to [13] had high variations in their correlation with PCC.
Finally, we have reached the conclusion that SDCC models outperform the other models even DCC models, and they are a suitable option for estimating dynamic FC in resting-state fMRI. We have witnessed differences in b Fig. 5 Shown are BOLD response time series of PCC and the other ROIs, and the estimation of their dynamic FC. The confidence intervals and static FCs are also presented estimated dynamic FC between distinctive regions and also different subjects. This outcome only seems reasonable since not only do we have rest fMRI which is task-free, but also different people do not have exactly the same results in real experiments.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.