Improved dynamic connection detection power in estimated dynamic functional connectivity considering multivariate dependencies between brain regions

Abstract To estimate dynamic functional connectivity (dFC), the conventional method of sliding window correlation (SWC) suffers from poor performance of dynamic connection detection. This stems from the equal weighting of observations, suboptimal time scale, nonsparse output, and the fact that it is bivariate. To overcome these limitations, we exploited the kernel‐reweighted logistic regression (KELLER) algorithm, a method that is common in genetic studies, to estimate dFC in resting state functional magnetic resonance imaging (rs‐fMRI) data. KELLER can estimate dFC through estimating both spatial and temporal patterns of functional connectivity between brain regions. This paper compares the performance of the proposed KELLER method with current methods (SWC and tapered‐SWC (T‐SWC) with different window lengths) based on both simulated and real rs‐fMRI data. Estimated dFC networks were assessed for detecting dynamically connected brain region pairs with hypothesis testing. Simulation results revealed that KELLER can detect dynamic connections with a statistical power of 87.35% compared with 70.17% and 58.54% associated with T‐SWC (p‐value = .001) and SWC (p‐value <.001), respectively. Results of these different methods applied on real rs‐fMRI data were investigated for two aspects: calculating the similarity between identified mean dynamic pattern and identifying dynamic pattern in default mode network (DMN). In 68% of subjects, the results of T‐SWC with window length of 100 s, among different window lengths, demonstrated the highest similarity to those of KELLER. With regards to DMN, KELLER estimated previously reported dynamic connection pairs between dorsal and ventral DMN while SWC‐based method was unable to detect these dynamic connections.

