1 Introduction
During their lifetime, 60.7% of men and 51.2% of women experience at least one potentially traumatic event (Kessler et al., 1995). Of those experiencing traumatic events, 1040% develop psychiatric symptoms of clinical relevance (Breslau, 2009) such as posttraumatic stress disorder (PTSD). PTSD is one of the most common mental disorder in USA and results in significant impairments of psychological and physical health (Kessler et al., 1995). There is therefore great interest in identifying neural mechanisms that contribute to some individuals’ resilience in the face of stressors such as trauma exposure, especially given that many individuals do not maintain impairing high levels of psychiatric symptoms.
Previous studies of risk for psychiatric disorders such as PTSD after trauma have identified functional abnormalities in several brain areas including amygdala, hippocampus, ventromedial prefrontal cortex (mPFC), anterior cingulate cortex (ACC) and insula (Brown et al., 2014), that are associated with emotion, memory and executive functions (Shalev et al., 2017). In individuals exposed to high levels of childhood and adulthood trauma, activation of the vmPFC (Stevens et al., 2016) and hippocampus (van Rooij et al., 2016) are associated with trait psychological resilience, and decreased risk for PTSD. However, it is well recognized that these areas do not act in isolation, and that the diseases severity can often be better accounted for by considering networks of interacting brain regions. Such interactions or coactivations can be captured via functional connectivity (FC) that encodes the temporal coherence between regions (Smith et al., 2011). Resilient individuals who experienced trauma but did not have chronic PTSD show lower restingstate functional connectivity (rsFC) between amygdala and insula (Rabinak et al., 2011), greater amygdalahippocampus rsFC (Sripada et al., 2012), and lower amygdalaACC/dorsal mPFC rsFC (Brown et al., 2014), relative to age and traumaseveritymatched groups with current PTSD. Lower rsFC in salience network but increased rsFC in default mode network were also observed in Koch et al. (2016).
The majority of connectivity studies of risk and resilience to traumarelated psychiatric symptoms so far have utilized static functional connectivity methods, which assume that functional connectivity between regions is stationary over time (Koch et al., 2016). However, static FC does not account for changes in the network over the scanning period that may arise due to nonneuronal related factors, as well as variations in the BOLD signal mean and variance over time (Hutchinson et al., 2013; Lindquist et al., 2014). Further, dynamic variations in rsFC are relevant to changes in vigilance (Thompson et al., 2013), arousal (Chang et al., 2013), emotional state (Cribben et al., 2012), and behavioral performance (Jia et al., 2014). Compared to stable traitlevel functional network alterations that may be apparent only when estimating static group level networks, one expects individual level networks to exhibit dynamically changing states reflective of psychological changes. A recent study showed that dynamic rsFC models allowed for higher accuracy than static rsFC in differentiating individuals with PTSD from traumaexposed controls (Jin et al., 2017), suggesting that the network alterations involved in risk and resilience to PTSD after trauma exposure are not static, but may also involve dynamic alterations in the networks that are engaged during rest.
Hence, there is a serious and timely need for statistical approaches for estimating dynamic brain functional networks, but there are several unmet challenges. Perhaps, the most commonly used strategy has been a sliding window approach (Allen et al., 2014; Handwerker et al., 2012). While the sliding window technique is a valuable tool for investigating temporal dynamics of functional brain networks, there are some known limitations associated with this approach (Lindquist et al., 2014) such as the choice of window length that is often difficult to determine. Alternate approaches involve hidden Markov models which are capable of detecting reoccurring quasistable connectivity states (Baker et al., 2014). These approaches assume that the brain network reorganizes itself into any one of a fixed number of states at each time point, but they may not ensure that the networks across consecutive time points are more similar, which is a natural assumption.
Another popular class of approaches is based on change point models that prioritize the network homogeneity across consecutive time points (Hutchinson et al., 2013; Cribben et al., 2012, 2013; Xu and Lindquist, 2015), with the network changes occurring as discrete jumps in the time series. Change point models visualize the dynamic brain network as a collection of state phases representing various modulations in the network, where the network remains constant across all time points within a state phase (defined by the period between consecutive change points), but varies between state phases. The discrete network changes are meant to serve as an approximation of a more realistic scenario of slower changes in the network consistent with the timescale of haemodynamic brain function. One potential interpretation of change points is that they essentially summarize transition periods spanning multiple consecutive time points via a representative discrete jump, which is often meant to speed up computations and ensure scalability for high dimensional fMRI applications. Discrete changes also facilitate greater interpretability and easier visualizations.
Although existing change point models have been somewhat successful in describing the temporal changes in the brain network, there are some unaddressed challenges which involve the detection of rapidly evolving brain organization. This is especially significant, given that brain networks can change within as little as 3060s (Shirer et al., 2012). Such rapid fluctuations will divide the scanning period into a collection of narrow state phases, each containing only a few time scans from which to estimate the brain network, thus creating difficulties in terms of estimating the change points (Kundu et al., 2018). A promising solution for estimating rapidly evolving network changes is to incorporate brain anatomical information in terms of structural connectivity or SC (computed via DTI data) to guide the estimation of dynamic brain functional networks. The hope is that incorporating supplemental SC strengths may help increase the power and accuracy of detecting dynamic network changes. There is some preliminary work on fiber centered functional connectivity estimation (Li et al., 2013), but these articles only focus on pairs of brain regions that are structurally connected, whereas our goal is to investigate the whole brain connectome.
Of course, the motivation for structurally guided dynamic functional connectivity comes from a wellestablished literature illustrating the relationship between static FC and SC. In particular, there is strong evidence that white matter fiber tracts regulate static FC (Damoiseaux and Greicius, 2009; Sporns, 2013). Based on such evidence, there has been some limited development of static FC approaches guided by SC knowledge (Xue et al., 2015; Hinne et al., 2014; Messe et. al., 2014). However, the vast majority of the existing approaches, including the above ones, use multisubject data which requires registration of images to a shared template under the assumption that the volumes are similar and can be matched, thus ignoring the variability in cortical anatomy and function (Zhu et al., 2012). Moreover such group level networks do not account for variations on the network across individuals that arise naturally in PTSD studies due to the underlying heterogeneity. Approaches using single subject data were recently proposed by Ng et al. (2012) and PinedaPardo et al. (2014). Their methods are based on adaptive graphical lasso (Wang, 2012) that specify the edgespecific shrinkage parameters as deterministic functions of the SC information. Although useful, such a deterministic specification may not adequately capture the complex underlying structurefunction relationships, or account for heterogeneity in FC for a given SC strength, and it may be more susceptible to misspecification of prior SC knowledge. To overcome these and other associated difficulties, Higgins et al. (2018) proposed a more flexible approach that models the edgespecific shrinkage parameters as a random function of SC, which can accommodate more diverse FCSC relationships. This approach was shown to have good accuracy and reproducibility in estimating the static brain network guided by brain SC information.
It is nontrivial to generalize the above approaches for static network estimation to our settings of interest involving the estimation of dynamic brain networks guided by SC information. Clearly, such an approach is appealing in terms of fusing multiple imaging modalities (i.e. fMRI and DTI) in order to estimate temporal FC changes more accurately. We propose a novel approach to solve this problem which relies on a Bayesian change point model that partitions the time course into a collection of nonoverlapping state phases characterized by distinct brain networks, guided by SC information. Both the location and the number of change points are not known in advance and estimated in an unsupervised manner. The SC information modulates the presence/absence of functional connections and is incorporated when modeling the unknown edgespecific shrinkage parameters for the inverse covariance matrix elements in a Gaussian graphical model (GGM), using a similar strategy in Higgins et al. (2018). However unlike that paper, our focus is on dynamic connectivity instead of static connectivity. Our change point estimation strategy leverages the dynamic connectivity regression (DCR) method proposed in Cribben et al. (2013) but is distinct in using the SC information to guide the estimation of the dynamic brain network, and in developing a novel subnetwork sampling scheme that is able to detect transition periods instead of discrete jumps. The subnetwork sampling strategy is appealing in terms of making the approach scalable to a large number of nodes (a valuable feature in whole brain connectivity studies) and being able to report approximate measures of uncertainty in change point estimation.
We develop a computationally efficient algorithm to obtain a maximumapriori (MAP) estimate of dynamic network. Extensive simulation studies illustrate that the proposed method is better able to detect true changes in the network and has a higher network estimation accuracy compared to existing dynamic connectivity approaches that do not incorporate SC knowledge, and is also superior to alternate network modeling approaches that assume the network is stationary. We apply our approach to estimate the dynamic brain network using resting state fMRI data, and DTI measurements for a sample of civilian participants who have experienced high levels of lifetime trauma and have high rates of related mood and anxiety disorders, drawn from the Grady Trauma Project. The proposed approach is able to detect important dynamic connections that are common across participants in regions of the brain previously implicated in PTSD risk and resilience. Higher temporal fluctuations for connections lying within these regions are shown to result in a decrease in resilience. Moreover, the dynamic functional network computed under our approach produces significantly better prediction for resilience (Ioannidis et al., 2016) in trauma exposed individuals compared to other existing approaches, which highlight its potential as a promising neuroimaging biomarker in trauma exposure studies, thus adding to the findings in Jin et al. (2017).
2 Methods
Our goal is to develop a novel change point approach based on Gaussian graphical models (GGM) for estimating dynamic connectivity while incorporating SC knowledge. GGMs assume that the fMRI measurements are normally distributed and are characterized by a sparse inverse covariance or precision matrix that has zero offdiagonals corresponding to absent edges in the network. Moreover, the nonzero elements of the precision matrix encode the strength of the important edges. The change point model specifies precision matrices that remain constant within a state phase, but changes between state phases, thereby reflecting dynamic connectivity. We note that in a dynamic GGM, the pattern of zeros in the precision matrix at each time point essentially provides all the necessary information about the timevarying network.
To fix notations, denote as the vector of spatially distributed fMRI measurements over voxels or regions of interest (ROI), at the th time point (
). Denote the SC probability corresponding to the edge
as where and denote the corresponding SC probability matrix as . These SC probabilities are obtained using probabilistic tractography based on DTI data, and are made symmetric (i.e. ) as in Higgins et al. (2018). We specify the following dynamic Gaussian graphical model (GGM)(1) 
where
refers to a multivariate Gaussian distribution with mean
and covariance matrix , denotes the inverse covariance or precision matrix at time point that depends on the timevarying network characterized by the vertex set and edge set , as well as brain SC information . The vertex set consists of a set of predefined voxels/ROIs or nodes, the edge set contains the set of all edges present in , and encodes the strength of the time varying FC. Our goal is to develop an algorithm that is able to learn the best fitting partition of the time course defined by the boundary points so that remains constant for consecutive time points except for discrete jumps at the change points. These are unknown for our problems of interest, and estimated in an unsupervised and dataadaptive manner. For conciseness, we denote as the constant precision matrix for the th state phase corresponding to the interval . throughout the article.2.1 Structurally Informed Precision Matrix Estimation
Suppose the time course is partitioned into prespecified nonoverlapping intervals , where contains time scans such that . Conditional on these intervals, the sample mean and the covariance matrix for the partition are given by and respectively. The precision matrix estimate for the th state phase is obtained as a MAP estimator under a Bayesian version of the GGM in (1) involving appropriate prior distributions as follows
(2) 
where represents the prior distribution, is the intractable normalizing constant for (Higgins et al., 2018), denotes vector of edgespecific shrinkage parameters in , is the overall penalty parameter (higher values imply greater network sparsity and viceversa),
denotes the exponential distribution, double exponential distribution, and lognormal distribution respectively,
is the indicator function, and denotes the collection of all symmetric and positive definite matrices. We assign a lognormal type prior on , which restricts the shrinkage parameters to nonnegative values and enables us to model the edgespecific shrinkage parameters in terms of the SC strengths. This prior includes unknown hyperparameters which represent edgespecific baseline effects that are independent of the given SC knowledge, andthat are positive random variables controlling the average effect of SC on FC across the different statephases. These hyperparameters are unknown and are modeled using priors
and respectively under a fully Bayesian specification ().The anatomically informed prior on the shrinkage parameters in (2) specifies a probabilistic relationship between the edge specific shrinkage parameters and the given SC knowledge via . In particular, increasing positive values of implies an increasing dependence of FC on the given SC, while small values of implies a negligible relationship. Moreover, the shrinkage parameters are stochastically monotonically decreasing with respect to the SC strength, under the restriction . This implies that as the SC strength for the edge is increased, the corresponding shrinkage parameter will take smaller values in probability, resulting in values which are away from zero. In other words, the presence (absence) of FC is encouraged for large (small) values of the corresponding SC strength, via the shrinkage parameters . The above model specifications are designed to respect the relationship between FC and SC commonly observed in literature (Higgins et al., 2018). Additionally, the baseline effects () corresponds to variations in neuronal activity that are independent of the brain anatomy. Overall, increasing (decreasing) absolute values of discourages (encourages) the presence of the corresponding edge, independent of the anatomical information. The proposed model enables (a) more flexibility in the FCSC relationship by allowing the possibility of strong FC corresponding to poor SC strengths, and viceversa; and (b) heterogeneity in FC across edges which possesses similar SC strength. We note that (2) adapts the approach in Higgins et al. (2018) to the case of dynamic FC involving state phase specific networks.
The posterior distribution for parameters corresponding to the state phase is given by
(3) 
where denotes the conditional distribution of given . One can optimize the logposterior distribution (obtained by taking the logtransform of the above posterior) to calculate the MAP estimate for the parameters as , where
(4)  
All the parameters in the objective function are updated iteratively until convergence for a fixed value of the sparsity parameter . The precision matrix is updated given other parameters using the existing graphical lasso algorithm (Friedman et al., 2008), whereas are updated via a closed form expressions and are updated via a NewtonRaphson step since a closed form solution does not exist. The objective function (4) is optimized over a range of penalty parameter values and we choose the value of the penalty parameter than optimizes some goodness of fit score such as the Bayesian Information Criteria (BIC) (Yuan and Lin, 2007) that guard against overfitting. For a given state phase, the update steps are described in Appendix A.
2.2 Structurally Informed Change point estimation
The above dynamic precision matrix estimation was conditional on prespecified change points that are unknown in practice. We now describe a dataadaptive algorithm for change point estimation, which is motivated by the DCR approach in Cribben et. al (2012), but has important differences. The approach involves a greedy partitioning scheme which begins by obtaining estimates for the precision matrix based on the entire time series obtained by minimizing equation (4) subject to no change points. The model is fit over a range of values and the minimum BIC score is recorded.
Upon completion of this step, the time course is split into two partitions and , with the understanding that any split of the time course that results in a BIC reduction is acceptable. For each partition, equation (4) is refit again so as to obtain a series of graphs corresponding to a path of values, and the optimal network is selected using BIC. This procedure is repeated along the entire time path, with the data partitioned into two subsets corresponding to split points ranging from to . In order to ensure reliable estimation of the network in each state phase we fix the minimum number of time scans between consecutive change points. The partition with the smallest combined BIC score is chosen and, if its value is less than the BIC score for the entire data set, the corresponding split point is identified as the first change point. The procedure continues by recursively applying the same method to each individual partition element until they can no longer be split any further, i.e. no additional reduction in BIC is seen. As the final output, the DCR algorithm will have split the entire time course into nonoverlapping partitions . The number and location of the change points are thus determined in a data adaptive manner and guided by SC knowledge, since equation (4) incorporates the given SC strengths. We denote the above change point estimation approach as the structurally informed connectivity change point detection (siCCPD) approach.
2.3 A subnetwork sampling scheme
In our experience, the change point estimation approach used described above (motivated by Cribben et al. (2012)) may not yield accurate change point detection results under rapid transitions of the network, and unfortunately it is not be scalable to a large number of nodes, as illustrated via extensive simulations in this article. To overcome such difficulties, we propose a novel heuristic subnetwork sampling strategy, which is based on the key observation that alterations in the network represented by change points may be detected using only a subset of the nodes in the network, as long as the edges connected to at least one of these subset of nodes undergoes temporal connectivity changes. One potential pitfall of this strategy is that the changes corresponding to the subset of nodes may not be strong enough compared to the overall changes in the network at a given time point, in which case it may not be detected. However, by repeatedly sampling random subset of nodes and computing the connectivity change points under each of these subsets, the hope is that the true underlying change points will be detected by a large proportion of these subnetworks. By drawing an adequate number of subnetworks and computing change points for each subnetwork, one can ensure that all important connectivity changes between subsets of nodes in the network are accounted for.
Our strategy is to apply the proposed approach in Section 2.2 to a randomly sampled subset of nodes , which yields a set of change points associated with the corresponding dynamic subnetwork. We then repeat this process over multiple randomly sampled subnetworks, thereby generating a collection of sets of change points. Then, the frequency of each time point being identified as a change point over the collection of sampled subnetworks is computed. Finally, a systematic thresholding approach (described below) is proposed to detect the important change points as those which show up most frequently across the subnetworks. We note that by increasing the number of subnetwork samples, the accuracy for change point detection is expected to increase although it comes at a cost of increased computation time. However, the approach can be parallelized over subnetworks, thereby resulting in computational speedups, when needed.
Suppose we sample subnetworks containing nodes each, where ROIs are selected randomly from for each subnetwork (the number of nodes may also be made different across subnetworks in principle). We then run the proposed siCCPD method for each subnetwork to obtain a set of change points for the th subnetwork as (), with the understanding that the number and location of change points may vary across subnetworks. The set of estimated change points across all the subnetworks is aggregated to obtain the set of all identified change points . Moreover, we also calculate the corresponding frequencies with which each time point was identified as a change point across the subnetworks. This is denoted by , with for those time points that are not identified as change points under any subnetworks. The frequencies can be interpreted as an approximate probability for each time point to be a change point, and any time point for which is treated as a potential change point. Given and , we then apply a grouping approach to find out representative change points within this set, by identifying clusters of change points. In particular, two estimated change points are grouped into one cluster if they are consecutive time points or spaced one time point apart. This grouping method results in distinct groups or clusters of change points (see Figure 1 for a visualization pertaining to our simulation study).
Each cluster of change points is representative of a transition period for the network, that is more consistent with the slower timescale of the haemodynamic activity. In order to eliminate false clusters of change points, we adopt the following thresholding mechanism. First, the overall frequency of a particular cluster of change points is taken as the sum of the frequencies for all the change points within that cluster. Subsequently, all clusters having a combined frequency below a certain threshold are identified as false positives and eliminated. Based on extensive empirical experiments, we propose a threshold of for the combined frequency for clusters, which proved to be a good choice in terms of identifying true change points and eliminating false positives. Alternatively, the number of clusters may also be determined in a data adaptive manner using some goodness of fit measure such as BIC. Once we get the final groups, one can designate a representative single change point for each cluster as the median time point for that cluster. Thus, one is able to obtain both discrete jumps in FC that is the hallmark of existing change point models, as well as estimate transition periods, which is an added novelty of the approach.
3 Numerical Studies
Data Generation
We conducted extensive simulation studies to evaluate the performance of the proposed approach, under different network structures. In the first set of simulations (Scenario I), we generated data from an underlying change point model with three change points (i.e. four state phases). A network and the corresponding precision matrix was constructed at each time point (described below), and these were constant within each state phase. The measurements were generated under a Gaussian distribution characterized by the timevarying precision matrix. We generated data for regions and with time points. In a second set of simulations (Scenario II), we allow the network to change more slowly over time, that is more consistent with the timescale of the haemodynamic activity. In particular, instead of three change points as in the first scenario, we now have three transition periods, each comprising seven consecutive time points. A certain percentage of the edges are flipped from the network at the previous time point to obtain the modified network for the next time point within each transition period. The network is assumed to be constant between two consecutive transition periods. This scenario is more challenging since the network changes multiple times over the course of the experiment. The goal of this experiment is to investigate if the proposed approach can detect the transition periods and whether it can approximate the true dynamic network sufficiently well when the underlying assumptions of the proposed model may not hold.
We generated the functional network of the first bin using three different network structures: (a) ErdosRenyi random graph (Erdos and Renyi, 1960) that has the same probability for all connections; (b) scalefree graph that uses the preferential attachment model of Barabasi and Albert (1999); and (c) smallworld random graph that was obtained using the Watts and Strogatz (1998) model, and which is motivated by the characteristics of real life brain networks derived from fMRI data. Two different average densities for the network were used, 0.15 and 0.3, which indicate sparse connectivity patterns in brain organization (Eavani et al., 2015). Once the functional network for the first state phase was generated, the networks of the other three bins are generated by flipping a certain proportion of the edges of the first bin independently. This implies a fixed proportion of edges in the functional network for the first bin were changed to be absent in the second functional network, whereas a fixed proportion of pairs of nodes that were not connected in the first bin became connected in the second bin. This procedure is repeated independently for all the bins to generate distinct networks for different state phases that have similar densities and also share common patterns, while also exhibiting inherent differences.
Once the networks are generated, the corresponding precision matrix at each time point was constructed by randomly generating offdiagonals from a uniform(1,1) distribution corresponding to the edges in the network at that time point, while fixing those offdiagonals corresponding to absent edges to be zero. The diagonal elements for the matrix are then recalibrated by summing the absolute values of the offdiagonals in the corresponding row/column, and adding a positive constant to this sum, so as to obtain a diagonally dominant matrix that is positive definite. Under a GGM, the fMRI observation at a particular time point was generated from a Gaussian distribution having the timedependent precision matrix constructed above. Moreover, we generated the SC information based on the true FC, as in Higgins et al. (2018). We specify that most edges that had strong FC () in first bin also have a strong SC (), while edges with weak SC (0.5) could have either strong or weak FC. We note that since the networks in the latter bins were obtained by flipping edges in the network corresponding to the first bin, it is very likely to have a nontrivial number of edges with small SC but large FC and large SC with weak FC, which violates the modeling assumptions and presents a challenging scenario.
In addition to Scenarios I and II, we also reported results under ErdosRenyi networks for 100 nodes with three true change points, to test the performance in higher dimensions, and investigated the scenario involving a large number of change points (10) with 50 nodes. These challenging settings help us evaluate the performance of the methods for high dimensions and large number of fluctuations in the dynamic network. For each simulation setting, 25 replicates were used.
Comparison and Assessment
We evaluate the performance of the proposed method with respect to two aspects: the ability to estimate true change points and accurately estimate the network. For the first scenario involving discrete jumps in connectivity, the estimated change point is said to be a true positive when the temporal distance between the true and estimated change point is less or equal than 2. For the second scenario involving transition periods of seven consecutive change points, we denote an estimated change point to be true positive if it lies within the from the center of the transition period. We report the proportion of true change points detected, and the number of false change points detected, averaged over all replicates. The performance of the graph estimation was assessed by comparing the estimated network with the true network at each time point, using the area under the curve (AUC) for the receiver operating characteristic (ROC) curve over different sparsity levels of the network obtained by varying the tuning parameter in (4). We also look at sensitivity and specificity for the optimal network obtained by minimizing the BIC. Sensitivity measures the power to detect true connections and is equal to the proportion of true edges that are successfully detected in the estimated network. Specificity represents the proportion of the absent edges that were successfully detected as absent, and is an indicator of how well the false discoveries are controlled for. Except for false change points detected, higher values of the other metrics imply a more accurate network estimation. All metrics are reported after averaging over all time points.
We compare our method with two competing approaches: dynamic connectivity regression (DCR) for single participants proposed by Cribben et al.(2013) and hierarchical Bayesian structurally informed Gaussian graphical model (siGGM) by Higgins et al. (2018). DCR estimates the dynamic brain network without incorporating SC knowledge, whereas the siGGM approach estimates the static network over the experimental time course while accounting for the SC information. We used default setup in the Cribben et al. (2012) paper to implement the DCR method, i.e. minimum block size of 30, significance level for permutation test of 0.05 and the number of bootstrap is set to be 1000. In contrast the minimum block size under the proposed approach () was set to 510 for greater flexibility, and the subnetwork size was assumed to be 10. While both the change point detection and the network estimation performance were reported under the DCR approach, the siGGM only reports the network estimation performance, since it does not accounting for dynamic changes. We note that the DCR approach was not feasible for the 100 node example due to an unrealistic computational burden.
Results
Table 1 reports the results under the first scenario with discrete jumps, and Table 2 for the second scenario involving transition periods. The results suggest that the proposed approach is able to detect essentially almost all the true change points under the first scenario and adequately detect the transition periods in the second scenario. Moreover, the proportion of false change points detected is close to 0 or negligible. In contrast, the detection of false change points is much higher under DCR, and it has poor performance in terms of detecting the true change points. Additional experiments (not presented here) reveal that the performance of DCR improves when the total number of time points in the experiment, along with the distance between consecutive change points is increased. However, for practical experiments with a few hundred time points, the DCR approach seems to fail in terms of change point estimation. The proposed approach, which incorporates SC knowledge and espouses a novel subnetwork sampling scheme, performs considerably better in terms of detecting the true change points while incurring minimal false positives. Figure 1(a)(b) provides a visual illustration of the change point detection performance under the first scenario, where the peaks in frequency are seen to concentrate around the true jumps.
In terms of graph estimation, the proposed approach has highest AUC compared with the other two methods. In addition, it has higher sensitivity and comparable or higher specificity than siGGM in most cases. Moreover, compared to DCR, the proposed method consistently has a significantly higher AUC value for network estimation. In addition, it has a higher or comparable sensitivity and consistently higher specificity, implying lower false positives in graph estimation and suitable power to detect true positives. In several cases, both the sensitivity and specificity under siCCPD were higher compared to DCR and siGGM. These results imply that the proposed method has much stronger power to detect true connections and is able to better control for false discoveries compared to siGGM for the vast majority of cases, which is expected when the true underlying network is dynamic, since siGGM is only designed to estimate static networks. Moreover, the difficulties inherent in the DCR approach are evident in experiments involving closely spaced change points and where the number of regions in the brain are not restricted to be small. In fact it was computationally infeasible to apply the DCR approach for the nodes example (reported below). Interestingly, all methods register a drop in specificity and an increase in sensitivity when the number of regions is increased from to , which suggests an increase in false positives under all approaches for higher dimensions.
Figure 1(b) presents the change point detection results for the highdimensional case of for data simulated under the ErdosRenyi network with three jump points. The Figure clearly depicts high peaks around the true change points, thereby suggesting that the proposed approach could successfully detect true change points based on the subnetwork sampling mechanism. Finally, the results for the case involving a higher number of true change points (10) with regions and time points, is presented in Figure 1(c), which clearly shows peaks under the subnetwork sampling scheme around all of the 10 true change points, thereby indicating the power of the proposed approach in detecting discrete jumps. In order to accurately detect all the true change points, we increased the number of subnetwork samples to 100 for this case. For this case, the DCR approach fails to detect an overwhelming majority of the jumps (results not reported).
In addition, we assessed the performance of the proposed approach when both the number of change points and ROIs increases. Simulation results (not presented here) shows that as long as the number of subnetworks is large enough and true change points are not exceedingly close together, our method successfully detects peaks around the true change points, under the subnetwork sampling scheme. We conjecture that an increasing number of subnetworks will be required for a good performance as the number of nodes in the network, as well as the number of true change points is increased. While this may likely increase the overall computation time, one can resort to parallel computing schemes across subnetwork samples for more efficient computations. Figure 1(d) provides the computation time per subsample as the subnetwork size is varied. In practice, subnetworks with 510 nodes is adequate for our approach.
In summary, using a combination of a powerful subsampling scheme and incorporating prior SC knowledge, the proposed method is shown to provide vast improvements over existing approaches in literature, and it is also shown to be scalable to a large number of nodes and change points, which illustrates its ability to adapt to wide variety of scenarios.
4 Analysis of Posttraumatic Stress Disorder Data
4.1 Description
Our study involves female AfricanAmerican participants from the Grady Trauma Project (GTP). These participants were recruited from primary care clinics at Grady Health System (GHS), a publicly funded, tertiary care center serving a predominantly socioeconomically disadvantaged innercity population in Atlanta, Georgia. A majority of these participants have experienced significant psychological trauma of various types (Gillespe et al., 2009). Imaging modalities including resting state functional magnetic resonance imaging (RsfMRI) and diffusion tensor imaging (DTI) data were collected for each individual. We focused on a smaller subset of 19 individuals who were aged below 30 years and hence we expected to have more homogeneous brain function and structure. We preprocessed DTI data for this subgroup of individuals and obtained SC strengths (described in Section 4.2). We use a whole brain parcellation corresponding to the time courses from the 264 ROIs under the Power system (Power et al., 2011) to perform our network analysis. Although the network analysis involved 264 nodes, we further group these ROIs into 10 functional modules as identified by Cole et al. (2013), which better characterize resting state functional networks. These modules included sensory/somatomotor, cinguloopercular, salience, auditory, subcortical, default mode, visual, frontoparietal, ventralattention, and dorsalattention. The coordinates for the ROIs and their allocation to these modules in presented in a Table in Appendix B.1. We note that 37 ROIs were excluded from our analysis based on Cole et al. (2013), either because they were located in cerebellum or they were identified with unknown functionality.
Clinical and Exposure Variables: Information was also collected on clinical and exposure variables such as childhood maltreatment using Childhood Trauma Questionnaire or CTQ (Bernstein et al., 1998) and adult trauma (via Trauma Exposure Index or TEI), as well as clinical measures such as DSMIV (which is a PTSD symptom scale or PSS (Foa, et al., 1993)). Our target clinical outcome is the ConnorDavidson Resilience Scale (CDRISC), which quantifies resilience (Connor and Davidson, 2003), i.e., the individual‘s ability to thrive when facing adversity. Resilience is a transdiagnostic indicator of mental health in the face of adversity, and is highly relevant to groups of individuals who have experienced high levels of trauma exposure and other forms of stress, as in the Grady Trauma Project sample. The CDRISC10 score has been shown to display good internal consistency and construct validity, and hence is a robust clinical measure of resilience. We investigate the potential of the resting state dynamic FC as a neuroimaging biomarker in terms of how well it can predict the CDRISC10 score, and compare its performance with siGGM and a SC naive version of our approach. Because of the high number of ROIs, DCR could not be directly applied here, but instead, comparisons were made with a SC naive version of the proposed approach that has similarities with DCR.
4.2 Data Acquisition and Preprocessing
Restingstate fMRI collected via a Siemens 3.0 Tesla Magnetom Trio TIM wholebody MR scanner(Siemens, Malvern, PA) included 146 brain volumes for each participant. The voxel size for the resting state fMRI acquisition was 3x3x3 mm with full brain coverage. The other related parameters were: TR/TE=2000/30ms, flip=, FOV=200x200 mm, slice thickness/gap = 4mm/mm. We applied several standard preprocessing steps for restingstate fMRI data. It included despiking, slice timing correction, motion correction, MNI(2mm) standard space registration, percent signal change normalization, linear trend removal, movement removal, bandpass filtering from 0.009 to 0.08. and spatial smoothing with a 6mm FWHM Gaussian Kernel. Finally, Power Atlas System was applied to aggregate the data into 264 ROIs. The average measurement for all voxels in each ROI was used as a representative measure for that ROI.
The DTI data was acquired on a Siemens 3T TIMTrio scanner at Emory University Hospital with the following acquisition parameters: 39x2.5mm thick axial slices, matrix=128x128, FOV=220x220 mm with voxel size = 1.72x1.72x2.5mm, with 60 directions, and a series of 4 b0 images. Standard preprocessing procedures, such as eddy current correction and biasfield correction were applied to the diffusion weighted data. Subsequently, we use the FSL functions bedpostx and probtracx2 to estimate the distribution of fiber tracts at each voxel and the count of white matter fibers tracts connecting all pairs of brain regions, respectively. In order to obtain the SC scores, we computed the average of the number of tracts reaching from the first to the second region, and from the second to the first region, divided by the total number of tracts sent out. Fiber tracks passing through gray matter or cerebrospinal fluid were discarded. These SC scores can be interpreted as the probability of structural connectivity between two ROIs, which is often refer to as SC strength.
4.3 Analysis
Our primary objective is to calculate dynamic FC while accounting for SC information for each subject, to investigate functional connections in the brain that undergo large dynamic changes in trauma exposed individuals. We calculated networks which were approximately around 15 percent density for all participants, which seem to reflect an acceptable sparsity level. The edges with the largest temporal fluctuations are identified as those that consistently switch over the different state phases, and are reported in Figure 3
. This Figure depicts the proportion of times each edge flips (changes from present to absent from one time point to the next, and viceversa) over the course of the scanning session, where a higher proportion implies greater temporal variability. We also performed a regression analysis to identify those edges whose temporal variability directly influence the resilience, with the candidate edges restricted to those functional modules that display the largest temporal fluctuations, as identified via Figure
3. The regression analysis involved resilience as the outcome and the explanatory variables were standard deviations of the edge strengths, computed in terms of the temporal variability of the edgespecific partial correlations. Due to a large number of edges involved, an univariate analysis was performed for one edge at a time, and significant effects were identified after multiplicity adjustment for pvalues using BenjaminiHochberg corrections (Benjamini and Hochberg, 1995). Edges having a significantly positive association under this analysis will imply connections where greater temporal variability enhances resilience and viceversa, whereas significant negative associations will imply connections where greater temporal changes leads to decrease in resilience and viceversa.
We also perform a secondary analysis to investigate the ability of the dynamic FC to differentiate individuals with distinct resilience levels. Instead of using edge level features, we use more global network metrics such as global efficiency and global clustering coefficients, as well as local clustering coefficient and local efficiency in the visual, salience, subcortical, ventral attention, and dorsal attention functional modules. We chose these modules since they were identified as regions with the highest FC changes for trauma exposed individuals under our dynamic network analysis (see Figure 3). We note that all the network metrics change across state phases under the dynamic network, and hence are timedependent. The network metrics represent the efficiency of information transmission at a global or local level in the brain, and were computed at each time point and for each subject using the Matlab toolbox Brain Connectivity Toolbox (Rubinov and Sporns, 2010a).
In order to predict resilience using these timevarying network summaries, we perform a scalaronfunction regression (Morris, 2015), which allows the associations between resilience and the timevarying features to also change over time. Hence, this model predicts resilience as a weighted sum of the timevarying features where the weights or regression effects are estimated from the data. We use the R package FDBoost (Brockhaus and Rgamer, 2016) to implement this approach. The performance of the scalaronfunction regression was assessed using goodness of fit ( values, higher is better) and predictive accuracy (out of sample mean squared error, lower is better), and compared with the SC naive version of the method as well as the siGGM approach. Since all the individuals included in the analysis were AfricanAmerican females aged between 1930 years, we did not include gender, race or age in our regression model.
4.4 Result
Figure 2 provides the histogram for the number of change points for the proposed approach and the SC naive version of the method. Our method detected 5 change points on average across all participants, with the number of change points ranging from 3 to 7. On the other hand, the SC naive version of the method registers only 1 change point for a large majority of participants, and only one subject has 3 change points. Given 146 brain volumes in the fMRI time series, the number of change points under the SC naive version seem to be unrealistic whereas the number of change points under the proposed method appears more practical. From Figure 3, we observe that connectivity in salience (SAL), subcortical (SCOR), frontoparietal (FPR), visual (VIS), ventral attention network (VAN), and dorsal attention network (DAN), consistently see the highest dynamic fluctuations across subjects, whereas the connections in other regions (including default mode network) see minimal or no temporal changes. These networks have been previously implicated in PTSD and trauma exposure studies (see Discussion). Interestingly, the highest temporal fluctuations in connectivity occur between the above functional networks, whereas the temporal fluctuations within most of these networks is somewhat limited.
The final panel in Figure 3 illustrate which dynamic FC within the SAL, SCOR, FPR, VIS, VAN, and DAN functional modules, have significant temporal variability that is related to resilience. The overwhelming majority of the connections whose temporal variability drive resilience are located between the FPR and VIS networks, the FPR and SAL networks and between the DAN and other networks. On the other hand, the temporal variability for connections within and between other modules have lesser direct impact on resilience. Moreover almost all the associations are negative, which means increased temporal fluctuations for these connections lead to decreased resilience, and viceversa. This is a novel finding that is unique to our dynamic connectivity analysis.
For our second objective, the boxplots for MSE and corresponding to the different network metrics are presented in Figure 4 in the article, and Figure 5 in Appendix B.2. respectively. From the Figures, it is clear that the proposed approach has a higher and lower MSE compared to the SC naive version and the siGGM approach for the almost all cases. A permutation test revealed that the improvements in MSE and were significantly better under siCCPD for all cases except panels (a), (l) and (o) in Figure 3, that correspond to the case of predicting resilience using global efficiency, local efficiency in the VAN, and local efficiency in the VAN and VIS modules combined, respectively. This suggests that while global and local efficiency may not always be able to differentiate the different approaches in terms of predicting resilience, the global and local clustering coefficient computed using dynamic resting state connectivity via the proposed approach is better in predicting resilience levels. The greatest gains in prediction accuracy are seen when predicting resilience using the local clustering coefficient in the VAN module, and the local efficiency in the VIS and DAN modules combined. These findings may indicate modules in the brain where dynamic connectivity changes have the greatest differentiating power with respect to resilience.
5. Discussion
In the current study we developed a novel method (siCCPD) for estimating dynamic changes in fMRI resting state functional connectivity guided by SC information in the brain. We demonstrated via extensive simulations, the clear advantages of the approach in estimating the dynamic network compared to alternate methods that do not account for SC information, and approaches that assume static connectivity while accounting for brain SC knowledge. The use of SC information to estimate dynamic changes in the brain network, results in greater power to detect true connectivity changes while reducing false discoveries. The approach can be scaled to a large number of change points as well as nodes, and is applicable to diverse settings.
We applied the method to the GTP study, to discover resting state network based alterations among participants exposed to varying degrees of trauma, and to predict psychological resilience based on the dynamic FC. We were able to accurately predict resilience using global and local clustering coefficients computed from the dynamic network, that may emerge as a potentially promising networkbased biomarker in trauma exposure neuroimaging studies. The prediction and goodness of fit results illustrate the clear advantages of siCCPD over the SC naive version that does not incorporate prior knowledge, and the siGGM method that incorporates prior knowledge to estimate a stationary network. The findings suggest that 1) including metrics of dynamic change in resting networks will improve models for predicting psychiatric risk and resilience to trauma and stress, and 2) estimating dynamic network connectivity models can be improved with the addition of DTIbased structural constraints. Regions of the brain that have the greatest predictive power for resilience were also identified, that may have implications in clinical research.
Regions in the SAL, SCOR, FPR, VIS, VAN, and DAN functional modules emerged as having the highest dynamic network changes across all participants. The SAL and SCOR modules have previously demonstrated altered resting state connectivity (Rabinak et al., 2011; Brown et al., 2014), and grey matter alterations were also discovered in subcortical areas for PTSD individuals (O’Doherty et al., 2017). Moreover, resting state connectivity differences in prefrontal cortex (Kennis et al., 2015) and hypoactivity in medial prefrontal cortex (Koenigs et al., 2009) have been observed. Finally, increased resting state FC has been discovered in dorsal and ventral attention networks (Block et. al, 2017), and differences in FC was also observed in the VAN when individuals were subject to noninvasive transcranial stimulation (Etkin et al., 2019). Hence our analysis is able to discover regions in the brain with large fluectuations in resting state FC that have been implicated previously in static connectivity studies. Unlike many of the previous approaches that used static connectivity, our analysis uses dynamic FC to assess these changes, and hence provides a unique perspective. Another key finding was that an increase in temporal fluctuations for a overwhelming proportion of connections within regions of the brain that display the highest temporal variability across participants, result in a decrease in resilience.
Future work may include deriving functional connectivity approaches that incorporate other ways of measuring structural connectivity in the brain, that may improve the power to detect dynamic changes even further. Another possibility is to propose models to better understand the relationships between SC and dynamic FC in the brain. It is important to note that although the proposed approach was applied to PTSD example, it is equally applicable to any resting state fMRI experiment where DTI data is also available for each participant.
6. Acknowledgements
This work was supported by National Institutes of Mental Health grant (R01MH120299, R01 MH071537, R21 MH098212); Georgia CTSA grant ULITR002378; National Institute of Child Health and Development (K12 HD085850); NIH National Centers for Research Resources (M01RR00039), National Center for Advancing Translational Sciences of the NIH (UL1TR000454).
References

Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., & Calhoun, V. D. (2014). Tracking wholebrain connectivity dynamics in the resting state. Cerebral cortex, 24(3), 663676.

