Inference of Dynamic Graph Changes for Functional Connectome

TL;DR

This paper introduces a new statistical method to detect and analyze dynamic changes in brain networks over time using high-dimensional graphical models, enhancing understanding of brain responses to stimuli.

Contribution

A novel hypothesis testing framework using bootstrap statistics for detecting change points in dynamic brain connectivity networks.

Findings

Successfully identified change points in simulated data.

Detected meaningful temporal brain network changes during stimuli.

Correlated network changes with emotional annotations.

Abstract

Dynamic functional connectivity is an effective measure for the brain's responses to continuous stimuli. We propose an inferential method to detect the dynamic changes of brain networks based on time-varying graphical models. Whereas most existing methods focus on testing the existence of change points, the dynamics in the brain network offer more signals in many neuroscience studies. We propose a novel method to conduct hypothesis testing on changes in dynamic brain networks. We introduce a bootstrap statistic to approximate the supreme of the high-dimensional empirical processes over dynamically changing edges. Our simulations show that this framework can capture the change points with changed connectivity. Finally, we apply our method to a brain imaging dataset under a natural audio-video stimulus and illustrate that we are able to detect temporal changes in brain networks. The functions of the identified regions are consistent with specific emotional annotations, which are closely associated with changes inferred by our method.

Figures

Tables

Table 1. Table 1: Type I error rate

	$N_T$=$700$	$800$	$900$	$1000$
$p$=50	0.04	0.04	0.04	0.05
100	0.04	0.05	0.04	0.04

Table 2. Table 2: Power of detecting change points with α = 0.2 𝛼 0.2 \alpha=0.2

	$N_T$=1000	1500	2000	2500	3000
dGCI	0.40	0.59	0.76	0.89	0.95
DCD	0.04	0.27	0.59	0.85	0.91
covcp	0.26	0.17	0.19	0.24	0.37

Table 3. Table 3: Average senstivity

	$N_T$=1000	1500	2000	2500	3000
a=0.2	0.06	0.08	0.09	0.12	0.15
0.3	0.23	0.31	0.40	0.46	0.52
0.4	0.28	0.36	0.44	0.50	0.54

Table 4. Table 4: Average false positive rate

	$N_T$=1000	1500	2000	2500	3000
a=0.2	0.18	0.09	0.04	0.02	0.02
0.3	0.05	0.03	0.02	0.01	0.01
0.4	0.06	0.03	0.03	0.02	0.02

Equations

$$\hat{\mathbf{\Sigma }}(t)=\frac{{\sum }_{i\in [n]}{K}_h(t_i-t){\text{X}}_i{\text{X}}_i^T}{{\sum }_{i\in [n]}{K}_h(t_i-t)}$$

$$K(u)=0.75\cdot (1-u^2)\text{mathbbm}{1}_{|u|\le 1},$$

$$\displaystyle\widehat{\mathbf{\Theta}}_{jk}^{d}(t)=\widehat{\mathbf{\Theta}}_{jk}(t)-\frac{\big{(}\widehat{\mathbf{\Theta}}_{j}(t)\big{)}^{T}\big{[}\widehat{\mathbf{\Sigma}}(t)\widehat{\mathbf{\Theta}}_{k}(t)-\textbf{e}_{k}\big{]}}{\big{(}\widehat{\mathbf{\Theta}}_{j}(t)\big{)}^{T}\widehat{\mathbf{\Sigma}}_{j}(t)},$$

$$U_E$$

$$\begin{array}{rl}& H_0:\forall t\in (0,1),G^{\text{+}}(t)=G^{\text{-}}(t),\\ & H_1:\exists t_0\in (0,1),(j,k)\in E(t)\text{ s.t. }{G}_{jk}^{\text{+}}(t_0)\ne G_{jk}^{\text{-}}(t_0),\end{array}$$

$$\displaystyle\widehat{\mathbf{\Theta}}_{jk}^{d\pm}(t)=\widehat{\mathbf{\Theta}}_{jk}^{\pm}(t)-\frac{\big{(}\widehat{\mathbf{\Theta}}_{j}^{\pm}(t)\big{)}^{T}\big{[}\widehat{\mathbf{\Sigma}}^{\pm}(t)\widehat{\mathbf{\Theta}}^{\pm}_{k}(t)-\textbf{e}_{k}\big{]}}{\big{(}\widehat{\mathbf{\Theta}}^{\pm}_{j}(t)\big{)}^{T}\widehat{\mathbf{\Sigma}}^{\pm}_{j}(t)},$$

$$\displaystyle\begin{split}U_{jk}^{\pm}(t)&=\widehat{\mathbf{\Theta}}_{jk}^{d\pm}(t)-\mathbf{\Theta}_{jk}^{\pm}(t),\\ U_{E}&=\sup_{t\in(0,1)}\max_{(j,k)\in E(t)}\sqrt{nh}\big{|}U_{jk}^{\texttt{+}}(t)-U_{jk}^{\texttt{-}}(t)\big{|},\end{split}$$

$$\displaystyle M_{ijk}^{B\pm}\!(t)\!=\!K_{h^{\pm}}(t_{i}\!-\!t)\!\cdot\!\big{(}\!\widehat{\mathbf{\Theta}}_{j}^{\pm}(t)\!\big{)}^{T}\!\big{[}\mathbf{X}_{i}\mathbf{X}_{i}^{T}\!\widehat{\mathbf{\Theta}}^{\pm}_{k}(t)\!-\!\mathbf{e}_{k}\big{]}$$