While dFC studies have recently drawn increasing attention, statistical assessment of the results to capture the underlying dynamic pattern from rs-fMRI data is of great importance. The spurious fluctuations due to inherent noise in the rs-fMRI data, low signal-to-noise ratio, and physiological artifacts can easily result in false dynamic connections, which are not originated form neural interactions. In addition, how dFC is estimated is influential in the detection of statistically significant dynamic connections (Hindriks et al., 2016;Lindquist et al., 2014;Savva, Mitsis, & Matsopoulos, 2019). For example , Hindriks et al., (2016) have claimed that it is impossible to detect dynamic connections using the SWC method in individual sessions through simulation studies and validated this claim by using both the rs-fMRI data of the human and the macaque; however, they have reported that averaging the statistical measure across subjects/sessions could increase the power of detecting dynamic connections. In another study, Savva et al., (2019) have shown that mutual information and variation of information yield most consistent results by achieving high reliability with respect to dFC estimations for different window sizes in comparison with correlation metrics such as Pearson correlation, Spearman and Kendall correlation, and Pearson and Spearman partial correlation.
Thus, their findings suggested that how dFC is estimated, greatly affects the power of detecting dynamic connections. In consequence, it has recently become critical to determine whether the estimated dFC is in fact due to neuronal interactions or random noise (Hindriks et al., 2016;Kudela, Harezlak, & Lindquist, 2017;Leonardi & Van De Ville, 2015;Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014). Therefore, detection of statistically significant dynamic connections is essential for dFC studies.
An appropriate statistical framework is required to determine whether the observed variation in the FC time series can be characterized as dynamic pattern or whether it is due to statistical uncertainty (Hindriks et al., 2016;Sakoglu et al., 2010). To this end, a commonly used approach is to calculate a test measure that characterizes the fluctuation in the FC time series and subsequently applying a statistical hypothesis test. In this framework, the null hypothesis states that the estimated dFC time series is static and is evaluated on the basis of the distribution of the calculated test measure. In the literature, several test measures have been proposed to test the presence of dynamicity in the estimated dFC time series, including the variance of the dFC time series (Hindriks et al., 2016;Sakoglu et al., 2010), a linear measure based on the dFC time series' Fourier transform (Handwerker et al., 2012), and a nonlinear measure (Zalesky et al., 2014). Since the null distributions of the measures cannot be analytically derived, surrogate data, produced based on the statistical properties of the observed data, are used and dynamic connectivity tested based on a test statistics measure that reflects the dynamicity of the estimated dFC time series (Pereda, Quiroga, & Bhattacharya, 2005). Considering the variance of the dFC time series (σ 2 ) as the test measure, in the absence of dynamicity (null hypothesis), this measure is expected to be only due to statistical uncertainties and remain relatively small over time. On the other hand, in the presence of dynamicity (alternative hypothesis), this measure will not be only due to statistical uncertainties and becomes relatively large. In other words, the variance under the null hypothesis is positive but statistically smaller than that under the alternative hypothesis. Consequently, if the variance of the estimated dFC is located in the upper five percentile of the null distribution, the null hypothesis will be rejected with p-value <.05. This is an evidence for the presence of dynamicity in the estimated dFC time series (Chang & Glover, 2010;Hindriks et al., 2016;Zalesky et al., 2014).
Since one of the main factors that affects the statistical power in the detection of dynamic connections is how dFC is estimated, developing a powerful method to estimate dFC with high accuracy is of critical importance. The SWC-based methods, as the conventional methods to estimate dFC, suffer from some limitations that can impact the interpretation of the findings (Hindriks et al., 2016;Savva et al., 2019) as follows: 1. Low detection power: SWC method uses equal weights across all observations within a window (Lindquist et al., 2014), which in turn leads to variations in the estimation results (Hindriks et al., 2016;Kudela et al., 2017;Lindquist et al., 2014). In consequence, spurious fluctuations caused by noise can easily show up as dynamic changes in the estimated dFC. Hence, the quality of the estimated dFC has an important effect on the power of detecting dynamic connections. Furthermore, the bivariate nature of SWC, which only captures the strength of association between pairs of brain regions, might be another explanation for this limitation. This is because multiple brain regions are engaged in cognitive tasks and resting state conditions (Anzellotti, Caramazza, & Saxe, 2017;Gallagher & Frith, 2003) and using bivariate measures to estimate interactions between these regions may not explore the neural bases of behavior or cognition. Recent exploration of uncertainty in estimation of dFC has reported issues due to stationarity and statistical testing of dFC (Liegeois, Laumann, Snyder, Zhou, & Yeo, 2017). Parametric approaches show greater power in detecting dFC changes. For example, Liegeois et al., (2017) have suggested that Autoregressive models are powerful tools for exploring the dynamical properties of rs-fMRI. They also explored different frameworks including phase randomization and autoregressive randomization for generating surrogate data for statistical testing of dFC. Their findings showed that bivariate autoregressive randomization approach is prone to false-positives compared with phase randomization and multivariate autoregressive randomization approaches.
2. Appropriate window length: Setting the length of time window is very critical and affects the connectivity results (Shakil, Lee, & Keilholz, 2016). Using a long window risks to miss fast changes in FC evolution over time, whereas a short window will reduce effective sample size and make the estimation procedure unreliable (Hutchison, Womelsdorf, Allen, et al., 2013). However, some efforts have been made to address the algorithmic selection of the window length to explore dFC (Vergara, Abrol, & Calhoun, 2019). Vergara et al., (2019) proposed to use an averaged SWC, which requires a window length smaller than that of SWC. This is important because shorter windows allow for more accurate estimation of transient dynamicity of FC. Including an averaging step in the processing of SWC as proposed in (Vergara et al., 2019) provides a method for eliminating artifact fluctuations due to windowing compared with the common SWC. In this way, the averaged-SWC identifies dFC fluctuations better than the common SWC.
3. Sparsity of dFC networks: The dFC networks resulting from the SWC method are fully dense because of the presence of noise and other nonneuronal sources that contribute to the acquired BOLD signals, whereas numerous studies have reported that brain networks are of parsimonious structure which enable brain to process and transfer information with high efficiency and low redundancy (Bullmore & Sporns, 2009). Thus, estimation of the dFC network structure is a problem (A. Liu et al., 2015;Monti et al., 2014), which is inherently relevant in the estimation of the network's sparsity.
In this work, we adopt the kernel-reweighted logistic regression (KELLER) from genetic studies (Song, Kolar, & Xing, 2009), which provides us with a comprehensive framework to estimate dFC networks while overcoming limitations of SWC. We change the definitions of the variables in order to match them with those of the rs-fMRI data.
To address the first mentioned limitation, we model dependencies of brain regions by considering multivariate relations between BOLD signals of pairs of brain regions in KELLER. Moreover, kernel-reweighted feature of KELLER weighs observations unevenly in each window in a way that the adjacent observations have stronger contributions to the estimation process than the distant observations. Thus, we hypothesize that this modeling leads to a more accurate estimation of the brain dynamic interactions. To overcome the second limitation, KEL-LER utilizes a sophisticated parameter selection technique based on the Akaike information criterion (AIC). Finally, KELLER ensures the sparsity of the estimated dFC networks by applying an ℓ1-penalized term in the loss function which effectively yields a sparse network.
To utilize KELLER in estimating dFC from rs-fMRI data, we define a new time-varying network model based on temporal modeling of rs-fMRI time series in which we model the multivariate probability density function (pdf) of the rs-fMRI time series of all brain regions at each time point by using a pairwise Markov Random Field (MRF) model. In this model, dFC between each pair of brain regions indicates the strength of undirected interaction between them. The MRF model has been appealing in the analysis of complex dependence structures (Kaiser, 2007). We reformulate the multivariate pdf of the rs-fMRI time series of brain regions into a product of conditional pdfs. As a result, the problem of estimating dFC networks is decomposed into one of estimating a series of distinct and static FC networks. Moreover, this step provides an opportunity to estimate multivariate relations between brain regions by estimating functional interaction pattern of a brain region and other regions simultaneously at each time point, using a neighborhood selection procedure (Song et al., 2009;Wainwright, 2006). The resulting functional pattern of each brain region over time is referred to as dynamic neighborhood vector in the rest of the paper. In this way, we can estimate dFC networks by putting together all estimated dynamic neighborhood vectors related to all brain region with a temporal resolution equal to the sampling rate of the BOLD signal. Subsequently, null hypothesis significance testing based on the amplitude-adjusted phase randomization procedure surrogate data generation (Betzel, Fukushima, He, Zuo, & Sporns, 2016) is employed for detecting dynamic connections (Chen, Rubinov, & Chang, 2017). One possible approach towards obtaining such a statistical assessment is to use a statistic measure that characterizes the changes in the estimated dFC time series (Hindriks et al., 2016;Savva et al., 2019). Then, the estimated statistic measure is tested through its null distribution approximated using surrogate data to verify the presence of dynamic pattern in the estimated dFC time series. Thus, the second hypothesis in this paper is that utilizing KEL-LER for estimating dFC networks increases dynamic pattern detectability in estimated dFC time series due to modeling the multivariate relations between BOLD signals of brain regions. Moreover, the ability of KELLER to automatically estimate true sparse structure of dFC network at each time point increases the accuracy of the estimated dFC networks.
As mentioned above, SWC has low power of dynamic connection detection because of weighting all observations equally in each window (Lindquist et al., 2014). This limitation of SWC can be overcome by using tapered windows. Thus, we evaluate the performance of KELLER on a series of simulation studies and real rs-fMRI data in comparison with the SWC-based methods including SWC and Tapered-SWC (T-SWC) (see Section 2.5). Moreover, since the SWC-based approaches can only obtain time-varying estimates of the covariance matrices, for fair comparison, we apply the graphical lasso (Friedman, Hastie, & Tibshirani, 2008) subsequently to learn true sparsity structure in the dFC networks. The combined methods are referred to as SWCGL (SWC and graphical lasso) and T-SWCGL (T-SWC and graphical lasso) in the rest of the paper.
Another important issue in the SWC method is the effect of window size on the estimated dFC time series. In the literature, a suitable window for dFC is suggested to be between 30 and 100 seconds (s) (Wilson et al., 2015). On the other hand, a rule of thumb has been proposed by Leonardi and Van De Ville for removing spurious fluctuations due to inappropriate window length. They set window length to 1/f min s or larger in SWC, where f min corresponds to the smallest frequency in the spectrum (Leonardi & Van De Ville, 2015;Zalesky & Breakspear, 2015). The spectrum of fMRI BOLD signals has been proposed to start at 0.01 Hz after studying frequencies dominated by neuronal activity and away from physiological noise such as cardiac and respiratory activities (Chen & Glover, 2015). Moreover, Zalesky and Breakspear (Zalesky & Breakspear, 2015) have provided further statistical support for this rule of thumb, but suggested that if fMRI data has a moderate SNR, the window length of 1/f min s may be overly conservative and in this case, dFC can in theory be detected with much shorter window lengths (e.g., 40 s). They have also suggested that statistical testing and appropriate surrogate data is crucial in this respect. Thus, in this work, to minimize the effect of window length on the capability of SWC-based methods, we use different window lengths from 20 to 140 s in 20 s steps.
The remainder of the paper is organized as follows. In the next section, we will first introduce the KELLER algorithm in detail for computing dFC networks from rs-fMRI data. Next, in the Materials and Methods Section, we explain simulation studies as well as real rs-fMRI data. Then, we describe how to estimate dFC time series by utilizing KELLER and detecting dynamic connections. In the Results Section, we present the results of the simulation studies to evaluate the performance of KELLER in comparison with SWC-based methods. We also present the findings of applying KELLER on the rs-fMRI data to study the dynamic behavior of the whole brain in healthy subjects.

