Dynamic Functional Connectivity Change-Point Detection With Random Matrix Theory Inference

To study the dynamic nature of brain activity, functional magnetic resonance imaging (fMRI) data is useful including some temporal dependencies between the corresponding neural activity estimates. Recent studies have shown that the functional connectivity (FC) varies according to time and location which should be incorporated into the model. Modeling this dynamic FC (DFC) requires time-varying measures of spatial region of interest (ROI) sets. To know about the DFC, change-point detection in FC is of particular interest. In this paper, we propose a method of detecting a change-point based on the maximum of eigenvalues via random matrix theory (RMT). From covariance matrices for FC of all ROI's, the temporal change-point of FC is decided by an RMT approach. Simulation results show that our proposed method can detect meaningful FC change-points. We also illustrate the effectiveness of our FC detection approach by applying our method to epilepsy data where change-points detected are explained by the changes in memory capacity. Our study shows the possibility of RMT based approach in DFC change-point problem and in studying the complex dynamic pattern of functional brain interactions.


INTRODUCTION
Our brain is a complex network which consists of spatially distributed, but functionally linked regions that continuously share information with each other. Recent advanced functional neuroimaging studies have provided new tools to measure and to examine functional interactions between brain regions. These make the examination of DFC in the human brain and lead a development of theory associated with the complex FC structures. The fMRI data contains information about brain activity with complex spatial-temporal correlations. FC in the brain may be represented by these spatial dependence patterns and may reveal important characteristics of brain function and individual variations in cognition and behavior. FC is fundamentally a statistical concept, and is generally defined and assessed using statistical measures such as correlation (Biswall et al., 1995), and cross-coherence (Sun et al., 2004) in the temporal evolution of neural activity between distinct brain locations. Hutchison et al. (2013a,b) researched promise, issues, and interpretations as related with DFC.
Resting-state fMRI measures spontaneous fluctuations in the blood-oxygenation-leveldependent (BOLD) signal in gray matter regions. The fMRI connectivity approaches estimate statistical dependencies between gray matter activity arising from different regions (Biswall et al., 1995;Sato et al., 2006;Friston, 2011). For example, FC alterations are linked to depression, schizophrenia, and other illnesses (Mayberg et al., 1999;Jafria et al., 2008). Koike et al. (2011) provided that the default mode network (DMN) and subsystems might serve to integrate brain regions with performing function specific to each level of arousal.
Functional connectivity modelings are addressed with independent component analysis (ICA) (McKeown et al., 1998;Calhoun et al., 2001), pairwise correlation analysis (Supekar et al., 2008), and partial correlation analysis (Salvador et al., 2005). A key assumption of the aforementioned approaches is that functional connectivity remains stationary throughout the scanning session. However, many recent researches provide evidence of the dynamic nature of the brain's functional organization (Hellyer et al., 2014). Functional networks appear to fluctuate on a time scale and to reorganize under tasks during the scanning session (Chang and Glover, 2010).
Time-varying and dynamic nature of connectivity is considered in many recent neuroimaging studies. A well-known method for exploring time-varying connectivity is applying a sliding window of pre-specified length. This method calculates the correlations between distinct locations throughout the duration of the window at each step. Limitations with this sliding window approach, that the window length strongly influences the temporal smoothness or variability observed in the resulting correlations are discussed (Lindquist et al., 2014).
Identifying change-points is an important subject in neuroimaging study as it shows properties of brain networks as they relate to experimental stimuli and disease processes. Many change-point methods for fMRI signals have been proposed (Lindquist et al., 2007;Robinsona et al., 2010;Cribben et al., 2013;Xu and Lindquist, 2015). Xu and Lindquist (2015) assume that the timing of a subject's activation onset and duration are random variables, and their distributions are used to approximate the probability that a voxel/region is activated as a function of time with multi-subject change-point modeling. Bayesian approaches Kim and Cheon, 2010) and a non-parametric Fourier functional approach (Kim and Hart, 2011) are proposed to estimate change-points with distributional parameter changes which can be extended to the covariance change-point problems. Shao and Zhang (2010) study tests for change-point in time series, which is applicable to each individual fMRI series. Recently change-point detection for the multivariate covariance matrix are proposed for brain networks (Kim, 2019;Chenouri et al., 2020).
There have been many studies on testing for distributional structure changes for a sequence of independent and identically distributed (IID) random variables, including Csörgő and Horváth (1997). However test statistics developed for changepoint detection in the IID context may not work with the time series observed in neuroimaging as they do not take into account the temporal dependence in the data. Aue et al. (2009) proposed break detection in covariance based on the vech operator to summarize information in the covariance matrix. Here vech is the vectorization operator which converts the matrix into a column vector. Aue and Horváth (2013) reviewed break points estimation in covariance structure in time series. Inclán and Tiao (1994) proposed variance change detection with cumulative sums for independent series. Petersen and Muller (2016) studied functional data covariance and applied to Alzheimer fMRI data. For detecting correlation change, Wied and Krömer (2012) considered the cumulated sums of empirical correlations to make a change test.
In the present study, we aim to deal with spatial time series data to detect FC change-points based on random matrix theory. Our approach does not make any parametric assumptions. It is especially designed for detecting one change-point but is extendable to multiple change-points in FC. In addition, we applied our developed method to fMRI data of post-surgical epilepsy patients. Previous studies have shown that activity and FC of the medial temporal lobe (MTL) during the post-encoding offline period are closely related to the memory consolidation (Tambini et al., 2010;Staresina et al., 2013;Tambini and Davachi, 2013). However, how the brain supports the memory consolidation process in the absence of one MTL structure has not been investigated. Here, we applied our developed method to the post-encoding rest period to evaluate its usefulness in characterization of memory consolidation networks in MTL resected brains.
This paper is organized as follows. In section 2, we begin by looking at the covariance change-point problem with the proposed statistics. Simulation studies are shown in section 3. Application of the proposed method to the experimental epilepsy data is shown in section 4. The paper concludes with a discussion.