$$U_{jk}^{B\pm }(t)$$

$$U_E^B$$

$$\begin{array}{r}\inf \{q\in \mathbb{R}|P(U_E^B\le q|\{(t_i,{\mathbf{X}}_i){\}}_{i\in [n]})\ge 1-\alpha \}.\end{array}$$

$$\displaystyle\mathcal{R}(t)\!=\!\big{\{}(j,k)\!\in\!E(t)\big{|}\sqrt{nh}\cdot\big{|}\Delta\widehat{\mathbf{\Theta}}_{jk}(t)\big{|}>c(1-\alpha,\mathbf{E})\big{\}},$$

$$\displaystyle\frac{\sum_{i=1}^{N}K_{h^{\pm}}(t_{i}-t)\Big{(}\big{(}\widehat{\mathbf{\Theta}}_{j}^{\pm}(t)\big{)}^{T}\big{[}\mathbf{X}_{i}\mathbf{X}_{i}^{T}\widehat{\mathbf{\Theta}}^{\pm}_{k}(t)-\mathbf{e}_{k}\big{]}\Big{)}^{2}}{\sum_{i=1}^{N}K_{h^{\pm}}(t_{i}-t)}.$$

$$U_{E}=\sup_{t\in(0,1)}\max_{(j,k)\in E(t)}\sqrt{nh}\cdot\Big{|}\delta\widehat{\mathbf{\Theta}}_{jk}(t)\Big{|}$$

$$U_{E}^{B}\!=\!\sup_{t\in\!(0,1)}\!\max_{(j,k)\in\!E(t)}\!\sqrt{nh}\cdot\tilde{\sigma}^{-1}_{jk}\!(t)\bigg{|}U_{jk}^{B\texttt{+}}\!(t)-U_{jk}^{B\texttt{-}}\!(t)\bigg{|}$$

$$\begin{array}{rl}{\mathcal{U}}_{s,m,M}& =\{\Theta \in {\mathbb{R}}_{p\times p}|\Theta \succ 0,{\lambda }_1(\Theta )>m,\\ & \Vert \Theta {\Vert }_2\le \rho ,\underset{j\in [p]}{\max }\Vert {\Theta }_j{\Vert }_0\le s,\underset{j\in [p]}{\max }\Vert {\Theta }_j{\Vert }_1\le M\},\end{array}$$

$$\begin{array}{rl}\underset{t\in [0,1]}{\sup }\underset{j,k\in [p]}{\max }\max \{& {\Sigma }_{jk}^{+}(t),{\dot{\Sigma }}_{jk}^{+}(t),{\ddot{\Sigma }}_{jk}^{+}(t),\\ & {\Sigma }_{jk}^{-}(t),{\dot{\Sigma }}_{jk}^{-}(t),{\ddot{\Sigma }}_{jk}^{-}(t)\}\le M_{\sigma }.\end{array}$$

$$\underset{t\in [a,b]}{\sup }\Vert {\hat{\mathbf{\Sigma }}}^{\pm }(t)-{\mathbf{\Sigma }}^{\pm }(t){\Vert }_{max}\le C\cdot \left(h+\sqrt{\frac{\log (np/\sqrt{h})}{nh}}\right)$$

$$\underset{n\to \infty }{\lim }\underset{\mathbf{\Theta }(\cdot )\in {\mathcal{U}}_{s,m,M}}{\sup }{P}_{\mathbf{\Theta }(\cdot )}\left(U_E\ge c\left(1-\alpha ,E\right)\right)\le \alpha .$$

Full text

1 INTRODUCTION

The time series of multivariate variables with underlying time-varying graphical structures are commonly studied in various areas such as transcriptomics (Rodius et al.,, 2016; Willsey et al.,, 2013), finance (Silva et al.,, 2015; Isogai,, 2017), and neuroscience (Kabbara et al.,, 2017; Soreq et al.,, 2019). With time-evolving gene co-expression networks, dynamic stock price networks, and dynamic brain connectivities in cognitive processes, these studies have shown the necessity of investigating dynamic changes in graph structures, especially when temporal changes or specific turning events are expected.

1.1 Motivation

The turning events during a time course can be abrupt, and bring sudden perturbations to the system. For example, bifurcation exists in important biological regulatory networks such as the p53 system (Hat et al.,, 2016), where a small change or perturbation can cause a qualitative change of the state in the system (Champneys and Tsaneva-Atanasova,, 2013). Structural breaks may be observed in other data types with a similar design or sudden changes caused by unexpected events, such as financial data during bursts of emergencies (Pepelyshev and Polunchenko,, 2015).

Brain imaging is also a typical example of signals with graphical variety during a certain external stimulus, and the brain transits to distinct network states according to task demands (Gonzalez-Castillo and Bandettini,, 2018). The complex network map of the human brain, which is defined as the human connectome by Sporns et al., (2005), is a collection of connections in the brain. It consists of scattered voxels or regions of interest (ROIs) and their entangled spatiotemporal relationships, including structural and functional connectivity (Contreras et al.,, 2015). Functional connectivity (FC) is inferred using neural activity signals such as BOLD sequences (Smith et al.,, 2013). Specific statistical measures such as correlation and mutual information are commonly used to characterize the brain network structure (Jeong et al.,, 2001), where a stationary interaction structure is assumed in FC calculation. However, fluctuations of both BOLD sequences and FC have become of great interest with the advancement of technology and improved experimental designs (Preti et al.,, 2017). Calhoun et al., (2014) extended the connectome concept with a dynamic perspective, namely “dynamic connectome”, to describe a spectrum of approaches aiming to identify time-varying properties of dynamic functional connectivity (dFC).