| KELLER ALGORITHM
In this section, we introduce the KELLER algorithm (Song et al., 2009) for computing dFC networks from rs-fMRI data. First, a dFC network model based on temporal modeling of rs-fMRI time series is described.
Afterwards, we explain the core of KELLER for estimating dFC networks as the estimation of a dynamic neighborhood vector. Finally, we discuss how parameters of KELLER are set. In the following, matrices, vectors, and scalars are denoted by boldface capital letters, boldface lowercase letters, and lowercase letters, respectively.

| Dynamic FC network model based on temporal modeling of rs-fMRI time series
Let us consider the representative rs-fMRI time series of p regions of interest (ROIs) at a given time point t as a random vector denoted by We suppose that T observations of y are available denoted by YR T×p . Representative rs-fMRI time series of each ROI is extracted as the mean rs-fMRI time series of all voxels in that ROI.
Prior biological knowledge of rs-fMRI data allows us to hypothesize that there may be a meaningful correlation at a given time point t between each pair of y t ð Þ m and y t ð Þ n , m, n1 : p. Our objective is to estimate dFC matrices {θ (t) , t = 1 : T} ≔ {θ (1) , …, θ (T) } between all pairs of y (t) R p over time, that is, an FC matrix for every time point of the rs-fMRI measurement. To consider multivariate interactions between brain regions in estimating dFC matrices, we used the KELLER algorithm originally proposed in the genetic studies framework as a generative model based on a pair-wise MRF which represents the multivariate pdf of the random vector y (t) at a given time point t. To translate KELLER from genetic to rs-fMRI, we need to define a dichot- 1 , …, y t ð Þ p R p at a given time point t will be within the interval (0,1) and then can be dichotomized to "High" or "Low" level of func- gby setting a fixed threshold of 0.5. Note that this thresholding is a kind of adaptive thresholding based on the variation of functional activity in a given ROI, because the representative rs-fMRI time series of the p regions of interest (ROIs) were separately normalized between zero and one before thresholding. Thus, the generative model can be defined based on pair-wise MRF which represents the joint probability of measured BOLD signal of all ROIs at time point t, y (t) as follows: where θ t ð Þ mn encodes the undirected interaction between the rs-fMRI signals of each pair of ROIs (m and n) at time point t θ t ð Þ . In other words, dFC between each pair of brain regions is modeled as mn which indicates the strength of undirected interaction between them. The partition function Z(θ (t) ) in the MRF model normalizes Equation (1) to a probability function. In MRF modeling, the dependencies of a set of random variables are represented by an undirected graph. Therefore, dFC matrices {θ (t) } are expected to be undirected, weighted, and symmetric. The presence of a non-zero entry in the dFC matrix θ (t) means that the fluctuations of the rs-fMRI time series of the corresponding ROIs are functionally related at a given time point t. Changes in the value of parameter θ t ð Þ mn over time is considered as the dynamic interaction between ROIs (m and n). In other words, if the value of θ t ð Þ mn increases (or decreases), it reflects that the interaction between ROIs (m and n) becomes stronger (or weaker) or connection between them appeared (or disappeared). Accurate estimation of a series of dFC matrices {θ (t) , t = 1 : T} ≔ {θ (1) , …, θ (T) } is the main focus of this study.
From the perspective of the brain functional organization, we impose the following two properties on the estimation procedure of the dFC matrices {θ (t) }. Since the topology of FC networks changes smoothly over time (Lin et al., 2017), the first property is the smoothness of variation in dFC pattern over time. In mathematical terms, smooth changes of dFC networks means that the changes in is not practically feasible by maximizing the log-likelihood of the joint probability function in Equation (1),  (1) can be simplified as a set of conditional probabilities of the level of functional activity at each ROI based on the functional activity of the rest of ROIs at time point t. These neighborhood vectors reflect multivariate relations between a brain region and the rest of the brain regions at a given time point t. Afterwards, we join the corresponding dynamic neighborhood vectors to recover the overall dFC network at each time point. It is worth mentioning that dynamic neighborhood vector of each ROI at a given time point t defines the multivariate functional pattern of a particular brain region with other brain regions.