Statistical Problems for Dynamic Functional Connectivity
A natural unit of fMRI analysis is multi-voxel regions of the brain that change their level of activation. The fMRI data can be widely used to detect and delineate regions of the brain that change in response to specific stimuli and tasks. We represent a single scan from a subject as a 3-D rectangular lattice consisting of volume elements (voxels). From one scanning session, each voxel contains serial measures of localized brain activity called blood-oxygen-level dependent (BOLD) fMRI responses. Through the hemodynamic response (HR) process, blood releases oxygen to neurons at a greater rate than to inactive neurons. This causes a change of the relative levels of oxyhemoglobin and deoxyhemoglobin that is detectable due to magnetic susceptibility. Since the number of intracranial voxels is too large to estimate a global spatial correlation matrix including all voxel pairs, we partition the voxels into mutually exclusive neuroanatomical regions and select ROIs in order to make the ROI-level inference.
Let y ivt be an observation measured at voxel or region v from the serial fMRI BOLD responses for an ith subject. Let i = 1, . . . , N represent subjects, v = 1, . . . , d voxels (or ROIs), and t = 1, . . . , T time points. For the ith subject, d × T observed data matrix is obtained. It makes a multivariate time series for each subject and its covariance change is of interest. Tests that assess the structural stability of volatilities or cross-volatilities for multivariate non-linear time series such as fMRI data and change-point estimation is of great importance for understanding the underlying dynamics.
Let's fix the ith subject for the notational simplicity. Consider the observed y t is a d × 1 random vector at the time point t with E[ y t 2 ] < ∞. Here · denotes the Euclidean norm in R d . Based on the observations of the random vector y 1 , . . . , y T , consider a test to identify covariance matrix change for each subject.
The null hypothesis is which indicates the constancy of the covariances during the observational time period. Against this null hypothesis the alternative hypothesis with at most one change in the covariance is H A : Cov(y 1 ) = · · · = Cov(y τ ) = Cov(y τ +1 ) = · · · = Cov(y T ) (2) where τ is the change-point. Let t = Cov(y t ).