Baker, A. P., Brookes, M. J., Rezek, I. A., Smith, S. M., Behrens, T., Smith, P. J. P., & Woolrich, M. (2014). Fast transient networks in spontaneous human brain activity. Elife, 3, e01867.

Barabasi, A. L., Albert, R., & Jeong, H. (1999). Meanfield theory for scalefree random networks. Physica A: Statistical Mechanics and its Applications, 272(12), 173187.

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57: 289300.

Bernstein, D. P., & Fink, L. (1998). Childhood trauma questionnaire: A retrospective selfreport: Manual. Harcourt Brace & Company.

Block, S. R., King, A. P., Sripada, R. K., Weissman, D. H., Welsh, R., & Liberzon, I. (2017). Behavioral and neural correlates of disrupted orienting attention in posttraumatic stress disorder. Cognitive, Affective, & Behavioral Neuroscience, 17(2), 422436.

Breslau, N. (2009). The epidemiology of trauma, PTSD, and other posttrauma disorders. Trauma, Violence, & Abuse, 10(3), 198210.

Brockhaus, S., & Rügamer, D. (2016). FDboost: boosting functional regression models. R package version 0.0, 16.

Brown, V. M., LaBar, K. S., Haswell, C. C., Gold, A. L., Workgroup, M. A. M., Beall, S. K., … & Green, K. T. (2014). Altered restingstate functional connectivity of basolateral and centromedial amygdala complexes in posttraumatic stress disorder. Neuropsychopharmacology, 39(2), 361.