| Estimation of neighborhood vector
In this new framework, we employ neighborhood selection procedure (Song et al., 2009;Wainwright, 2006) to convert estimation of dFC matrices at each time point to estimation of a sequence of neighbor- In other words, estimating dFC network is equivalent to recovering the structure of interactions of each ROI (n) with the rest of ROIs. In fact, if we can correctly estimate the neighborhood vectors, it will lead to exact recovery of dFC networks. Therefore, the joint probability function in Equation (1) is decomposed into the product of conditional probability functions of y , which represents the conditional probability of the functional activity of ROI n at a time point t, given the measured BOLD signal of all ROIs except ROI n at time point t, Here, we justify how the joint probability function in Equation (1) is decomposed into the product of conditional probability functions.
As mentioned in Section (2.1), if two distinct ROIs m and n are con- (1) can be simplified as follows: where T v 2 is the inner product. We use n 0 to determine all ROIs excluding {n}, that is, n 0 ≔ {1, …, p} − {n}; n = 1 : p. Now, we can rewrite this joint probability function by using only conditional probabilities. Based on the chain rule, it is proven that and the above equation based on the conditional probabilities would be: Because of the upper-triangular property of the θ (t) matrix, each probability of the functional activity of each ROI at a time point t, given the measured BOLD signals of the rest of the ROIs at time point t, y t ð Þ n 0 . As described in Section 2.1, in the new framework adopted from the KELLER algorithm, we used a binary block to define the level of functional activity of each ROI at every time point and assumed a linear relationship between the functional activity level of ROI n and those of all other ROIs except ROI n at time point t and defined d t ð Þ n as a binary variable of y n is a binary variable, it would be mathematically possible to assume that Banerjee, 2007). In our case, the log-odds ratio follows the following equation: where ℓ is the log-odd, and hv 1 , v in our case, all conditional probabilities follow the logistic function as: In this model, θ It is not possible to estimate θ t ð Þ n 0 n o by directly minimizing Equation (7)  n o using Equation (7) and ensure that the estimated θ t ð Þ n 0 n o varies smoothly over time, we introduce the following kernel reweighted function w (t) (t*): where k h Á ð Þ = k Á h À Á is a nonnegative symmetric kernel and h is a bandwidth parameter that controls the kernel size. In this work, we define k h (Á) as a radial basis function Gaussian kernel as k h (t) = exp(−t 2 /h) in a way that the adjacent observations have stronger contributions to the estimation than the distant observations. It is noteworthy that weighting the observations is used in other methods such as shorttime Fourier transform to extract the transient frequency components (Ahmed & Xing, 2009;Song et al., 2009).
Finally, sparsity is introduced into the model by using an ℓ1-norm regularization term which assigns a large penalty to vectors with large absolute values. In this way, the penalty term shrinks elements to zero effectively. KELLER minimizes the following loss function: We use k.k 1 for the ℓ1-norm, v k k 1 = P p n = 1 v n j j. In Equation (10), δ ≥ 0 is a constraint to control the magnitude of the estimated dynamic neighborhood vectors and the sparsity of the dFC network defined by combining these neighborhood vectorŝ
In the estimation of dFC networks during resting state or a cognitive process, there are tens of subjects, hundreds of ROIs, and hundreds of time points, and hence about a million optimization problems. Therefore, it is crucial to choose an efficient algorithm for solving the minimizing problems defined in Equation (10) to minimize the overall computation cost. Here, we parallelized the optimization procedure across different ROIs and different time points by implementing the projected gradient (PG) method (Duchi et al., 2008) because of its simplicity and efficiency. Since ℓ1-regularized logistic regression loss function can be reformulated as a constrained minimization problem, we can rewrite Equation (10) as follows: where C δ is the upper bound of the first order norm of θ t ð Þ n 0 and determines the area (Ω) that contains all the estimated parameters. A oneto-one correspondence exists between the penalty parameter δ in Equation (10) and C δ in Equation (11). In the new formulation, the objective function L θ t ð Þ n 0 is a convex function and its derivative with respect to vector θ t ð Þ n 0 is obtained as follows: In the PG method, the parameters are updated in line with a negative gradient. Following an update, if the parameter is outside the Ω area, it is projected back into the Ω area. Otherwise, the algorithm goes to the next step. The basic step in this method, which guarantees its performance, is the projection of the parameter into the Ω area: where Π Ω (a) = argmin b {ka − bk | b Ω} is the Euclidean projection of vector a into the Ω area (Duchi et al., 2008). The implemented version of PG algorithm for the optimization problem in Equation (11) is described in the Supporting Information (S1). It should be noted that the PG method has several internal parameters, such as α, ε, and σ, which are adjusted in accordance with (Bertsekas, 2016).

| Parameter selection
The proposed KELLER method requires two input parameters h and δ, which can be adjusted using the available data. The parameter h is the width of the Gaussian kernel. This parameter is the most important factor in controlling the temporal smoothness of the estimated dFC networks. A large kernel size allows for more observations to be used in the estimation of the dFC networks while increases the possibility of losing rapid changes in the dFC network. On the other hand, a small kernel size increases the sensitivity to rapid changes, while, the estimation variance increases due to a drastic drop in the number of observations used for the estimation.
The parameter δ controls the sparsity. In particular, a large δ will result in a network that has a high degree of sparsity. Therefore, determination of both h and δ parameters is very important. We h and δ, we define AIC as: where N is the estimated number of degrees of freedom and equals It is worth mentioning that the Bayesian information criterion (BIC) has been used to tune hyper-parameters (Song et al., 2009).
However, BIC selects the correct model if an infinite amount of data are available (Burnham & Anderson, 2002) or the correct model is among a set of candidates (Olofsen & Dahan, 2015). Since in our application, there is no guarantee that the correct model belongs to a set of candidates, we use AIC.