1.2 Related Works

The applications of a number of dynamic graph estimation methods developed under different assumptions have successfully identified biologically meaningful dynamic patterns for different brain states (Betzel et al.,, 2016; Gonzalez-Castillo et al.,, 2015; Wang et al.,, 2013). A sliding window approach is commonly adopted to link functional brain dynamics and cognitive abilities where the functional connectivity at a specific stage is constructed with observations in the corresponding time window (Elton and Gao,, 2015; Fong et al.,, 2019). Using the sliding window approach, Sakoğlu et al., (2010) revealed different dynamic connectivities in schizophrenia in terms of task-modulation. Although the sliding window approach is straightforward and easy to interpret, the choice of the window length is still subject to debate (Hutchison et al.,, 2013; Preti et al.,, 2017). Lindquist et al., (2014) proposed multivariate volatility models to refine the sliding window methods with exponential weights and used dynamic conditional correlation to measure dFC (Preti et al.,, 2017). Chang and Glover, (2010) applied time-frequency analysis to circumvent the window limitation, and extended the dynamic connectivity to multiple frequencies.

Furthermore, with the advancement of change point detection methods, BOLD sequences can be further divided into segments for a better understanding of network structure within each segment (Xu and Lindquist,, 2015). A general comparison between time segments can be achieved by a summary statistic, for example, a likelihood ratio statistic (Barnett and Onnela,, 2016), where a multivariate Gaussian distribution is assumed to facilitate statistical inference. However, it is hard to impose specific hypotheses on the heterogeneity of the graphs along a time course because changes can be continuous, which invalidates the time segments proposed by the likelihood ratio test. In (Park et al.,, 2018), a two-sample test based on covariance matrix estimators on different time windows was proposed to detect the change time points of brain connectivity. Avanesov et al., (2018) also addressed the multi-scale properties in their testing framework, which directly tests the difference of covariance structures. Though the testing framework is similar, our method is different from the existing literature in three aspects: 1. We assume continuity instead of structural invariance for the null hypothesis, with a focus on discontinuity with sudden changes. 2. We consider the temporal dependence of the covariance structures by applying a time kernel. 3. We accordingly apply the Gaussian multiplier bootstrap to approximate the distribution of the test statistic to detect potential changes over the whole time course, and we prove that the testing procedure appropriately controls the type I error rate.

1.3 Time-varying Gaussian Graphical Models

Different testing statistics may be applied under similar frameworks. Given a time series X in a temporal scale observed on $[0,T]$, which is commonly assumed as Gaussian at a specific time point, a test statistic can be constructed for the null hypothesis where there is no change in the distribution or the covariance structure over the time course. The null hypothesis is then rejected when the test statistic derived from the observed time series reaches certain statistical significance. We adapt the inter-subject time-varying Gaussian graphical models proposed by Tan et al., (2019) to consider individuals instead of two representative time series for the estimation. The following is the setting of the inferential framework. Let $\textit{{X}}\in\mathbb{R}^{p}$ be a $p$-dimensional random variable and define a time variable $\text{T}\in[0,1]$. In the time-varying Gaussian graphical model, we assume $\textit{{X}}|\text{T}=t\sim N_{p}(\mathbf{0},\mathbf{\Sigma}(t)),$ where the covariance matrix $\mathbf{\Sigma}(t)$ is defined at $\text{T}=t$ and the corresponding precision or inverse covariance is $\mathbf{\Theta}(t)=(\mathbf{\Sigma}(t))^{-1}$ which governs the conditional dependencies for the graph. The inverse covariance matrix $\mathbf{\Theta}(t)$ is regarded as a representation of the time-varying graph $G(t)$. Thus, graphs across time can be retrieved by estimating $\mathbf{\Theta}(t)$. A kernel smoothed covariance estimator (Yin et al.,, 2010; Zhou et al.,, 2010) is used to obtain $\widehat{\mathbf{\Sigma}}(t)$:

[TABLE]

where $K_{h}(t_{i}-t)=K((t_{i}-t)/h)/h$ with $h>0$ as the bandwidth parameter. For simplicity, we use the Epanechenikov kernel described in Tan et al., (2019), i.e.,

[TABLE]

which is commonly used in non-parametric estimation (Epanechnikov,, 1969). The precision matrix can then be estimated using the CLIME estimator (Cai et al.,, 2011), i.e., $\widehat{\mathbf{\Theta}}_{j}(t)=\operatorname*{argmin}_{\bm{\theta}\in R^{p}}||\bm{\theta}||_{1}$, subject to $||\widehat{\mathbf{\Sigma}}(t)\bm{\theta}-\textbf{e}_{j}||_{\infty}\leq\lambda$, where $\textbf{e}_{j}$ is the $j$-th canonical basis in $R^{p}$ for each $j\in\{1,...,p\}$ and $\lambda>0$ is a tuning parameter to control the sparsity of $\widehat{\mathbf{\Theta}}_{j}$. Instead of taking precision as the parameter of interest, Neykov et al., (2018) derived a de-biased estimator for the CLIME estimator