Chang, C., Metzger, C. D., Glover, G. H., Duyn, J. H., Heinze, H. J., & Walter, M. (2013). Association between heart rate variability and fluctuations in restingstate functional connectivity. Neuroimage, 68, 93104.

Cole, M. W., Reynolds, J. R., Power, J. D., Repovs, G., Anticevic, A., & Braver, T. S. (2013). Multitask connectivity reveals flexible hubs for adaptive task control. Nature neuroscience, 16(9), 1348.

Connor, K. M., & Davidson, J. R. (2003). Development of a new resilience scale: The Connor‐Davidson resilience scale (CD‐RISC). Depression and anxiety, 18(2), 7682.

Cribben, I., Wager, T., & Lindquist, M. (2013). Detecting functional connectivity change points for singlesubject fMRI data. Frontiers in computational neuroscience, 7, 143.

Cribben, I., Haraldsdottir, R., Atlas, L. Y., Wager, T. D., & Lindquist, M. A. (2012). Dynamic connectivity regression: determining staterelated changes in brain connectivity. Neuroimage, 61(4), 907920.

Damoiseaux, J. S., & Greicius, M. D. (2009). Greater than the sum of its parts: a review of studies combining structural connectivity and restingstate functional connectivity. Brain structure and function, 213(6), 525533.