| Tapered sliding window correlation (T-SWC)
T-SWC is identical to SWC but it uses weighted Pearson cross-correlation. As mentioned previously, SWC uses equal weights for all observations in a window (Lindquist et al., 2014), which in turn leads to variations in the estimation results (Hindriks et al., 2016;Kudela et al., 2017;Lindquist et al., 2014 fair comparison, we also apply the graphical lasso (Friedman et al., 2008) subsequently on the results of T-SWC to learn the sparsity structure in the estimated dFC networks (see Section 2.5.3). This combined method is referred to as T-SWCGL in the rest of the paper.

| How the graphical lasso was applied on the results of SWC and T-SWC
To apply graphical lasso on the results of SWC and T-SWC, we first convert the calculated correlation matrix at the t th window to the covariance matrix by: Here, the ij th element of the covariance matrix is related to the corresponding element of the correlation matrix by the above formula where σ i and σ j are the standard deviation (SD) of the i th and j th variables at the t th window. Then, the corresponding precision (inverse covariance) matrix at the t th window is estimated while considering sparsity in its structure by using sparse inverse covariance estimation with the graphical lasso proposed by Friedman et al., (2008). The point which should be noted is that the sparsity in KELLER is inherited in the algorithm however, for the SWC-based methods, we do it as a post-processing step.

| Simulated data generation and analysis
In this section, we evaluate the performance of KELLER in estimating the dFC networks in comparison with T-SWCGL and SWCGL. To generate simulated data, we consider well-known properties of functional brain organization such as high positive temporal autocorrelation of BOLD signals (B. Biswal et al., 1995;Friston, 2011) and self-organization and scale-free characteristics of brain networks (Eguiluz, Chialvo, Cecchi, Baliki, & Apkarian, 2005;Lee et al., 2010;X. Liu, Ward, Binder, Li, & Hudetz, 2014). Since, in the literature (Liegeois et al., 2019;Monti et al., 2014;Rogers, Katwal, Morgan, Asplund, & Gore, 2010;Valdes-Sosa et al., 2005), the first-order Vector Autoregressive (VAR) processes are used to evaluate dFC networks, we generated simulated rs-fMRI data based on the first order VAR process. The VAR process is well suited to the task of producing auto-correlated multivariate time series as they are capable of encoding autocorrelations within components as well as crosscorrelations across components (Cribben et al., 2012;Monti et al., 2014). In order to evaluate the performance of different methods in estimating dFC networks, simulated rs-fMRI datasets were generated based on a first order VAR model with pre-defined temporal autocorrelation structures and modulation (Deler & Nelson, 2001). We studied the performance of the proposed algorithm by using two types of random graphs as the structure of a pre-defined autocorrelation network: Erd} os-Rényi random graphs (Erdos & Renyi, 1959) and scalefree random graphs obtained by using the preferential attachment model of Barabási and Albert model (Barabási & Albert, 1999). Erd} os-Rényi random graphs are the simplest and most widely studied type of random networks while the use of scale-free networks is motivated by the fact that they resemble some aspects of fMRI networks. For example, previous studies (Eguiluz et al., 2005;van den Heuvel, Stam, Boersma, & Hulshoff Pol, 2008) suggested that the degree distribution of the resting state functional brain organization follows a power law.
Moreover, it has been shown that the self-organization property of the functional brain organization is linked with the power-law (scalefree) scaling property of functional brain organization (Gisiger, 2001).
In the case of scale-free networks, the power of preferential attachment (new connection) was set to unity on the rs-fMRI networks, similar to previous studies (Monti et al., 2014). Additionally, we used pre-defined temporal modulation of autocorrelation matrices in VAR process to ensure that the dynamicity is inherited in the simulated data with gradual changes within a state and that changes from one state to another state are smooth and thus without any abrupt changes.
Schematic overview of how we generated simulated rs-fMRI data are illustrated in Figure 1a. We generated simulated dataset based on the first order VAR process while the temporal structure of network and the temporal modulation of connections were predefined. For the simulated data, we considered three states and the correlation structure of each state was randomly generated by Erd} os-Rényi random networks (Erdos & Renyi, 1959) or scale-free random networks (Barabási & Albert, 1999)

| Performance measure to evaluate methods
The goal of this study is to estimate the dFC networks at a sampling rate that leads to a correct estimation of the nonzero elements in the estimated dFC matrices,θ t ð Þ ,t = 1,…, T , so we evaluate the performance of the estimation procedures using an F1 score (Chinchor, 1992) (Figure 1b). All of the nonzero entries inθ percentage of true edges estimated by the algorithm. The precision and recall are calculated as follows: The F1 score attempts to balance between the Pre and Rec as a prevalent metric of the performance measure. Finally, the F1 (t) score is defined as High performance in retrieving the true structure of the network depends on having a high F1 (t) score, which in turn requires that both the precision and recall are also high.