[TABLE]

which is proven to be asymptotically normal (Neykov et al.,, 2018). The construction of the de-biased precision estimator is able to circumvent the asymptotic non-normality of traditional test statistics such as score test (Neykov et al.,, 2019).

Based on (2), Tan et al., (2019) proposed a combinatoric inferential framework to test monotone graph property for time-varying graphs. They introduced a test statistic and a bootstrap method for hypothesis testing with an appropriate type I error control. In their work, Tan et al., (2019) defined the test statistic for the graph as:

[TABLE]

Based on the maximal statistic as defined in (3), we replace precision components with the difference of limits to build the test statistic for hypothesis testing. The bootstrap statistic is adapted accordingly, and the testing procedure remains the same as that in Tan et al., (2019) except for the skip-down operation. Finally, we introduce a testing framework called dynamic Graphical Change Inference (dGCI), which aims to test whether there is any sudden change in the covariance structure in the observed time series.

1.4 Our Contributions

The major contributions of our current work are:

1. We extend the null hypothesis to allow continuous changing patterns.
2. We propose a generalized testing framework for change point detection of dynamic covariance structure.
3. Besides change point detection, we introduce a powerful tool to identify connectivity changes with adequate sample sizes.

2 METHODS

In this section, we introduce the time-varying Gaussian graphical models and propose testing procedures for change point detection. We adapt model settings in Tan et al., (2019) and design a novel test statistic for detecting change points. We also propose a normalized estimate of the difference between precision matrices to enhance the performance of our method.

2.1 Test Statistic for Graph Change

Here, we formulate the testing problem and propose a corresponding bootstrap method for statistical inference adapted from (Tan et al.,, 2019). Furthermore, we apply normalization to de-biased estimates and use normalized estimates to derive the test statistic.

2.1.1 Inference of Graph Change

Given the precision matrix $\mathbf{\Theta}(t)$, we note that its left and right values at $t_{0}$ are defined as $\mathbf{\Theta}^{\texttt{-}}(t_{0})\!\coloneqq\!\lim_{t\rightarrow t_{0}-}\mathbf{\Theta}(t)$ and $\mathbf{\Theta}^{\texttt{+}}(t_{0})\!\coloneqq\!\lim_{t\rightarrow t_{0}+}\mathbf{\Theta}(t)$. We define $\mathbf{\Sigma}^{\texttt{+}}(t)$, $\mathbf{\Sigma}^{\texttt{-}}(t)$, $G^{\texttt{+}}(t)$ and $G^{\texttt{-}}(t)$ based on $\mathbf{\Sigma}(t)$ and $G(t)$ similarly. We say edge $(j,k)$ is changed at time point $t_{0}$ if the left and right values of $\mathbf{\Theta}_{jk}(t)$ at $t_{0}$ are different, i.e., $\mathbf{\Theta}_{jk}^{\texttt{+}}(t)\neq\mathbf{\Theta}_{jk}^{\texttt{-}}(t)$. Here, we do not assume temporal segments of constant precisions, which is different from the common null hypothesis with a static covariance structure (Avanesov et al.,, 2018; Xu and Lindquist,, 2015).​ For $G(t)\!=\!(V,E(t))$, the hypothesis testing of graph changes can be formulated as:

[TABLE]

Unlike a general graph property $\mathcal{P}$ defined by Neykov et al., (2019) and Tan et al., (2019), this specific null hypothesis does not require a skip-down procedure because the rejection of a single edge will reject the null at the same time. It simplifies the inference although we can still perform the skip-down approach for specific changes of properties. In order to observe changes of $\mathbf{\Theta}(t)$ at $\text{T}=t$, we need to consider both sides of $t$, i.e., $\mathbf{\Theta}^{\texttt{+}}(t)$ and $\mathbf{\Theta}^{\texttt{-}}(t)$, which represent the right and left limits, respectively. If there is no change in the current graph, it is expected that $|\mathbf{\Theta}^{\texttt{+}}(t)-\mathbf{\Theta}^{\texttt{-}}(t)|=0$.

To perform the statistical test at $\text{T}=t$, first, we generalize (1) into right and left approximated estimates. We define the right-sided covariance estimator as (1) with $K_{h}$ = $K_{h^{+}}$, where $K_{h^{+}}(u)=2K_{h}(u)\mathbbm{1}\{u>0\}$ to ensure $\int K_{h^{+}}(u)du=1$. We also define the left-sided covariance estimator $\widehat{\mathbf{\Sigma}}^{-}(t)$ similarly with $K_{h^{-}}(u)=2K_{h}(u)\mathbbm{1}\{u<0\}$.

Following the same procedure above, $\widehat{\mathbf{\Theta}}^{\texttt{+}}(t)$ and $\widehat{\mathbf{\Theta}}^{\texttt{-}}(t)$ are estimated by GLASSO (Friedman et al.,, 2008) or CLIME (Cai et al.,, 2011) using $\widehat{\mathbf{\Sigma}}^{\texttt{+}}(t)$ and $\widehat{\mathbf{\Sigma}}^{\texttt{-}}(t)$, respectively. Similar to (2), the de-biased estimators of approximated estimates for the two sides are