Dynamic Functional Connectivity Test With the Largest Eigenvalues
We summarize the covariance structure using maximum eigenvalues, defined as the largest eigenvalues of the covariance matrix as a functional representative measure of the matrix.
Let λ 1 (t) = max λ 1 (t), . . . , λ d (t) denote the largest eigenvalues of the matrix t . The null hypothesis in (1) can be considered as that which indicates the constancy of the maximum eigenvalues of covariances during the observation time period. The alternative hypothesis can be written with the maximum eigenvalue measure as where τ is denoted as the change-point. Suppose that the data are from d−dimensional Gaussian distribution with mean zero and positive definite covariance matrix . The sample covariance matrix (SCM) formed with these T samples as has the (central) Wishart distribution (Goodman, 1963;Raj Rao et al., 2008).
Estimating the eigenvalues of a population covariance matrix from a sample covariance matrix is a problem of fundamental importance in multivariate statistics since the eigenvalues of covariance matrices play a key role in many widely used techniques with a suitable theory (Goodman, 1963;Okamoto, 1973;Gupta and Nagar, 1999;Mohsen, 2013). The amount of variance explained is measured by the eigenvalues of the population covariance, d×d . Random matrix theory predicts that the eigenvalues of a sample covariance matrix S T are not good estimators of the eigenvalues of the population covariance when the sample size T and the number of variables d are large. El Karoui (2008) proposed a better estimator for the eigenvalues of large dimensional covariance matrices. Under the weak assumption that the Marčenko and Pastur (1967) equation holds, the El Karoui (2008) estimator can be thought as a shrinkage in a non-linear fashion.
Let λ 1 ≥ λ 2 ≥ · · · ≥ λ d be the eigenvalues of the population covariance matrix. Accordingly let l 1 ≥ l 2 ≥ · · · ≥ l d be the eigenvalues of the sample covariance matrix. Anderson (1963) showed that the case where = I d×d with all the population eigenvalues one, under some moment conditions, the maximum eigenvalue is not a consistent estimator of λ 1 and can be extremely biased when T and d are both large. To de-bias extreme sample eigenvalues, we use population spectral distribution, a probability measure that characterizes the population eigenvalues (El Karoui, 2007). The limiting distribution of the sample eigenvalues are dependent on data dependency and their covariance matrix. Johnstone (2001) derived the limiting distribution of the largest sample eigenvalue with independent Gaussian case and Bai and Silverstein (2010) dealt with this problem the sub-Gaussian case. In terms of the spectral decomposition of a covariance matrix, an estimator is rotation-equivariant if and only if the eigenvectors are the same as those of the sample covariance matrix (Mohsen, 2013). Thus, the differences between two such estimators appear only in their eigenvalues.

Dynamic Functional Connectivity Test With Tracy-Widom Distributions
Assume that y u 's are independent from d−dimensional distribution, u = 1, . . . , T. Let the sample covariance matrix up to t and after t be are the sample mean vectors separated by t. If y u 's are independent from d−dimensional normal distribution, A t and B t are d × d matrices following Wishart distributions denoted by The scale matrix has no effect on the distribution of the eigenvalues, and so consider the largest eigenvalue λ 1 (t) of (A t + B t ) −1 A t as the greatest root statistic.
If there exists change and change-point occurs at time point t, there is change in A t and B t . This change is reflected in the test statistic G t .
Under H 0 in (3), it holds that λ 1 (t * ) = · · · = λ 1 (T − t * ). For example, we can set t * = 30 to have enough degrees of freedom. We consider the following statistic to provide information about change and change-points in the covariance matrix.
Johnstone (2009) provided the appropriate centering and scaling, and the logit transform The distribution F 1 was found by Tracy and Widom (1996) as the limiting law of the largest eigenvalue of d×d Gaussian symmetric matrix. The centering and scaling parameters are given by where the angle parameters ϕ = ϕ t , γ = γ t are defined by We apply Tracy and Widom scaling and centering for each maximum eigenvalues of (A t +B t ) −1 A t and then max eigen Therefore our test is modified by Tracy and Widom correction to de-bias of the eigenvalues as follows (14) Here max eigen TW denotes that the maximum eigenvalue is modified by Tracy and Widom correction. We propose the test statistic for change as which rejects H 0 if T is large. Accordingly the proposed changepoint estimator isτ Permutation tests have the flexible class of tests for both non-parametric and parametric (model-based). The proposed statistics can be used regardless of the data distributions using this permutation approach. The threshold calculation based on block permutation can be applied for testing for change-points (Strasser and Weber, 1999;Husková, 2004). Exchangeability of the errors might be too strong of an assumption for time series applications due to the dependence structure of the observations. Block permutation is an applicable alternative. Block permutation principles are suitable for testing autoregressive series and they are refined (Kirch and Steinebach, 2006;Kirch, 2007;Zeileis and Hothorn, 2013). Approximations of critical values of the change test statistics can be obtained through block permutation methods (Husková, 2004;Luger, 2006;Kirch, 2007).
To do our testing procedure, the approximate critical values are required. For each subject with dependent time series data, each critical value can be obtained according to the block permutation tests. This block permutation tests provide asymptotically valid results in the presence of dependency.