| Defining dynamic connections in the benchmark model in the simulation study
Although we simulated dynamicity in the temporal modulation of functional connections as a predefined pattern by gradual emergence or disappearance of connections, we applied a statistical hypothesis analysis on the simulated dFC time series to define connections whose simulated fluctuations in their temporal modulation pattern were due to their dynamic nature. In this analysis, the distribution of the calculated test measure is constructed based on the null hypothesis which corresponds to temporal modulation pattern of the simulated connection being static while the alternative hypothesis corresponds to being dynamic. As illustrated in Figure 1c, we first calculated the variance of the predefined temporal modulation of connections (or equivalently, their SD) as the statistical measure (the observed measure). Subsequently, we constructed 500 surrogate sets for the simulated rs-fMRI data set based on the amplitude-adjusted phase randomization procedure (Betzel et al., 2016

| Real data
We used open access data from the imaging center of the

| Preprocessing of real data
The rs-fMRI and anatomical data were preprocessed using Statistical

| Statistical assessment of estimated dFC time series
To detect true dynamic connections, statistical assessment of the estimated dFC time series is essential for all dFC studies. Therefore, a proper statistical framework should be applied to determine whether the observed variation in the estimated dFC time series can be characterized as dynamic pattern or it is due to statistical uncertainty (Hindriks et al., 2016;Sakoglu et al., 2010). To this end, a commonly used approach is to calculate a test measure that characterizes the fluctuation in the estimated dFC time series by applying a statistical hypothesis test with a null hypothesis which is constructed based on the distribution of the calculated test measure. The null hypothesis states that the estimated dFC time series is static.

| Hypothesis testing and statistic measure
In this study, we focus on the variance of the dFC time series Thus, D contained the estimated dFC time series of all ROI pairs. Now, the test measure is represented by the SD of the estimated dFC time series for the n th ROI pair as: whereθ n =θ n with a p-value <.05 and conclude that the estimated dFC time series for the n th ROI pair is dynamic. To estimate the distribution of σ under the null hypothesis, we use randomized data, known as surrogate data, similar to the analysis of nonstationary time series and dFC (Betzel et al., 2016;Chang & Glover, 2010;Prichard & Theiler, 1994;Zalesky et al., 2014).

| Surrogate data generation
To approximate the null distribution, we construct 1,000 surrogate sets of BOLD time series for the 112 ROIs of the Harvard-Oxford atlas for each rs-fMRI scan in the data set, using an approach similar to that presented in (Betzel et al., 2016). This approach is based on the amplitude-adjusted phase randomization procedure in which surrogate BOLD time series are generated with randomized phase, but with the same amplitude distribution so as to preserve the static FC pattern of the real data.

| Estimating dFC time series and detecting dynamic connections
To assess the power of KELLER in detecting dynamic connections and compare with those of SWC based method, we applied KELLER and T-SWCGL methods on real rs-fMRI data set. A graphical summery of processing steps on the real rs-fMRI data based on KELLER method is illustrated in Figure 2. After preprocessing, we estimated individual dFC matrices from 112 extracted time series based on Harvard-Oxford atlas (Bohland, Bokil, Allen, & Mitra, 2009)   with p-value <.001 (Table 1).
In Figure 3d and   Note: Simulated network with the following setup: network topology = Erd} os-Rényi random graph; number of nodes = 5; Detailed information is given in Figure 3. Abbreviations: KELLER, kernel-reweighted logistic regression; SD, standard deviation; SWCGL, sliding window graphical lasso; T-SWCGL, tapered sliding window graphical lasso.
Additionally, we evaluated whether an increase/decrease in window length influenced the results of T-SWCGL and SWCGL. We ran this simulation study 100 times with 5-nodes and random network structure. We also applied paired-samples t test to assess if the "average" values obtained by each window length was "significantly" better than those of other window lengths. The results showed that there was no significant difference between the obtained results for different window lengths. However, both T-SWCGL and SWCGL led to comparable performance with regards to F1-score and dynamic connection detection power for a window length of 100 sec (Tables 2   and 3).  Tables S2 and S3. It is evident that for the number of nodes larger than 40, this ratio decreases considerably compared with the number of nodes smaller than 40. In other words, if the number of nodes goes to infinity, the ratio of the false positive to false negative values tends to unity. Thus, the performance of methods for the number of nodes larger than 40 can be considered as the steady state performance. As indicated in the Section 2.5.3 the sparsity estimation for the SWCbased methods was done as a post-processing step by applying the graphical lasso, and we also compared the estimated sparsity of SWC and T-SWC with those of KELLER. We reported the mean ± SD of the resulting sparsity in the simulated networks by different methods in the simulation study with over 100 runs (Tables S4 and S5). Subsequently, we applied paired-samples t test to assess whether there were any significant (p-value <.05) differences between the results.
However, no significant differences were found, meaning that the performance of KELLER, with inherited property of forcing the estimated dFC to be sparse, is similar to the performance of SWCGL and the combination of SWC-based methods with the graphical lasso.
Moreover, for each of the different methods, we assessed the percentage of the observed fluctuations in the estimated dFC time series that were due to its dynamic nature and not statistical uncertainty in comparison with the benchmark model. In Figure 5, we illustrate the mean percentage of statistically significant dynamic T A B L E 2 Performance of SWCGL method in estimating dynamic connections in terms of F1 score, dynamic detection power as the window length was increased from 20 to 140 s with 5-nodes Erd} os-Rényi random network topology structures.    Table S6. Results in Figure 5 show that their patterns are similar to those obtained in the performance of methods in estimating true structure of dFC networks illustrated in Figure 4. In both simulation studies of random and scalefree networks, KELLER worked more efficiently in detecting dynamic connections since this method estimated dFC time series more accurately than the other methods as demonstrated in Figure 4.

| Experimental results
We applied the KELLER and T-SWCGL methods to real rs-fMRI data and compared the results. We did not apply the SWCGL method to real rs-fMRI data since this method was outperformed in estimating dFC networks over time by the other two methods in the simulation studies.