[TABLE]

where the notation $\pm$ means we use the left or right estimates respectively in the definition of $\widehat{\mathbf{\Theta}}_{jk}^{d+}(t)$ or $\widehat{\mathbf{\Theta}}_{jk}^{d-}(t)$. Here, we extend (3) with left- and right-sided estimates for characterizing graph changes

[TABLE]

where $\mathbf{\Theta}_{jk}^{\texttt{-}}(t)$ and $\mathbf{\Theta}_{jk}^{\texttt{+}}(t)$ represent the left- and right-sided values of $\mathbf{\Theta}_{jk}(t)$, respectively. We expand the Gaussian multiplier bootstrapper as

[TABLE]

where $t\in(0,1)$, $i\in[n]$ and $(j,k)\in E(t)$. Similar to (3) and (5), we can replace the estimates of the bootstrap statistics defined by Tan et al., (2019), which are not side-specific, with the differences of the side-specific ones as shown in (7). With two components defined in (6), the bootstrap statistic is defined as

[TABLE]

where $\xi_{i},...,\xi_{n}\overset{i.i.d.}{\sim}N(0,1)$. $c(1-\alpha,\mathbf{E})$, the conditional $(1-\alpha)$-quantile of $U_{E}^{B}$ given $\{(t_{i},\mathbf{X}_{i})\}_{i\in[n]}$ is defined as:

[TABLE]

In our implementation, we uniformly sample $t\in(0,1)$ and approximate $U_{E}^{B}$ by taking the supreme over sampled $\mathbf{t}_{B}$ instead of $(0,1)$. $U_{E}^{B}$ is sampled with $\bm{\xi}$ as defined in (7) and (8). Similar to Tan et al., (2019), the rejected edge set at $t$ is defined as

[TABLE]

where \Delta\widehat{\mathbf{\Theta}}_{jk}(t)=\big{(}{\widehat{\mathbf{\Theta}}_{jk}^{d\texttt{+}}(t)-\widehat{\mathbf{\Theta}}_{jk}^{d\texttt{-}}(t)}\big{)}. The null hypothesis will be rejected if there exists a $t_{0}\in(0,1)$ such that $\mathcal{R}(t_{0})\neq\varnothing$.

2.1.2 Approximated Variance of Test Statistic

With arbitrary changes of the graph, the variance of the de-biased estimator can vary. To avoid the dominance of unwanted variance, we normalize the estimates with approximated variance. We approximate the variance of the de-biased estimator \tilde{\sigma}^{2}\big{(}\widehat{\Theta}_{jk}^{d\pm}(t)\big{)} by taking the kernel weighted average of squares:

[TABLE]

The normalization term $\tilde{\sigma}_{jk}(t)$ for $(j,k)\in\text{V}$ is defined as \tilde{\sigma}_{jk}(t)=\sqrt{\tilde{\sigma}^{2}\big{(}\widehat{\Theta}_{jk}^{d\texttt{+}}(t)\big{)}+\tilde{\sigma}^{2}\big{(}\widehat{\Theta}_{jk}^{d\texttt{-}}(t)\big{)}}. With \delta\widehat{\mathbf{\Theta}}_{jk}(t)=\big{(}{\widehat{\mathbf{\Theta}}_{jk}^{d\texttt{+}}(t)-\widehat{\mathbf{\Theta}}_{jk}^{d\texttt{-}}(t)}\big{)}/{\tilde{\sigma}_{jk}(t)}, we normalize the test statistic

[TABLE]

and the bootstrapper as

[TABLE]

The definition of the quantile $c(1-\alpha,\mathbf{E})$ defined in (9) remains the same.

2.2 Hypothesis Testing

We summarize the bootstrap procedure and the hypothesis test in Algorithm 1 and Algorithm 2.

3 THEORETICAL RESULTS

In this section, we establish the uniform rates of convergence for the kernel smoothed covariance matrix estimator $\widehat{\mathbf{\Sigma}}^{\pm}(t)$ and the inverse covariance estimator $\widehat{\mathbf{\Theta}}(t)$ in CLIME. In addition, we show that our testing procedure in Algorithm 2 is a uniformly valid test. We study the asymptotic regime in which both $n$ and $p$ are allowed to increase. Our theoretical results are based on the results in Tan et al., (2019), and we address the differences here and in the supplementary notes. To detect “sudden changes”, we don’t assume that the covariance matrix $\Sigma\left(\cdot\right)$ is continuous. Consequently, we obtain different rates of convergence from those in Tan et al., (2019). Throughout our proof, we consider the parameter space

[TABLE]

where $0<\lambda_{1}(\Theta)\leq\lambda_{2}(\Theta)\leq\ldots\leq\lambda_{p}(\Theta)$ are the eigenvalues of $\Theta$. A similar class of matrices were considered in the literature on inverse covariance matrix estimation (Cai et al.,, 2011). We allow $s$ to increase with $n$ and $p$. Next, we impose some conditions on the regularity and smoothness of the marginal density of $t$, $f_{T}(\cdot)$ and covariance matrix $\Sigma(\cdot)$.

Assumption 1