Simulation Study
In this section we set up the simulation in order to assess and compare the performance of the proposed estimator for changepoint estimation. The simulation study is focused on changepoint estimation with multi-subject data from the multivariate vector autoregressive (MVAR) models. The objective of each simulation is to detect the functional connectivity change-point.
Before change-point estimation, testing for change is done. Since it is not possible to obtain the exact distribution of the test statistic T analytically, it is determined by simulation. The data are synthetically generated from multivariate normal distribution N(0, d ) for the null distribution with the given dimension d as the number of ROIs, and the number of time points T. Several types of d are designed. For the comparison we also consider the Aue et al. (2009) method. Aue et al. (2009) proposed the test for covariance change based on Their change test statistic is whereˆ T is a consistent covariance estimator. Aue et al. (2009) change-point estimator iŝ We useˆ T as the sampled version of the moment estimator. The asymptotic critical values of the test statistics are computed from the block permutation based empirical distribution.
We simulate data according to MVAR(1) and MVAR(2) models. with the regional dimension d = 5 and T = 200 time points in 1,000 repetitions. Since the brain network is highdimensional, we also consider d = 20. The empirical critical values under the null distribution are obtained from no changepoint model using block permutations for each subject. The fraction of one change-point is considered as θ = τ/T = 0.5 and 0.7 in the simulated data. We calculate the performance measure such as mean, median, square root of MSE (mean squared error), and the proportion around the true change-point as p5 =P(|τ − τ | < 5) with 95 % confidence level.
Data are generated from the following model, y u ∼ MVAR(p) as where B k is kth parameter coefficient matrix of MVAR(p), k = 1, . . . , p. And 1 and 2 are covariance matrix of ǫ u 's before and after the change-point respectively. We consider MVAR(1) models in cases (i) and (ii), and MVAR(2) models in (iii) and (iv) in the followings.

I. MVAR(1) with coefficient matrix B
where I is the identity matrix and J is a matrix whose components are all 1's.
II. MVAR(2) with coefficient matrices B 1 and B 2 B 1 = 11 = [(0.2)1, 0, · · · , 0], u ≤ τ 12 = 11 + (0.1)I, u > τ , To illustrate the potential to detect multiple change-points, we can apply the procedure with the subsegment after the changepoint. Each time we reject H 0 , we re-apply the method to the subsamples splitted by the proposed change-point. We considered MVAR(1) and MVAR(2) cases with the different covariance structure after the change-point at θ = 0.5, 0.7. Tables 1, 2 provide MVAR(1) cases in which our estimators have smaller MSE's and higher matching proportions than Aue estimators. Tables 3, 4 give the change-point estimation results in MVAR(2) cases. Our method performs better than the Aue method in Tables 3, 4. Aue method depends on the covariance estimatorˆ T in the statistic in (19), and the change-point estimation results depend on the covariance structure accordingly. But our change-point estimator gives stable simulation results in the performance based upon mean and MSE. Overall in this simulation the proposed estimator performs very well.