Eavani, H., Satterthwaite, T. D., Filipovych, R., Gur, R. E., Gur, R. C., & Davatzikos, C. (2015). Identifying sparse connectivity patterns in the brain using restingstate fMRI. Neuroimage, 105, 286299.

Etkin, A., MaronKatz, A., Wu, W., Fonzo, G. A., Huemer, J., Vértes, P. E., … & RamosCejudo, J. (2019). Using fMRI connectivity to define a treatmentresistant form of posttraumatic stress disorder. Science translational medicine, 11(486), eaal3236.

Erdos, P., & Renyi, A. (1960). On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5(1), 1760.

Foa, E. B., Riggs, D. S., Dancu, C. V., & Rothbaum, B. O. (1993). Reliability and validity of a brief instrument for assessing post‐traumatic stress disorder. Journal of traumatic stress, 6(4), 459473.

Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432441.

Gillespie, C. F., Bradley, B., Mercer, K., Smith, A. K., Conneely, K., Gapen, M., … & Ressler, K. J. (2009). Trauma exposure and stressrelated disorders in inner city primary care patients. General hospital psychiatry, 31(6), 505514.

Handwerker, D. A., Roopchansingh, V., GonzalezCastillo, J., & Bandettini, P. A. (2012). Periodic changes in fMRI connectivity. Neuroimage, 63(3), 17121719.