Assume that there exists a constant $\underline{f}_{T}$ such that $\inf_{t\in[0,1]}f_{T}(t)\geq\underline{f}_{T}>0$. Furthermore, assume that $f_{T}$ is twice continuously differentiable and that there exists a constant $\overline{f}_{T}<\infty$ such that $\max\{\|f_{T}\|_{\infty},\|\dot{f}_{T}\|_{\infty},\|\ddot{f}_{T}\|_{\infty}\}\leq\overline{f}_{T}$.

Now, we introduce Assumption 2. Different from Assumption 2 in Tan et al., (2019), to detect graph changes, we assume $\Sigma_{jk}(t)$ is right continuous in $[0,1]$.

Assumption 2

Assume $\Sigma_{jk}(t)$ is right continuous and $\Sigma_{jk}^{\pm}(t)$, $\dot{\Sigma}_{jk}^{\pm}(t)$ and $\ddot{\Sigma}_{jk}^{\pm}(t)$ exist for $\forall t\in(0,1)$, $j,k\in[p]$. $\Sigma^{+}_{jk}(0)=\Sigma_{jk}(0)$, $\Sigma^{-}_{jk}(1)=\Sigma_{jk}(1)$ and there is a finite number of discontinuities of $\Sigma(t)$ for $t\in[0,1]$. There exists a constant $M_{\sigma}$ such that

[TABLE]

Details of corollaries and proofs can be found in supplementary notes. Here we only state the theorems. The following Theorem 1 establishes the uniform rates of convergence for the kernel smoothed covariance matrix estimator under the maximum norm. The proof of Theorem 1 is a generalization of Theorem 1 in Tan et al., (2019). Right and left Epanechnikov kernels are used in our paper to detect change points, which makes our rates of convergence different from those in Tan et al., (2019). Corollary 1 in the supplementary notes gives the uniform rates of convergence of $\widehat{\mathbf{\Theta}}^{\pm}(t)$ defined in CLIME.

Theorem 1

Assume that $h\!=\!o(1)$ and that $\log^{2}\!np\cdot\log(np/\sqrt{h})/nh\!=\!o(1)$. For any $0\!<\!a\!<\!b\!<\!1$, under Assumptions 1-2, there exists a positive constant $C$ such that

[TABLE]

for a sufficiently large $n$.

Recall from (9) the definition of $c(1-\alpha,\textbf{E})$. The following Theorem 2 shows that the Gaussian multiplier bootstrap after normalization can be approximately dominated by the quantile of the test statistic $U_{E}$. So it is valid by effectively controlling type I error. Similar to Theorem 2 in Tan et al., (2019), to prove that $U_{E}^{B}$ is a good approximation of $U_{E}$, we define a series of intermediate processes, and we show both $U_{E}^{B}$ and $U_{E}$ can be approximated well by these processes with applications of Theorems A.1 and A.2 in Chernozhukov et al., (2014).

Theorem 2

Assume that $\sqrt{nh^{3}}=o(1)$. In addition, assume that $\operatorname{poly}(s)\cdot\sqrt{\log^{4}(np/\sqrt{h})/nh^{2}}+\operatorname{poly}(s)\cdot\log^{8}(p/h)\cdot\log^{2}(ns)/(nh)=o(1)$, where $\operatorname{poly}(s)$ is a polynomial of $s$. Under the same conditions in Corollary 2 of supplementary notes, we have

[TABLE]

Remark. By integrating the scaling conditions that $\sqrt{nh^{3}}=o(1)$ and $1/\left(nh^{2}\right)=o(1)$, we conclude that the rate of $h$ is $n^{-1/2}\precsim h\precsim n^{-1/3}$, where we say $a_{n}\precsim b_{n}$ when $a_{n},b_{n}\geq 0$ and there exists a constant $C$ such that $C\cdot a_{n}\leq b_{n}$ for any $n=1,2,\ldots$. In practice, we set $h\sim n^{-0.4}$, which satisfies the scaling conditions.

4 SIMULATION

We evaluate the performance of our proposed method through simulations under various scenarios. Each simulation consists of 500 repeated runs with specific parameters including the number of observed time points $N_{T}$, the number of time segments with different precisions and the number of nodes $p$. Then, we uniformly generate $N_{T}$ time points $t\in(0,1)$ along with $\mathbf{\Theta}(t)$. For the estimation of precisions, we use the CLIME estimator at each time point. In all simulations, we set $h\!=\!C_{1}\cdot n^{-0.4}$ and $\lambda\!=\!C_{2}\cdot(h\!+\!\sqrt{log(np/\sqrt{h})/(nh)})$, where $C_{1}$ and $C_{2}$ are the tuning parameters. We uniformly take 50 points in $(0,1)$ as $\mathbf{t}_{B}$ for bootstrapping, and set $\alpha\!=\!0.05$ for hypothesis testing.

4.1 Changing Precision without Sudden Breaks