| Detecting dynamic connections in healthy subjects
The first row in Figure 6, depicts the results achieved by KELLER while the following rows show the results of T-SWCGL with different window lengths from 20 to 140 s, spaced 20 s apart. In Figure 6 (column a), we demonstrated the SD of the estimated dFC time series for all ROI pairs of each subject, as well as the averaged case. The related p-values of all those ROI pairs were obtained from statistical hypothesis testing and adjusted for multiple comparisons using Bonferroni correction. These results surpassed the 5% significance threshold for each individual subject as well as in the averaged case are illustrated in Figure 6 (column b). The histogram of the SD value of the estimated dFC time series for all ROI pairs for all subjects is depicted in Figure 6 (column c). Figure 6 (column d) shows the number of connections with statistically significant dynamic pattern. We found that in most plots shown in Figure 6 (column d), for an individual, a few connections have significant dynamic pattern while for the averaged case, this number increases remarkably. As previously mentioned, in the averaged case, we averaged the observed and the surrogate test measure values for each ROI pair connection across subjects. We then applied null hypothesis testing to define which connections were dynamic. In fact, the averaging process naturally increased the statistical power of the null hypothesis testing in detecting dynamic connections. However, findings of KELLER and T-SWCGL with window length of 100 s showed different patterns. The number of significant dynamic connection for individuals in both results obtained by KELLER and T-SWCGL with window length of 100 s are more than those of T-SWCGL with window lengths of 20, 40, 60, 80, 120, and 140 s. Figure   6 (column E) demonstrates the mean dFC network for the averaged T A B L E 3 Performance of T-SWCGL method in estimating dynamic connections in terms of F1 score, dynamic detection power as the window length was increased from 20 to 140 s with 5-nodes Erd} os-Rényi random network topology structures.  Patterns are similar to those in Figure 4. In both network categories, KELLER is more efficient than the other methods in detecting dynamic fluctuations. Detailed information is provided in the Supporting Information (Table S4). dFC, dynamic functional connectivity; KELLER: kernelreweighted logistic regression; T-SWCGL: tapered sliding window correlation + graphical lasso; SWCGL: sliding window correlation + graphical lasso case across subjects with statistically significant dynamic connections adjusted for multiple comparisons using Bonferroni correction.
To compare the results obtained by KELLER with those of T-SWCGL considering different window lengths, the similarity between the mean estimated dFC matrix by both methods with statistically significant dynamic connections was calculated for each subject ( Table   4). The results revealed that in 68% of the subjects, T-SWCGL with window length of 100 s had the highest similarity to KELLER. The range of similarity varied from r = 0.37 to r = 0.77 for different window lengths. The findings are in line with the rule of thumb proposed by Leonardi and Van De Ville (2015) which suggests that appropriate SWC window length should be set to 1/f min s or larger, where f min corresponds to the smallest frequency in the spectrum (Leonardi & Van De Ville, 2015;Zalesky & Breakspear, 2015). In this study, f min equals 0.01 Hz and the window length for SWC should be equal to or larger than 100 s. In the averaged case, the highest similarity between KELLER and T-SWCGL was achieved for a window length of 60 s (r = 0.78), followed by window lengths of 40 and 100 s (r = 0.71).  Table S7. These ROIs are selected in association with seven resting state networks of functionally coupled parcellated regions across the cerebral cortex (Yeo et al., 2011) as well as subcortical regions including amygdala, hippocampus, and para-hippocampal regions. In Figure 7, we illustrated the dynamic connections in DMN yielded by different methods. We also listed the common dynamically connected ROI pairs identified by KELLER and T-SWCGL in Table 5. Results showed that most of the commonly identified dynamic connections in DMN were between subcortical, temporal, and frontal regions. On the other hand, the list of dynamic connections in DMN identified only by specific methods are presented in Table 6. F I G U R E 6 ROI-pairs analysis of dFC time series by KELLER (first row) and T-SWCGL with different window lengths including 20 up to 140 s in the following rows. SD value of the estimated dFC time series for all ROI pairs for each subject, as well as the averaged case (column a); Corresponding calculated p-values of all those ROI pairs by statistical hypothesis test and adjusted for multiple comparisons using Bonferroni correction, crossed the 5% significance threshold, for each individual subject as well as the averaged case (column b). Histogram of the SD value of the estimated dFC time series for all ROI pairs for all subjects (column c). Number of connections with statistically significant dynamic pattern (column d). We found that in the most plots in column d, for an individual, very few connections have significant dynamic pattern while for the averaged case this number increases substantially. However, findings of KELLER and T-SWCGL with window length of 100 s showed different pattern. In these cases, the number of significant dynamic connection for individuals are more than those of T-SWCGL with window length of 20, 40, 60, 80, 120, and 140 s. Mean dFC network for the averaged case across subjects with statistically significant dynamic connections adjusted for multiple comparisons using Bonferroni correction (column e). dFC, dynamic functional connectivity; KELLER, kernel-reweighted logistic regression; T-SWCGL, tapered sliding window correlation + graphical lasso 5 | DISCUSSION

| Overview of the current study
The main goal of this study was to improve the power of detecting dynamic connections in estimated dFC by conventional methods.
Thus, we introduced a framework developed in the gene regulatory studies, the KELLER algorithm (Song et al., 2009), and utilized it in the dFC network studies. KELLER allowed for the retrieval of the dFC pattern of brain networks in terms of the network's structure and modulation over time.
A series of simulation studies were done to measure the capability of KELLER in retrieving the underlying structure of dFC network as well as the power of KELLER in detecting dynamic connections in comparison with T-SWCGL and SWCGL. We also employed the proposed method in estimating whole brain dFC networks from rs-fMRI data of 31 healthy subjects to detect statistically significant dynamically connected brain region pairs. To achieve this, we performed proper statistical tests on the SD of dFC time series as a test statistic measure using appropriate surrogate data.
We demonstrated the following key results via simulation studies:

| Previous studies in assessing dFC in rs-fMRI by using SWC technique and comparison with the present study
The effect of different sliding window parameters such as window type, length, and step, as well as different FC metrics on the detection of dynamic connections or brain states have been investigated (Hindriks et al., 2016;Savva et al., 2019;Shakil et al., 2016). In (Shakil et al., 2016), segmented real BOLD time series were mixed to form a simulated setting where the switching between brain states were known. Their findings implied that window length and size affected the identification of brain state switching (Shakil et al., 2016). Similar to our findings on experimental fMRI data showed that window length influenced the results of T-SWCGL method. In another study, they tried to answer similar questions by applying different FC metrics to estimated dFC network and using experimental rs-fMRI data in two separate groups for test-retest validation (Savva et al., 2019). The authors showed that the obtained results based on mutual information and variation of information with a window length larger than 120 s yielded the most consistent results by applying test-retest analysis. Moreover, in (Hindriks et al., 2016) it was claimed that it was impossible to detect dynamic connections by using SWC in individual sessions through simulation studies and this claim was validated using both rs-fMRI data of humans and macaques (Hindriks et al., 2016). Recently, the possible impact of global signal regression on the estimation of dynamic connection and brain states was evaluated by H. Xu et al., (2018).  Figure S3). It is noteworthy that most of the previous dFC studies including (Hindriks et al., 2016;Shakil et al., 2016) have not considered the impact of global signal regression on their findings.
T A B L E 5 Common dynamic connections in DMN identified by KELLER and T-SWCGL with different window length (w = 40, 60, and 100 s).