Hsieh, C. J., Dhillon, I. S., Ravikumar, P. K., & Sustik, M. A. (2011). Sparse inverse covariance matrix estimation using quadratic approximation. In Advances in neural information processing systems (pp. 23302338).

Higgins, I. A., Kundu, S., & Guo, Y. (2018). Integrative Bayesian analysis of brain functional networks incorporating anatomical knowledge. NeuroImage, 181, 263278.

Hinne, M., Ambrogioni, L., Janssen, R. J., Heskes, T., & van Gerven, M. A. (2014). Structurallyinformed Bayesian functional connectivity analysis. NeuroImage, 86, 294305.

Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., … & Handwerker, D. A. (2013). Dynamic functional connectivity: promise, issues, and interpretations. Neuroimage, 80, 360378.

Ioannidis, K., Askelund, A. D., Kievit, R., & Van Harmelen, A. L. (2016). The complex neurobiology of resilient functioning after child maltreatment.

Jia, H., Hu, X., & Deshpande, G. (2014). Behavioral relevance of the dynamics of the functional brain connectome. Brain connectivity, 4(9), 741759.

Jin, C., Jia, H., Lanka, P., Rangaprakash, D., Li, L., Liu, T., … & Deshpande, G. (2017). Dynamic brain connectivity is a better predictor of PTSD than static connectivity. Human brain mapping, 38(9), 44794496.