To test if the method can control type I error rate when there is no sudden change, we simulate time series with smoothly changing precisions. We first generate $N_{U}=100$ sparse integer vectors $\mathbf{U}$ using a multinomial distribution with $p=\{50,100\}$ elements with equal probability, and set the number of trials $m=3$. The graph $\mathbf{\Theta}$ is generated as $\sum_{i\in[{N_{U}}]}\beta_{i}\cdot\mathbf{U}_{i}\mathbf{U}_{i}^{T}$ where $\bm{\beta}$ is the set of coefficients drawn from $\text{Unif}(0,1)^{N_{U}}$. We randomly flip the signs of entries in $\mathbf{\Theta}$ to introduce negative entries. To introduce continuously and smoothly changing pattern of precisions, we generate $\{\mathbf{\Theta}_{1},\mathbf{\Theta}_{2}\}$, and $\mathbf{\Theta}(t)$ is defined as the linear interpolation of these two precision matrices given $t\in(0,1)$. We set $C_{1}=1$ when $p=50$, and $C_{1}=2$ when $p=100$. $C_{2}$ is set to $0.2$ which can be selected by cross-validation for CLIME.

As shown in Table 1, type I error is controlled under $\alpha=0.05$. We also compare our method with covcp introduced by Avanesov et al., (2018) and Dynamic Connectivity Detection (DCD) proposed by Xu and Lindquist, (2015). We use the same parameter settings described in Avanesov et al., (2018); Xu and Lindquist, (2015) for testing, in which the precision estimator with the same parameters used in our simulations is applied for covcp. Different from our approach, both these methods assume a constant $\mathbf{\Theta}$ throughout the time course. Since the covariance structure is changing through the time course, both covcp and DCD rejected the null hypothesis over 90% of all simulations.

4.2 Changing Precision with Sudden Breaks

Within a diagonal matrix $G_{p}=\text{diag}(g_{1},...,g_{p})$, $\mathbf{\Theta}_{ij}$ in $M$ entries are assigned with $a\cdot g_{i}g_{j}$ in which $0<a<1$. The signs are randomly set to introduce negative entries. After the initialization of $\mathbf{\Theta}$, we set $\mathbf{\Theta}_{jj}=\mathbf{\Theta}_{jj}-\min(\Lambda_{min},0)+0.05$ where $\Lambda_{min}$ is the minimum of eigenvalues of $\mathbf{\Theta}$. We generate one matrix with random entries and randomly flip the signs to form $(\mathbf{\Theta}_{1}^{(1)},\mathbf{\Theta}_{1}^{(2)})$, $(\mathbf{\Theta}_{2}^{(1)},\mathbf{\Theta}_{2}^{(2)})$ and $(\mathbf{\Theta}_{3}^{(1)},\mathbf{\Theta}_{3}^{(2)})$. The precision is a linear function of $T$ for each segment, which is an interpolation of each matrix pair. The diagonal precision entries are randomly initialized with $1$ or $9$ with the ratio 9:1. We set $p=50$; $M=50$; $N_{T}=\{1000,1500,2000,2500,3000\}$; $a=0.2$; $C_{1}=2$; $C_{2}=0.4$. A change point is considered successfully identified if it is localized within half of the bandwidth. We compared the power of covcp and DCD with our method for identifying both change points. As shown in Table 2, dGCI has the best power of both change points with relatively small $N$ and is slightly better than DCD with larger $N$. covcp cannot effectively identify both change points even with an increasing sample size when $a=0.2$. All three methods can identify both change points with more than 85% power when $a=\{0.3,0.4\}$.

Additionally, we tested if our approach can locate changed edges besides change points. Table 3 and Table 4 show the increasing ability of the algorithm to distinguish sudden changes in precisions given more samples and larger effect sizes. Here, sensitivity is the ratio of successfully identified edges over changed ones, and the false discovery rate is the ratio of falsely identified edges over all the detected.

5 APPLICATION

5.1 Data Description

A high-resolution fMRI dataset was collected from 15 scanned participants with an audio-video presentation of “Forrest Gump” (Hanke et al.,, 2014, 2016). All the subjects had listened to an audio presentation of the movie, and only one subject had never seen the movie before the experiment (Hanke et al.,, 2016). The original movie was cut into eight parts. Descriptions of semantic conflicts were extended from simple annotations of the cuts and scenes (Hanke and Ibe,, 2016). With the help of scene descriptions, we can locate the scenes of BOLDs. These extended descriptions are references for FC changes. The fMRI data were processed with standard methods and parcellated into 268 nodes using a whole-brain, functional atlas defined in a separate sample (Shen et al.,, 2013). Time series within the same node were averaged, and a sequence of 3539 frames was generated per node, in which TR (repetition time) is 2s. With part of the cerebellum (9 at both hemispheres) and brainstem (1 at the left) missing in some individuals, the total number of ROI is 249. This dataset, along with detailed preprocessing methods, can be attained from “http://studyforrest.org”.

5.2 Objectives

Individuals may share similar reactions when watching the same movie segment. We aim to capture the individual similarity in various contexts of interest. For example, certain regions of the brain may be more sensitive to specific cuts or conflicts. The changed edges identified using our inferential framework could be strengthened, weakened, and reversed in direction. The biological meaning of the identified edges is that the corresponding brain network structure may have a sudden change within a specific time or event, which could be the response of the brain to the external stimuli. Given an fMRI scan from an individual, we perform hypothesis testing, as described in section 2. We are interested to see 1: during the movie, which ROI has the most altered edges, 2: which edge/region is mostly altered across all tests, and 3: which edges are similar among individuals for a specific scene of interest.

5.3 Results