dFC connections in DMN
Common dFC connection in DMN between KELLER and T-SWCGL (w = 40, 60, and 100 s) In the present study, we used KELLER (Song et al., 2009) to estimate dFC networks at each time point of BOLD signal measurement.
The important feature of KELLER is that it considers multivariate dependencies between brain regions to estimate dFC networks through estimating functional pattern of each brain region spatially and temporally. SWC-based methods only capture bivariate linear association between brain regions. Additionally, KELLER also has the potential to define a proper window length to extract dFC time series, utilizing AIC as a model selection approach. In this study, for KELLER, we defined a range of values (from 10 to 150 with 10 s steps) for the window length parameter, h, Then, we employed an extensive gridsearch in the given range of the parameters h and δ. We used AIC to estimate the in-sample prediction error for each choice of parameters h and δ, allowing for a direct comparison across different values of each parameter. Finally, a pair of parameters that minimizes AIC is chosen as the optimal values. Moreover, the structure of dFC networks is automatically estimated at each time point because of the ℓ1-penalized term in the loss function optimized by KELLER which in turn yields a sparse network effectively.
To evaluate the performance of the proposed approach, we compared it with T-SWCGL and SWCGL in terms of their abilities to recover structure of dFC networks over time as well as detecting dynamically connected brain region pairs. T-SWCGL and SWCGL model dFC networks using correlation analysis at the i th window of observations. For a fair comparison, we applied graphical lasso to consider the sparsity of the estimated dFC networks. Moreover, to minimize the effect of window length in the capability of the SWCbased methods, we used various window lengths, from 20 to 140 s in 20 s steps. The simulation results suggested that KELLER was superior to T-SWCGL and SWCGL in terms of the mean F 1 score. It is notable to mention that the capability of considering multivariate dependencies between brain regions in KELLER results in the accurate estimation of dFC networks, which in turn may be the main reason of improving the detection of dynamic connections in the estimated dFC network by KELLER. As expected, modeling multivariate dependencies in estimating dFC networks increased the number of dynamic connections.

| "
The higher number of detected dynamic connections, the more statistical power of method to estimate dFC matrices," is this correct?
Answering this question is not straightforward. Actually, if there were a ground truth for dynamic behavior of brain region pairs, then this statement could be evaluated. In fact, the absence of ground truth in human brain network analysis warns us that the higher number of dynamically connected regions detected should not be interpreted as higher statistical power of a specific method in estimating dFC matrices (Savva et al., 2019;Shakil et al., 2016). However, to determine which method is the most appropriate to be applied in dFC analysis, there are two approaches: (a) designing simulation studies that provide ground truth for evaluation; and (b) comparing the results obtained on real data with those reported in the literature. Since there is no comprehensive study that investigated and reported dynamically connected brain region pairs in the whole brain, we only focused on the detected dynamically connected regions with posterior cingulate cortex (PCC) in DMN and compared the results with those of the previous studies (Chang & Glover, 2010;Savva et al., 2019). Specifically, in Savva et al., (2019), the presence of dynamic connections between brain regions comprising DMN was examined by employing various window lengths and FC metrics. They found that PCC and bilateral parietal lobes were involved in most dynamic connections as the hub regions of dorsal and ventral DMN, respectively (Savva et al., 2019).
Moreover, in Chang & Glover, (2010), dynamic connections between PCC and those brain regions that are correlated or anti-correlated with PCC were examined, utilizing wavelet transform. The authors found that phase coupling between these regions and PCC were dynamic (Chang & Glover, 2010). In Table 5, the results are shown for the common dynamic connections detected by KELLER and T-SWCGL (w = 40, 60 and 100 s), focusing on DMN (Shirer, Ryali, Rykhlevskaia, Menon, & Greicius, 2012). Most of those commonly detected dynamic connections are located between subcortical, temporal and frontal regions of DMN. On the other hand, in Table 6, we reported specific dynamic connections detected only by KELLER or T-SWCGL. For KELLER, these connections are specially located at T A B L E 6 Specific dynamic connection identified in DMN only by KELLER or T-SWCGL with different window length (w = 40, 60, and 100 s). dorsal DMN connection pairs (PCC-frontal cortex, PCC-precuneus, and precuneus-frontal cortex) which were commonly detected in the previous studies (Chang & Glover, 2010;Savva et al., 2019). As can be seen in Table 6

| Study limitations and future work
We mentioned some limitations of SWC-based methods in this study and introduced KELLER to overcome them. However, KELLER also has its own limitations. For instance, it is very time consuming due to computational complexity of the optimization algorithm and the model selection approach involved in KELLER. This is especially true when the number of involved brain regions is large. On the other hand, SWC-based methods are time efficient in estimating dFC networks. However, since these methods have low sensitivity in detecting previously reported dynamic connections (Tables 5 and 6), developing a powerful method to estimate dFC time series with high accuracy is of critical importance.
Assessment of the estimated dFC networks to detect statistically significant dynamically connected brain region pairs is highly suggested for future work in this vibrant area of neuroimaging research. Future work can also apply KELLER for estimating dFC from rs-fMRI data of neurodegenerative and neuropsychiatric disorders to investigate the abnormal dynamic connectivity patterns and compare their findings with those of the conventional methods.

| CONCLUSION
In this study, a comprehensive analysis of dFC in the whole brain using rs-fMRI data were performed by employing the newly introduced KELLER method considering multivariate dependencies between brain regions to estimate dFC network. In order to evaluate the proposed dFC estimation method in comparison with SWC, a series of simulation studies was implemented providing ground truth, and then a hypothesis testing framework was applied to detect dynamically connected region pairs. The simulation results showed that KELLER outperformed T-SWCGL and SWCGL approaches in retrieving dFC pattern of brain networks in terms of the network's structure and modulation over time as well as in detecting dynamically connected brain regions. Experimental results showed that KELLER was able to detect previously reported dynamically connected brain region pairs within DMN. Overall, dFC network analysis has a promising capability for a better understanding of the brain as well as contributing to the development of new biomarkers for diagnosis or prognosis of mental disorders. However, statistical assessment of the estimated dFC should be done as a primary step to ensure that the detected dFC patterns are due to the dynamic nature of the cognitive process or resting state condition in the brain.