Kennis, M., Rademaker, A. R., van Rooij, S. J., Kahn, R. S., & Geuze, E. (2015). Resting state functional connectivity of the anterior cingulate cortex in veterans with and without post‐traumatic stress disorder. Human brain mapping, 36(1), 99109.

Kessler, R. C., Sonnega, A., Bromet, E., Hughes, M., & Nelson, C. B. (1995). Posttraumatic stress disorder in the National Comorbidity Survey. Archives of general psychiatry, 52(12), 10481060.

Koch, S. B., van Zuiden, M., Nawijn, L., Frijling, J. L., Veltman, D. J., & Olff, M. (2016). Aberrant resting‐state brain activity in posttraumatic stress disorder: A meta‐analysis and systematic review. Depression and anxiety, 33(7), 592605.

Koenigs, M., & Grafman, J. (2009). Posttraumatic stress disorder: the role of medial prefrontal cortex and amygdala. The Neuroscientist, 15(5), 540548.

Kundu, S., Ming, J., Pierce, J., McDowell, J., & Guo, Y. (2018). Estimating dynamic brain functional networks using multisubject fMRI data. NeuroImage, 183, 635649.

Li X, Lim C, Li K, Guo L, Liu T. (2013). Detecting brain state changes via fibercentered functional connectivity analysis. Neuroinformatics, 11(2), 193–210.

Lindquist, M. A., Xu, Y., Nebel, M. B., & Caffo, B. S. (2014). Evaluating dynamic bivariate correlations in restingstate fMRI: a comparison study and a new approach. NeuroImage, 101, 531546.

