Inference of Dynamic Graph Changes for Functional Connectome
Dynamic functional connectivity is an effective measure to characterize the brain’s responses to continuous stimulus. We propose an inferential method to detect the dynamic changes of brain networks based on the time-varying graphical model. Whereas most existing methods focus on testing the existence of change points, the fluctuations in the brain network structure represent more signals in many neuroscience studies. We propose a novel method to conduct hypothesis tests on the changes in the dynamic brain networks. We introduce a bootstrap statistic to approximate the supreme of the high-dimensional empirical processes over edges. Our simulations show that this framework can successfully capture change points with changed connectivity. Finally, we apply our method to a brain imaging dataset under the natural audio-video stimulus and illustrate that our approach is able to detect the 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 the framework.
mO#1\ThisStyle\stackengine0.2\SavedStyle#1 \addvbuffer [-0.05ex -1.3ex] OcFTS
The human connectome, which is a collection of brain connectivity, consists of scattered voxels or regions of interests (ROIs) with their entangled spatiotemporal relationship, and functional connectivity (FC) is constructed using BOLD sequences Sporns et al. . Specific statistical measures such as correlation and mutual information are commonly used to characterize the brain network structure Jeong et al. , where a stationary interaction structure is assumed in FC calculation. However, fluctuations of both BOLD sequences and FC have drawn great attention in recent years with the advancing technology and experimental designs Preti et al. . Calhoun et al.  extended the concept of connectome into a dynamic perspective, namely "dynamic connectome", to describe a spectrum of approaches aiming to identify time-varying properties of dynamic functional connectivity (dFC).
A number of dynamic graph estimation frameworks developed under different assumptions have successfully identified interesting dynamic patterns for different brain states Betzel et al. , Gonzalez-Castillo et al. , Wang et al. . A sliding window approach is commonly adopted to link functional brain dynamics and cognitive abilities although the choice of the window length is often subject to debates. Lindquist et al.  proposed multivariate volatility models to refine sliding window methods with exponential weights and used dynamic conditional correlation to measure dFC Preti et al. . Chang and Glover  applied time-frequency analysis to completely circumvent the window limitation, which extends the dynamic connectivity to multiple frequencies.
Furthermore, with the advance of change point detection methods, BOLD sequences can be further divided into segments for a better understanding of network structure within components Xu and Lindquist [2015a]. A general comparison between time segments is achieved by a summary statistic, for example, a likelihood ratio statistic Barnett and Onnela  and a multivariate Gaussian distribution is assumed to facilitate the likelihood analysis. However, it is hard to impose specific hypotheses on the heterogeneity of graphs across the time course because changes could be continuous which invalidates the time segments proposed by likelihood ratio. In Park et al. , a two-sample test based on covariance matrix estimators on different time windows are proposed to detect the change time points of brain connectivity.
Our paper generalizes the previous approaches mentioned above [Xu and Lindquist, 2015a, Barnett and Onnela, 2016, Park et al., 2018] in three aspects. First, we consider a nonparametric time-varying brain network model. Particularly, most existing papers [Barnett and Onnela, 2016, Park et al., 2018] on functional connectivity change points imposed the assumption that the brain network is static along the time trajectory in the null hypotheses. Such assumption violates experimental observations that the brain network is time-varying in cognitive processes [Luscombe et al., 2004, Hagmann et al., 2008, Chen et al., 2016]. In comparison, our model assumes that the brain network varies continuously under both null and alternative hypotheses. Second, we propose a method to detect the local edge change points when the left and right values of precision matrix at some time points are different. Comparing to testing if the precision matrix is globally invariant, our test on the local edge variations can detect more sudden changes in cognitive processes. Third, we show theoretic properties on the validity of our edge number change point tests. In this paper, we propose a combinatoric inferential framework on graph changes, which is adapted from the method introduced by Neykov et al. . Instead of inferring the dynamic graph structure, the bootstrap statistic is redesigned to measure the change or discontinuity of graph given any point in the time sequence. The framework is also proved to effectively control type I error. We then applied this method to the Forrest data set Hanke et al. [2014, 2016] that consists of fMRI data of subjects during the stimuli with an audio-visual movie presentation of "Forrest Gump" and demonstrate the individual differences in brain connectivity change for movie scenes.
In this section, we introduce the time-varying Gaussian graphical model and propose testing procedures for the edge number change point detection.
2.1 Time-Varying Gaussian Graphical Models
Let be a -dimensional random variable and define a time variable . In the time-varying Gaussian graphical model, we assume covariance matrix is defined given and the precision or inverse covariance at is which governs the conditional dependencies for the graph. The inverse covariance matrix is regarded as a representation of the time-varying graph structure. Thus, graphs across time can be retrieved by estimating . A kernel smoothed covariance estimator Yin et al. , Zhou et al.  is used to obtain :
where with as the bandwidth parameter. For simplicity, we use the Epanechenikov kernel described in Lu et al. , i.e., . The precision matrix can then be estimated using the CLIME estimator Cai et al. , i.e., , subject to , where is the -th canonical basis in for each and is a tuning parameter to control the sparsity of . Instead of taking precision as the parameter of interest, Neykov et al.  derived a de-biased estimator for the CLIME estimator
which is proven to be asymptotically normal Neykov et al. . The construction of de-biased precision estimator is able to circumvent asymptotic non-normality of traditional test statistics such as score test Neykov et al.  and we build the test statistic for graphical change inference based on (2).
2.2 Test statistic for graphical change
We will formulate the testing problem and propose the corresponding bootstrapping method.
2.2.1 Graph inference to graph change
Given the precision matrix , we say one edge is changed at the time point if the left and right values of at are different, i.e., . For at , the hypothesis testing of graph changes can be formulated as follows:
Unlike a general graph property defined by Neykov et al.  and Lu et al. , 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 at , we need to consider both sides of , i.e., and , which represent the right and left limits respectively. If there is no change in the current graph, it is expected that .
To perform the statistical test at , first, we generalize (1) into right and left approximated estimates. We define the right-sided covariance estimator as
where to ensure . We also define the left-sided covariance estimator by replacing in (4) with
Following the same procedure above, and are estimated by GLASSO or CLIME using and respectively. Similar to (2), the de-biased estimator of the right and left approximated estimate is
where the notation means we use the left or right estimators respectively in the definition of or . Then, we extend the original test statistic in Lu et al.  to
where and represent the left and right sided values of respectively. We expand the Gaussian multiplier bootstrap for both components as
where , and . With two components defined in (5), the bootstrap statistic is expressed as
in which . The conditional -quantile of given is
In our implementation, we uniformly sample and use sampled to approximate by taking the supreme over instead of . The rejected edge set is then defined as
where . The null hypothesis will be rejected if there exists a such that .
2.2.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 by taking the kernel weighted average of squares:
The normalization term for is defined as . Define , and we apply the normalization term to the test statistic as
and the bootstrapper as
3 Theoretical Results
In this section, we establish the uniform rates of convergence for the kernel smoothed covariance matrix estimator in (4) and the inverse covariance estimator 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 and are allowed to increase. There are details of parameter space, assumptions, corollaries and proofs in the 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. In addition, corollary 1 in supplementary notes gives the uniform rates of convergence of defined in CLIME.
Assume that and that . For any , under some assumptons, there exists a positive constant such that
with probability at least , for sufficiently large .
Assume that . In addition, assume that , where is a polynomial of . Under the same conditions in Corollary 2 of supplementary notes, we have
where the definition of and parameter space can be found in supplementary notes.
We evaluate the efficacy of the framework with simulations under various scenarios. Each simulation consists of 500 repeated runs with specific parameters including number of observed time points , number of time segments with different precisions and number of nodes . We let and uniformly divide into segments. For the estimation of precisions, we use the CLIME estimator at each of the time points. In all simulation, we set and , where and are tuning parameters selected by cross-validation. To generate graphs, we use two methods as listed below
1. With sparse vectors:
We first generate sparse integer vectors using multinomial distribution with elements and size of . We randomly change the signs of entries in to introduce negative entries. The graph is generated as where is the set of coefficients drawn from . The matrices are normalized with unit diagonal.
2. With random entries:
Within a diagonal matrix , in entries are assigned with in which . The signs are randomly set to introduce negative entries. After the initialization of as , we set where is the minimum of eigenvalues of . The matrices are normalized with unit diagonal.
After the specification of graph structures given , , the data is generated by where is sampled from . For bootstrapping and inference, we uniformly sample 50 time points with a fixed interval as reference and testing points. We evaluate the performance under two scenarios with settings and parameters described below.
1. Unchanged random graph:
The single graph for time points is generated with sparse vectors. ; ; ; ; ; ; ; .
2. Changes for blocks in the network:
We generate one matrix with random entries and randomly flip the signs for three times to form , and , and segments are assigned with these three matrices as precisions respectively. The diagonal terms of and are randomly initialized with or with the ratio 9:1. The matrices are not normalized with unit diagonal terms in this case. ; ; ; ; ; .
We use the first and second scenarios to evaluate the type I error and power respectively. The second case serves to test if the framework can correctly detect the change points given enough observations. A change of graph is identified when there exists at least one edge rejected. In this simulation, we aim to test the existence of a point under (3), in which the graph changes.
We calculated type I error rate by counting the number of falsely rejected runs in the simulation with the unchanged random graph. With increasing sample size, the type I error rate decreases within different graphs and type I error rate in all simulations when as shown in Table 1. Thus, our algorithm can successfully control the type I error rate with enough sample size.
|= 200||= 300||= 400||= 500||= 600|
We also evaluated the efficacy of locating change points and changes of precision with changes for blocks in the network. Time points with non-empty rejected sets are regarded as identified change points if they are around the pre-specified change points within half the bandwidth. As shown in Figure 1, after normalization of test statistic, the algorithm successfully detected two change points in 56.8% of all simulations while only 13.8% identified two change points without normalization. We also applied DCD proposed by Xu and Lindquist [2015b] with type I error rate . DCD identified the first and second change points 233 and 237 times but it only detected 10 times of 500 for both of them at the same time. Additionally, we are interested in if our approach is able to locate changed edges besides change points. For each change point, Table 2 shows the ability of the algorithm to pinpoint the changes in the precision given reasonable size of sample and effect.
|change point 1||change point 2|
5.1 Data Description
A high-resolution fMRI dataset was collected from 15 scanned participants with the 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. . The original movie was cut into eight parts. The descriptions of semantic conflicts were extended from simple annotations of cuts and scenes Hanke and Ibe . With the help of scene descriptions, we can locate the scenes of BOLD sequences. These extended descriptions are references for FC changes. FC could change at any time and do not necessarily form clusters under this sophisticated setting.
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. . Time series within the same node were averaged and a sequence of 3539 frames were generated per node, in which TR is 2s. With part of cerebellum (9 at both hemispheres) and brainstem (1 at the left) missing in some individuals, the total number of regions of interest is 249. Data with details of preprocessing and more information can be obtained from "http://studyforrest.org".
Individuals may share similar reactions when they are 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. We are especially interested in studying the covert brain reactions in these subsets of individuals and brain regions of interests. In this study, at least one change in FC of reactions to the movie progression is expected, and the number of change points is assumed finite. These FC changes are located using the inferential framework. Given an fMRI scan from individual , we perform hypothesis testing as (3) and we are interested to see 1: during the movie, which part has the most altered edges; 2: which edge/region is mostly altered across all tests; 3: which edges are similar among individuals for a specific scene of interest.
We first scale 3539 consecutive frames into the range and we uniformly take 30 points in for our method. Besides these 30 points, 7 time points between 8 cuts are also included for hypothesis testing. We set , , , and number of bootstrap as 500.
The results from six females and nine males are aggregated as shown in Figure 2a. We highlight some scenes that stand out for the identified edge number. Interestingly, individual 03, who had never watched the movie, yields more changed edges compared with other individuals especially at the scene at Gump’s House with main character Jenny’s return. The top five scenes with the most number of identified edges 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 scene movements with the process of the main character recollecting the experiences so we only label the primary scene for a given time. These five major scenes are also identified as changed points for all the individuals in each sex.
To further investigate the function of identified changes, we use the definition of subnetworks and lobes with functional or anatomical annotations to group nodes. The ten regions are prefrontal(nodes number: 46), motor(21), insula(7), parietal(27), temporal(39), occipital(25), limbic(36), cerebellum(23), subcortex(17) and brainstem(8) and ten canonical brain networks include medial frontal(29), frontoparietal(28), default mode(18), motor(49), visual I(18), visual II(8), visual association(17), limbic(30), basal ganglia(29) and cerebellum(23). Smith et al. , Shen et al. , Finn et al. . We take the results of five interesting scenes from Figure 2a and group them into interactions among nodes in different sub-networks. As shown in Figure 2b, limbic, basal ganglia(BG) and cerebellum form a cluster at the important scenes. Basal ganglia network including thalamus and striatum also has a strong linkage with motor network and serves as a major hub network for cortical and subcortical regions.
To make sure the results are not dominated by an outlier with many interactions of basal ganglia network, we also plot the connectome (Figure 2c) with edges identified within at least three individuals at the same scene. The subcortical region, which is the most critical component in basal ganglia network definition, exhibits reproducible connections between different lobes from the whole brain. It was believed that the functional organization of basal ganglia is a loop, in which basal ganglia is the modulator or hub for cortical information flow Lanciego et al. . The previously targeted role of BG used to be mainly motor function, but recent studies have extended the scope to cognitive functions Helie et al. . The results indicate that basal ganglia may play a key role in processing complex stimuli such as movies.
|Change-Emotion Correlation||Emotion: "Love"|
It is also intriguing to relate the FC response with published emotional annotations Labs et al. . Hanke et al 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 Labs et al. . To investigate if changes of functional connectivity are correlated with any emotions, we first use the same kernel used in our framework to smooth the collected emotion measure within the bandwidth, and we calculate the Kendall’s between emotion measures and changed edge numbers for each individual and aggregated gender groups. Interestingly, we identified "Love", as explained as "Affection for another person", to be positively correlated with FC changes for almost all the individuals. As an emotion-related region to dopamine, the contribution of changes with 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 and Conclusion
Change point detection of dynamic functional connectome involves extensively multiple testing. By designing the supreme-maximum test statistic based on the de-biased estimator with asymptotic normality, the hypothesis testing procedure can identify both change points and strong changes at change points, which also controls the type I error rate. This will help to understand the details of the dynamic flow of functional connectivity given complex stimuli. The consistency of identified changes among individuals also shows the ability of this algorithm to identify psychologically meaningful changes in dynamic functional connectivity and the results suggest the potential of basal ganglia in complex cognitive processes.
However, this method is not based on the population level, which makes it difficult to aggregate the information of identified changes. Besides, with an extremely high frequency of observations, it exhibits the best performance given the assumption of finite change points with the discontinuity of connectome though the change can be continuous. Despite the simplified assumption, with normal sample size, the algorithm is still powerful to detect consistent changes for individuals in the population.
- Sporns et al.  O. Sporns, G. Tononi, and R. Kotter. The human connectome: A structural description of the human brain. PLoS Comput Biol, 1(4):e42, 2005.
- Jeong et al.  J. Jeong, J. C. Gore, and B. S. Peterson. Mutual information analysis of the eeg in patients with alzheimer’s disease. Clin Neurophysiol, 112(5):827–35, 2001.
- Preti et al.  M. G. Preti, T. A. Bolton, and D. Van De Ville. The dynamic functional connectome: State-of-the-art and perspectives. Neuroimage, 160:41–54, 2017.
- Calhoun et al.  V. D. Calhoun, R. Miller, G. Pearlson, and T. Adali. The chronnectome: time-varying connectivity networks as the next frontier in fmri data discovery. Neuron, 84(2):262–74, 2014.
- Betzel et al.  R. F. Betzel, M. Fukushima, Y. He, X. N. Zuo, and O. Sporns. Dynamic fluctuations coincide with periods of high and low modularity in resting-state functional brain networks. Neuroimage, 127:287–297, 2016.
- Gonzalez-Castillo et al.  J. Gonzalez-Castillo, C. W. Hoy, D. A. Handwerker, M. E. Robinson, L. C. Buchanan, Z. S. Saad, and P. A. Bandettini. Tracking ongoing cognition in individuals using brief, whole-brain functional connectivity patterns. Proceedings of the National Academy of Sciences of the United States of America, 112(28):8762–8767, 2015.
- Wang et al.  J. X. Wang, J. Bartolotti, L. A. Amaral, and J. R. Booth. Changes in task-related functional connectivity across multiple spatial scales are related to reading performance. PLoS One, 8(3):e59204, 2013.
- Lindquist et al.  M. A. Lindquist, Y. Xu, M. B. Nebel, and B. S. Caffo. Evaluating dynamic bivariate correlations in resting-state fmri: a comparison study and a new approach. Neuroimage, 101:531–46, 2014.
- Chang and Glover  C. Chang and G. H. Glover. Time-frequency dynamics of resting-state brain connectivity measured with fmri. Neuroimage, 50(1):81–98, 2010.
- Xu and Lindquist [2015a] Y. Xu and M. A. Lindquist. Dynamic connectivity detection: an algorithm for determining functional connectivity change points in fmri data. Front Neurosci, 9:285, 2015a.
- Barnett and Onnela  I. Barnett and J. P. Onnela. Change point detection in correlation networks. Scientific Reports, 6, 2016.
- Park et al.  Hae-Jeong Park, Karl J Friston, Chongwon Pae, Bumhee Park, and Adeel Razi. Dynamic effective connectivity in resting state fmri. Neuroimage, 180:594–608, 2018.
- Luscombe et al.  Nicholas M Luscombe, M Madan Babu, Haiyuan Yu, Michael Snyder, Sarah A Teichmann, and Mark Gerstein. Genomic analysis of regulatory network dynamics reveals large topological changes. Nature, 431(7006):308–312, 2004.
- Hagmann et al.  Patric Hagmann, Leila Cammoun, Xavier Gigandet, Reto Meuli, Christopher J Honey, Van J Wedeen, and Olaf Sporns. Mapping the structural core of human cerebral cortex. PLoS Biol, 6(7):e159, 2008.
- Chen et al.  Janice Chen, Yuan Chang Leong, Kenneth A Norman, and Uri Hasson. Shared experience, shared memory: a common structure for brain activity during naturalistic recall. bioRxiv, 2016.
- Neykov et al.  M. Neykov, J. W. Lu, and H. Liu. Combinatorial inference for graphical models. Annals of Statistics, 47(2):795–827, 2019.
- Hanke et al.  Michael Hanke, Florian J Baumgartner, Pierre Ibe, Falko R Kaule, Stefan Pollmann, Oliver Speck, Wolf Zinke, and Jörg Stadler. A high-resolution 7-tesla fmri dataset from complex natural stimulation with an audio movie. Scientific data, 1:140003, 2014.
- Hanke et al.  Michael Hanke, Nico Adelhöfer, Daniel Kottke, Vittorio Iacovella, Ayan Sengupta, Falko R Kaule, Roland Nigbur, Alexander Q Waite, Florian Baumgartner, and Jörg Stadler. A studyforrest extension, simultaneous fmri and eye gaze recordings during prolonged natural stimulation. Scientific data, 3:160092, 2016.
- Yin et al.  J. Yin, Z. Geng, R. Li, and H. Wang. Nonparametric covariance model. Stat Sin, 20:469–479, 2010.
- Zhou et al.  S. H. Zhou, J. Lafferty, and L. Wasserman. Time varying undirected graphs. Machine Learning, 80(2-3):295–319, 2010.
- Lu et al.  J. W. Lu, M. Kolar, and H. Liu. Post-regularization inference for time-varying nonparanormal graphical models. Journal of Machine Learning Research, 18, 2018.
- Cai et al.  T. Cai, W. D. Liu, and X. Luo. A constrained l(1) minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106(494):594–607, 2011.
- Neykov et al.  M. Neykov, Y. Ning, J. S. Liu, and H. Liu. A unified theory of confidence regions and testing for high-dimensional estimating equations. Statistical Science, 33(3):427–443, 2018.
- Xu and Lindquist [2015b] Yuting Xu and Martin A Lindquist. Dynamic connectivity detection: an algorithm for determining functional connectivity change points in fmri data. Frontiers in neuroscience, 9:285, 2015b.
- Hanke and Ibe  Michael Hanke and Pierre Ibe. Lies, irony, and contradiction?an annotation of semantic conflict in the movie" forrest gump". F1000Research, 5, 2016.
- Shen et al.  Xilin Shen, Fuyuze Tokoglu, Xenios Papademetris, and R Todd Constable. Groupwise whole-brain parcellation from resting-state fmri data for network node identification. Neuroimage, 82:403–415, 2013.
- Smith et al.  Stephen M Smith, Peter T Fox, Karla L Miller, David C Glahn, P Mickle Fox, Clare E Mackay, Nicola Filippini, Kate E Watkins, Roberto Toro, Angela R Laird, et al. Correspondence of the brain’s functional architecture during activation and rest. Proceedings of the National Academy of Sciences, 106(31):13040–13045, 2009.
- Finn et al.  E. S. Finn, X. Shen, D. Scheinost, M. D. Rosenberg, J. Huang, M. M. Chun, X. Papademetris, and R. T. Constable. Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nat Neurosci, 18(11):1664–71, 2015.
- Lanciego et al.  José L Lanciego, Natasha Luquin, and José A Obeso. Functional neuroanatomy of the basal ganglia. Cold Spring Harbor perspectives in medicine, 2(12):a009621, 2012.
- Helie et al.  Sebastien Helie, Srinivasa Chakravarthy, and Ahmed A Moustafa. Exploring the cognitive and motor functions of the basal ganglia: an integrative review of computational cognitive neuroscience models. Frontiers in computational neuroscience, 7:174, 2013.
- Labs et al.  Annika Labs, Theresa Reich, Helene Schulenburg, Manuel Boennen, Gehrke Mareike, Madleen Golz, Benita Hartigs, Nico Hoffmann, Sebastian Keil, Malú Perlow, et al. Portrayed emotions in the movie" forrest gump". F1000Research, 4, 2015.