Epilepsy fMRI Data
The same subjects from our previous study are used in this study except for one patient who showed a severe movement artifact (Jeong et al., 2019). Patients who underwent unilateral MTL resection (MTLR) to treat medically intractable temporal lobe epilepsy at Seoul National University Hospital at least 1 year before recruitment and were between 19 and 50 years of age are retrospectively recruited. We only include subjects who showed at least a low average or a higher level (scores> 80) of memory capacity and general intelligence evaluated by the standard neuropsychological test. A total of 34 patients (16 left and 18 right; median age = 32.5 years) and 24 age-and educationyear-matched healthy controls (HC, median age = 32 years) are included in the present study. All subjects provided informed consent. This study was approved by the Institutional Review Board of Seoul National University Hospital.
The MR images were acquired on a research-dedicated 3T Magnetom Trio Tim Syngo (Siemens, Erlangen, Germany) using a 32-channel head coil. Five minutes of resting-state functional data (eyes open, fixation cross) were acquired both before (preencoding baseline) and after performing the in-scanner memory encoding paradigm of words and figures (post-encoding) using a T2 * -weighted gradient echo planar imaging sequence (36 axial slices, slice thickness = 3.4 mm, no gap, TR = 2,000 ms, TE = 30 ms, FOV the change-point 220 × 220 mm, flip angel = 80 • , voxel size = 3.4 × 3.4 × 3.4mm 3 , and interleaved). Whole-brain high-resolution anatomic T1-weighted images were obtained with the 3D TFL sequence (TR = 1,670 ms, TE = 1.89 ms, FOV = 250 × 250H mm, flip angle = 9 • , voxel size = 1.0 × 1.0 × 1.0mm 3 ). The resting-state functional data underwent a number of preprocessing steps including motion correction, slice time correction, co-registration, spatial normalization, and spatial smoothing (AFNI, version: 16.0.00, https://afni.nimh.nih. gov/afni/). More details about the memory encoding paradigm and data preprocessing have been described elsewhere (Jeong et al., 2018(Jeong et al., , 2019.  Since previous studies have shown that FC during the postencoding resting-state is related to the memory performance (Tambini et al., 2010;Staresina et al., 2013;Tambini and Davachi, 2013), we thought that data that showed evident FC changes of memory related areas rather than entire data of post-encoding resting-state may possibly provide the additional information of the memory consolidation network. In order to find the point that showed evident FC changes, we used change-point analysis.
For change-point estimation analysis, we select 20 ROIs within DMN, which is a well-known network that subserves memory function (Andrews-Hanna et al., 2010;Jeong et al., 2015). Table 5 shows the details of selected ROIs. In both preand post-encoding rest data, BOLD signals were extracted within each selected ROI (8 mm radius). For patients, one ROI, left hippocampus for left MTLR (LMTLR) and right hippocampus for right MTLR (RMTLR), was excluded from BOLD time  series extraction. The extracted time series are then used for change-point estimation of post-encoding rest data to detect changes of FC after memory encoding. Next, pair-wise Pearson correlation coefficients between each pair of ROIs are calculated and averaged for each ROI, separately for time series after change-point and pre-encoding. FC is additionally calculated for the whole-segment of post-encoding time series in order to evaluate the usefulness of the change-point estimation method in the detection of memory encoding-related FC changes.

Change Analysis for Epilepsy Data
Change analysis starts with an individual test for existence of change. Under the significance of change, change-point estimation can be accomplished for its identification. With the epilepsy data, we performed the block permutation method to obtain the critical values. For dependent data like fMRI data, block permutation is preferable and useful (Kirch and Steinebach, 2006;Kirch, 2007;Adolf et al., 2011). This block permutation procedure is done individually because the change-point should be estimated subject-wise. The block permutation critical values are obtained with the block size K = 5 in 1,000 repetitions. At the significance level α = 0.05 we reject the null hypothesis for each individual subject except hc22 and lt16 (significant at α = 0.10).
Since we focus on detectable change-point estimation, we get and use the change-points for all subjects. The change-point is estimable and informative for change, so we include these two subjects for the analysis. Therefore, we identify one change-point for each subject. Table 6 shows a summary of change-point estimation and Figure 1 shows the change-points of individual subjects. The experiments are conducted during the two resting states. Changepoint estimation is done during the second state because the first resting state is considered the baseline. Change-points are significantly different among groups [F (2,55) = 4.29, p < 0.05, ANOVA]. Bonferroni post-hoc tests reveal that the LMTLR group shows significantly higher change-point values than the HC (p < 0.05) and RMTLR groups (p < 0.05). There is no significant differences between the HC and RMTLR groups (p > 0.10).
Changes of FC after change-point and pre-encoding baseline are compared within each subject group by using the Wilcoxon signed rank test (Figure 2). FC without considering change-point estimation (whole-segment of post-encoding period) are also compared with baseline FC. FC values in areas that survived are presented in Table 7. In the HC group, the FC of the right parahippocampal cortex (PHC) shows a tendency to increase after the change-point (p = 0.067) and significantly increased in whole-segment of post-encoding rest (p < 0.05). In RMTLR group, none of ROIs shows significant FC differences. In the LMTLR group, although not statistically significant, the FC of the right temporal pole (TempP) shows a tendency to decrease after change-point (p = 0.088). In this ROI, wholesegment comparisons shows significant difference (p < 0.05). Interestingly, the FC of the right posterior inferior parietal lobule (pIPL) is significantly decreased after the change-point (p < 0.05) but shows no differences with the whole-segment of post-encoding rest (p = 0.642). Moreover, the FC difference (after change-point -baseline) in the right pIPL is negatively correlated with verbal immediate memory scores of the standard neuropsychological Rey-Kim memory test in LMTLR group (r = −0.548, p = 0.028; Spearman's correlation; Figure 2. There are no other significant clinical correlations of FC, to memory scores in either after-change point or whole-segment comparisons. The relationships between the behavioral accuracy (d-prime) of the in-scanner memory task and the FC values of postencoding rest after change-point are also analyzed. Our previous study describes behavioral accuracy in detail (Jeong et al., 2019). Figure 3 shows the significant correlations. The FC of the ventral

DISCUSSION
In this paper we develop a new procedure that can be used to determine a change-point in a individual FC system, using a function of maximum eigenvalues of suitably scaled covariance matrices. The proposed method is applicable to either resting-state or task fMRI studies providing individualized change-point estimation. The method yields individualized change-point estimation results, but it permits to group-level estimation. Our change-point estimation method generally works well in simulations and shows improved performance over the previous approach. Our work avoids the use of a sliding window usually applied to correlation measures, and reinforces the notion that time-varying techniques are required to study FC. Our results stress the importance of taking the change-point into account to characterize dynamic FC in order to describe how "information" is integrated and changed over time.
In this study, we applied change-point estimation analysis to post-encoding rest fMRI data of MTLR patients in order to evaluate the applicability and usefulness of the proposed method. We found that both segment after change-point and whole-segment of resting-state data after memory encoding detect similar FC changes. In HC, the MTL area (PHC) shows increased FC after memory encoding which is consistent with previous findings of memory consolidation studies (Tambini et al., 2010;Staresina et al., 2013). In LMTLR patients, the FC of the right pIPL and TempP decreased after memory encoding. Importantly, right pIPL is only detected when we use postencoding data after the change-point, but not detected with the whole-segment of data. Moreover, the magnitude of FC changes in the right pIPL is associated with verbal memory capacity of LMTLR patients. Patients who displayed more reduction in FC values have better verbal memory. Recent studies suggest that the posterior parietal cortex contributes to memory consolidation (Brodt et al., 2016(Brodt et al., , 2018. Therefore, we speculate that DMN areas other than damaged MTL, especially the right pIPL, may play a compensatory role for the verbal memory function of resected left MTL. In summary, by applying change-point analysis, we can reveal additional information on the memory consolidation network in an MTL resected brain. Further studies with other data sets are warranted to ensure the general use of this changepoint estimation analysis of fMRI data.
The other notable contribution of our work is that we develop a generalized multivariate extension to use an eigenvalue system based on RMT. While previous approaches have largely targeted univariate approaches examining a single time series or a single measure of association, our method use full information of connectivity via the eigenvalue system. Another interesting aspect of our change-point detection is the quantitative comparison before and after change-point. It reveals some properties of FC after change occurs. Different analytic approaches provides different perspective information on fMRI data, which necessitates integrating knowledge gained through the variety of models for the understanding of dynamic FC.
One limitation of the method is that the simulation results depend on the underlying covariance structure. When there are more change-points, simultaneous multiple change-points estimation procedure should be used. Further study with this extension is expected when there are multiple change-points in FC.

DATA AVAILABILITY STATEMENT
The datasets generated for this article are not readily available because: it can be available only if requested to the third author.
Requests to access the datasets should be directed to Chun Kee Chung, chungc@snu.ac.kr.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board of Seoul National University Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
Statistical method was developed and simulation work was done by JK. Epilepsy data were collected by CC. Epilepsy data analyses and interpretations were conducted by WJ. All authors contributed to the article and approved the submitted version.