Messe, A., Rudrauf, D., Benali, H., & Marrelec, G. (2014). Relating structure and function in the human brain: relative contributions of anatomy, stationary dynamics, and nonstationarities. PLoS computational biology, 10(3), e1003530.

Morris, J.S. Annual Review of Statistics and Its Application 2015 2:1, 321359.

Ng, B., Varoquaux, G., Poline, J. B., & Thirion, B. (2012, October). A novel sparse graphical approach for multimodal brain connectivity inference. In International Conference on Medical Image Computing and ComputerAssisted Intervention (pp. 707714). Springer, Berlin, Heidelberg.

O’Doherty, D. C., Tickell, A., …, Bennett, M. R., & Lagopoulos, J. (2017). Frontal and subcortical grey matter reductions in PTSD. Psychiatry Research: Neuroimaging, 266, 19.

PinedaPardo, J. A., Bruña, R., Woolrich, M., …, Maestú, F., & Vidaurre, D. (2014). Guiding functional connectivity estimation by structural connectivity in MEG: an application to discrimination of conditions of mild cognitive impairment. Neuroimage, 101, 765777.

Power, J. D., Cohen, A. L., Nelson, S. M., Wig, G. S., Barnes, K. A., Church, J. A., … & Petersen, S. E. (2011). Functional network organization of the human brain. Neuron, 72(4), 665678.

Rabinak, C. A., Angstadt, M., Welsh, R. C., Kennedy, A., …, Martis, B., & Phan, K. L. (2011). Altered amygdala restingstate functional connectivity in posttraumatic stress disorder. Frontiers in psychiatry, 2, 62.

Rubinov, M., & Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage, 52(3), 10591069.

Shalev, A., Liberzon, I., & Marmar, C. (2017). Posttraumatic stress disorder. New England Journal of Medicine, 376(25), 24592469.

Shirer, W. R., Ryali, S., Rykhlevskaia, E., Menon, V., & Greicius, M. D. (2012). Decoding subjectdriven cognitive states with wholebrain connectivity patterns. Cerebral cortex, 22(1), 158165.

Smith, S. M., Miller, K. L., SalimiKhorshidi, G., Webster, M., Beckmann, C. F., … & Woolrich, M. W. (2011). Network modelling methods for FMRI. Neuroimage, 54(2), 875891.

Sporns, O. (2013). The human connectome: origins and challenges. Neuroimage, 80, 5361.

Sripada, R. K., King, A. P., Garfinkel, S. N., Wang, X., Sripada, C. S., Welsh, R. C., & Liberzon, I. (2012). Altered restingstate amygdala functional connectivity in men with posttraumatic stress disorder. Journal of psychiatry & neuroscience: JPN, 37(4), 241.

Stevens, J. S., Ely, T. D., Sawamura, T., …, Ressler, K. J., & Jovanovic, T. (2016). Childhood maltreatment predicts reduced inhibition‐related activity in the rostral anterior cingulate in PTSD, but not trauma‐exposed controls. Depression and anxiety, 33(7), 614622.

Thompson, G.J., Magnuson, M.E., Merritt, M.D., Schwarb, H., Pan, W.J.J., McKinley, A., Tripp, L.D., Schumacher, E.H.,& Keilholz, S.D. (2013). Short‐time windows of correlation between large‐scale functional brain networks predict vigilance intraindividually and interindividually. Human brain mapping, 34(12), 32803298.

van Rooij, S. J., Stevens, J. S., Ely, T. D., Fani, N., Smith, A. K., Kerley, K. A., … & Jovanovic, T. (2016). Childhood trauma and COMT genotype interact to increase hippocampal activation in resilient individuals. Frontiers in psychiatry, 7, 156.

Wang, H. (2012). Bayesian Graphical Lasso Models and Efficient Posterior Computation. Bayesian Anal. 7, no. 4, 867–886. doi:10.1214/12BA729

Watts, D. J., & Strogatz, S. H. (1998). Collective dynamics of ‘smallworld’networks. nature, 393(6684), 440.

Xu, Y., & Lindquist, M. A. (2015). Dynamic connectivity detection: an algorithm for determining functional connectivity change points in fMRI data. Frontiers in neuroscience, 9, 285.

Xue, W., Bowman, F. D., Pileggi, A. V., & Mayer, A. R. (2015). A multimodal approach for determining brain networks by jointly modeling functional and structural connectivity. Frontiers in computational neuroscience, 9, 22.

Yuan, M., & Lin, Y. (2007). Model selection and estimation in the Gaussian graphical model. Biometrika, 94(1), 1935.

Zhu, D., Li, K., Guo, L., Jiang, X., … & Wee, C. Y. (2012). DICCCOL: dense individualized and common connectivitybased cortical landmarks. Cerebral cortex, 23(4), 786800.
5 Appendix
Appendix A: Computation for siGGM
We use the following update steps corresponding to the th state phase () iteratively till we achieve convergence in the objective function. We denote in what follows.

Update the precision matrices: We update the precision matrix for the th stage phase in the th step as
(5) where is the empirical covariance matrix for the th bin, that is defined at the start of Section 2.1. The above can be solved using a approximation solver, QUIC (Hsieh et al., 2011), available in R.

Update the baseline effects via the closed form expression

Update via the closed form expression below where

Update when we got and with the following formula:
As there is no closed form solution, we implement a Newton Raphson solver to find the optimal value. The above formula could to reexpressed as:
where , denotes the upper diagonal elements of the structural connectivity matrix . and denotes the upper diagonal elements of and respectively. is the element wise exponential for each element of . We could only focus upon the upper diagonal elements of because is symmetric and we do not shrink diagonal elements.
Based on step size , the Newton Raphson updating equation is:where ,
and . is a diagonal matrix with elements as the upper triangular elements of , and similarly for . Based on simple calculation, is diagonal matrix and could be easily inverted. The step size () is searched using a back tracking line search for each update of as in Change et.al (2017).
Appendix B.1: Coordinates for ROIs in Power Atlas
ROI  X  Y  Z  Modules  ROI  X  Y  Z  Modules 

1  25  98  12  Exclude  46  66  8  25  Sensory/somatomotor 
2  27  97  13  Exclude  47  3  2  53  Cinguloopercular 
3  24  32  18  Exclude  48  54  28  34  Cinguloopercular 
4  56  45  24  Exclude  49  19  8  64  Cinguloopercular 
5  8  41  24  Exclude  50  16  5  71  Cinguloopercular 
6  21  22  20  Exclude  51  10  2  42  Cinguloopercular 
7  17  28  17  Exclude  52  37  1  4  Cinguloopercular 
8  37  29  26  Exclude  53  13  1  70  Cinguloopercular 
9  65  24  19  Exclude  54  7  8  51  Cinguloopercular 
10  52  34  27  Exclude  55  45  0  9  Cinguloopercular 
11  55  31  17  Exclude  56  49  8  1  Cinguloopercular 
12  34  38  12  Exclude  57  34  3  4  Cinguloopercular 
13  7  52  61  Sensory/somatomotor  58  51  8  2  Cinguloopercular 
14  14  18  40  Sensory/somatomotor  59  5  18  34  Cinguloopercular 
15  0  15  47  Sensory/somatomotor  60  36  10  1  Cinguloopercular 
16  10  2  45  Sensory/somatomotor  61  32  26  13  Auditory 
17  7  21  65  Sensory/somatomotor  62  65  33  20  Auditory 
18  7  33  72  Sensory/somatomotor  63  58  16  7  Auditory 
19  13  33  75  Sensory/somatomotor  64  38  33  17  Auditory 
20  54  23  43  Sensory/somatomotor  65  60  25  14  Auditory 
21  29  17  71  Sensory/somatomotor  66  49  26  5  Auditory 
22  10  46  73  Sensory/somatomotor  67  43  23  20  Auditory 
23  23  30  72  Sensory/somatomotor  68  50  34  26  Auditory 
24  40  19  54  Sensory/somatomotor  69  53  22  23  Auditory 
25  29  39  59  Sensory/somatomotor  70  55  9  12  Auditory 
26  50  20  42  Sensory/somatomotor  71  56  5  13  Auditory 
27  38  27  69  Sensory/somatomotor  72  59  17  29  Auditory 
28  20  29  60  Sensory/somatomotor  73  30  27  12  Auditory 
29  44  8  57  Sensory/somatomotor  74  41  75  26  Defaultmode 
30  29  43  61  Sensory/somatomotor  75  6  67  4  Defaultmode 
31  10  17  74  Sensory/somatomotor  76  8  48  15  Defaultmode 
32  22  42  69  Sensory/somatomotor  77  13  40  1  Defaultmode 
33  45  32  47  Sensory/somatomotor  78  18  63  9  Defaultmode 
34  21  31  61  Sensory/somatomotor  79  46  61  21  Defaultmode 
35  13  17  75  Sensory/somatomotor  80  43  72  28  Defaultmode 
36  42  20  55  Sensory/somatomotor  81  44  12  34  Defaultmode 
37  38  15  69  Sensory/somatomotor  82  46  16  30  Defaultmode 
38  16  46  73  Sensory/somatomotor  83  68  23  16  Defaultmode 
39  2  28  60  Sensory/somatomotor  84  58  26  15  Exclude 
40  3  17  58  Sensory/somatomotor  85  27  16  17  Exclude 
41  38  17  45  Sensory/somatomotor  86  44  65  35  Defaultmode 
42  49  11  35  Sensory/somatomotor  87  39  75  44  Defaultmode 
43  36  9  14  Sensory/somatomotor  88  7  55  27  Defaultmode 
44  51  6  32  Sensory/somatomotor  89  6  59  35  Defaultmode 
45  53  10  24  Sensory/somatomotor  90  11  56  16  Defaultmode 
ROI  X  Y  Z  Modules  ROI  X  Y  Z  Modules 