We first scale 3539 consecutive frames into the $(0,1)$ range and we uniformly take 30 points in $(0,1)$ for bootstrapping. Besides these 30 points, 7 time points among 8 cuts are also included for hypothesis testing. We set $C_{1}=1.0$, $C_{2}=0.4$, $\alpha=0.05$, and the number of bootstrap at 500. To remove subject-specific intrinsic BOLDs fluctuations, we calculate the kernel smoothed BOLD sequences with the same $h$ as $\hat{\mu}(t)=\sum_{t^{\prime}=t_{1}}^{t_{T}}k_{h}(t^{\prime},t)\mathbf{X}_{t^{\prime}}/\sum_{t^{\prime}=t_{1}}^{t_{T}}k_{h}(t^{\prime},t)$, and subtract it from the original sequences.

The results from six females and nine males are aggregated as shown in Figure 1a, and the original indices for these individuals are “01”, “02”, “03”, “04”, “05”, “06”, “09”, “10”, “14”, “15”, “16”, “17”, “18”, “19”, and “20”. We highlight some scenes that stand out for the identified edge number. Notably, individual “03”, who had never watched the movie, yields more changed edges compared with the other individuals–especially during the scene at Gump’s House when the main character Jenny returns. The top five scenes with the largest number of identified edge changes for both sexes are labeled as “Doctor’s Office”, “Running”, “Vietnam”, “Bubba’s Grave” and “Jenny’s Grave”. “Forrest Gump” has a considerable number of scenes in which the main character recollects previous experiences, and we only label the primary scene for a given time. These five major scenes are also identified as changed points for all the individuals regardless of sex.

To investigate the functions of the identified changes, we group nodes with definitions of functional and anatomical annotations. The ten regions are prefrontal (PFC, nodes number: 46), motor (Mot, 21), insula (Ins, 7), parietal (Par, 27), temporal (Tem, 39), occipital (Occ, 25), limbic (Lim, 36), cerebellum (Cer, 23), subcortex (Sub, 17) and brainstem (Bsm, 8) and ten canonical brain networks include medial frontal (MF, 29), frontoparietal (FP, 28), default mode (DMN, 18), motor (Mot, 49), visual I (VI, 18), visual II (VII, 8), visual association (VAs, 17), limbic (Lim, 30), basal ganglia (BG, 29) and cerebellum (CBL, 23). (Smith et al.,, 2009; Shen et al.,, 2013; Finn et al.,, 2015). We group results of five interesting scenes from Figure 1a into interactions among nodes in sub-networks. As shown in Figure 1b, limbic, basal ganglia, and cerebellum form a cluster at the important scenes. The BG network, including thalamus and striatum, also has a strong linkage with the motor network and serves as a hub network for the cortical and the subcortical regions.

To make sure the results are not dominated by outliers, we plot the connectome (Figure 1c) with edges identified within at least three individuals at the same scene. The subcortical region, which is the most critical component of the BG network, exhibits reproducible connections with different lobes from the whole brain. It was believed that the basal ganglia structure is the primary modulator or hub for cortical information flow (Lanciego et al.,, 2012). Recent studies have extended the functional scope from motor to cognitive functions (Helie et al.,, 2013). The results indicate that the basal ganglia may play a key role in processing complex stimuli like movies.

It is also intriguing to relate the FC response with published emotional annotations. Labs et al., (2015) annotated the audio-visual stimulus of “Forrest Gump” with three major groups: dimensional emotion attributes (arousal, valence, direction), emotion categories (from admiration to shame) and emotion onset cues (audio, context, face, gesture, narrator, verbal) with detailed descriptions. To investigate if the changes of FC are correlated with any emotions, we use the same kernel in our framework to smooth the collected emotion measure within the bandwidth. We calculate Kendall’s $\tau$ between emotion measures and changed edge numbers for each individual and aggregated sex group. Interestingly, we identified “Love”, which is explained as “Affection for another person”, to be positively correlated with FC changes for almost all the individuals, and it is significant in average connection change for all subjects and female individuals. The contribution of changes with the basal ganglia system may be caused by its dopaminergic role as in neurotransmitter release. Further investigation is needed to relate emotion, emotion response in FC, and emotion-related neurotransmitters.

6 DISCUSSION

Change point detection of dynamic functional connectome involves extensive multiple hypothesis testing. By designing the supreme-maximum test statistic based on the de-biased estimator with asymptotic normality, the hypothesis testing procedure can identify change points along with strong edge changes with appropriate false positive control. This will improve our understanding of the dynamic flow of functional connectivity given complex stimuli. Our method is able to identify psychologically meaningful changes in dynamic functional connectivity, and our results suggest the involvement of the basal ganglia in complex cognitive processes.

We note that our method is not based on the population level modeling, which makes it difficult to aggregate the identified changes. Additionally, the method might lead to inflated type I error without enough observations when the change of connectivity is not smooth, which is not an issue of the validity of the theoretical results but an issue worth studying in applications. Despite the limitations, our algorithm can still detect consistent changes of individuals in real data.

ACKNOWLEDGEMENTS

Here we want to express our gratitude to Dr. Han Liu from the MAGICS (Modern Artificial General Intelligible and Computer Systems) lab at the Northestern University, who provided a valuable chance to connect the authors. We are also grateful for the help of Dr. Kean-Ming Tan in the Department of Statistics at the University of Michigan, whose work inspired and motivated ours. We want to thank Dr. Xilin Shen at Yale Magnetic Resonance Research Center for providing both technical details and insights into imaging processing and analyses.