91  3  49  13  Defaultmode  136  4  48  51  Exclude 
92  8  48  31  Defaultmode  137  46  31  13  Defaultmode 
93  15  63  26  Defaultmode  138  10  11  67  Ventralattention 
94  2  37  44  Defaultmode  139  49  35  12  Defaultmode 
95  11  54  17  Defaultmode  140  8  91  7  Exclude 
96  52  59  36  Defaultmode  141  17  91  14  Exclude 
97  23  33  48  Defaultmode  142  12  95  13  Exclude 
98  10  39  52  Defaultmode  143  18  47  10  Visual 
99  16  29  53  Defaultmode  144  40  72  14  Visual 
100  35  20  51  Defaultmode  145  8  72  11  Visual 
101  22  39  39  Defaultmode  146  8  81  7  Visual 
102  13  55  38  Defaultmode  147  28  79  19  Visual 
103  10  55  39  Defaultmode  148  20  66  2  Visual 
104  20  45  39  Defaultmode  149  24  91  19  Visual 
105  6  54  16  Defaultmode  150  27  59  9  Visual 
106  6  64  22  Defaultmode  151  15  72  8  Visual 
107  7  51  1  Defaultmode  152  18  68  5  Visual 
108  9  54  3  Defaultmode  153  43  78  12  Visual 
109  3  44  9  Defaultmode  154  47  76  10  Visual 
110  8  42  5  Defaultmode  155  14  91  31  Visual 
111  11  45  8  Defaultmode  156  15  87  37  Visual 
112  2  38  36  Defaultmode  157  29  77  25  Visual 
113  3  42  16  Defaultmode  158  20  86  2  Visual 
114  20  64  19  Defaultmode  159  15  77  31  Visual 
115  8  48  23  Defaultmode  160  16  52  1  Visual 
116  65  12  19  Defaultmode  161  42  66  8  Visual 
117  56  13  10  Defaultmode  162  24  87  24  Visual 
118  58  30  4  Defaultmode  163  6  72  24  Visual 
119  65  31  9  Defaultmode  164  42  74  0  Visual 
120  68  41  5  Defaultmode  165  26  79  16  Visual 
121  13  30  59  Defaultmode  166  16  77  34  Visual 
122  12  36  20  Defaultmode  167  3  81  21  Visual 
123  52  2  16  Defaultmode  168  40  88  6  Visual 
124  26  40  8  Defaultmode  169  37  84  13  Visual 
125  27  37  13  Defaultmode  170  6  81  6  Visual 
126  34  38  16  Defaultmode  171  26  90  3  Visual 
127  28  77  32  Defaultmode  172  33  79  13  Visual 
128  52  7  30  Defaultmode  173  37  81  1  Visual 
129  53  3  27  Defaultmode  174  44  2  46  Frontoparietal 
130  47  50  29  Defaultmode  175  48  25  27  Frontoparietal 
131  49  42  1  Defaultmode  176  47  11  23  Frontoparietal 
132  31  19  19  Exclude  177  53  49  43  Frontoparietal 
133  2  35  31  Exclude  178  23  11  64  Frontoparietal 
134  7  71  42  Exclude  179  58  53  14  Frontoparietal 
135  11  66  42  Exclude  180  24  45  15  Frontoparietal 
ROI  X  Y  Z  Modules  ROI  X  Y  Z  Modules 

181  34  54  13  Frontoparietal  223  2  13  12  Subcortical 
182  21  41  20  Exclude  224  10  18  7  Subcortical 
183  18  76  24  Exclude  225  12  17  8  Subcortical 
184  17  80  34  Exclude  226  5  28  4  Subcortical 
185  35  67  34  Exclude  227  22  7  5  Subcortical 
186  47  10  33  Frontoparietal  228  15  4  8  Subcortical 
187  41  6  33  Frontoparietal  229  31  14  2  Subcortical 
188  42  38  21  Frontoparietal  230  23  10  1  Subcortical 
189  38  43  15  Frontoparietal  231  29  1  4  Subcortical 
190  49  42  45  Frontoparietal  232  31  11  0  Subcortical 
191  28  58  48  Frontoparietal  233  15  5  7  Subcortical 
192  44  53  47  Frontoparietal  234  9  4  6  Subcortical 
193  32  14  56  Frontoparietal  235  54  43  22  Ventralattention 
194  37  65  40  Frontoparietal  236  56  50  10  Ventralattention 
195  42  55  45  Frontoparietal  237  55  40  14  Ventralattention 
196  40  18  40  Frontoparietal  238  52  33  8  Ventralattention 
197  34  55  4  Frontoparietal  239  51  29  4  Ventralattention 
198  42  45  2  Frontoparietal  240  56  46  11  Ventralattention 
199  33  53  44  Frontoparietal  241  53  33  1  Ventralattention 
200  43  49  2  Frontoparietal  242  49  25  1  Ventralattention 
201  42  25  30  Frontoparietal  243  16  65  20  Exclude 
202  3  26  44  Frontoparietal  244  32  55  25  Exclude 
203  11  39  50  Salience  245  22  58  23  Exclude 
204  55  45  37  Salience  246  1  62  18  Exclude 
205  42  0  47  Salience  247  33  12  34  Exclude 
206  31  33  26  Salience  248  31  10  36  Exclude 
207  48  22  10  Salience  249  49  3  38  Exclude 
208  35  20  0  Salience  250  50  7  39  Exclude 
209  36  22  3  Salience  251  10  62  61  Dorsalattention 
210  37  32  2  Salience  252  52  63  5  Dorsalattention 
211  34  16  8  Salience  253  47  51  21  Exclude 
212  11  26  25  Salience  254  46  47  17  Exclude 
213  1  15  44  Salience  255  47  30  49  Sensory/somatomotor 
214  28  52  21  Salience  256  22  65  48  Dorsalattention 
215  0  30  27  Salience  257  46  59  4  Dorsalattention 
216  5  23  37  Salience  258  25  58  60  Dorsalattention 
217  10  22  27  Salience  259  33  46  47  Dorsalattention 
218  31  56  14  Salience  260  27  71  37  Dorsalattention 
219  26  50  27  Salience  261  32  1  54  Dorsalattention 
220  39  51  17  Salience  262  42  60  9  Dorsalattention 
221  2  24  30  Exclude  263  17  59  64  Dorsalattention 
222  6  24  0  Subcortical  264  29  5  54  Dorsalattention 
Comments
There are no comments yet.