Path analysis: A method to estimate altered pathways in time-varying graphs of neuroimaging data

The human brain is one of the most complex networks, comprising about 86 billion neurons connected by about 150 trillion synapses allowing neurons to communicate and pass chemical or electrical signals to each other (Azevedo et al., 2009; Farahani, Karwowski, & Lighthall, 2019; Fornito, Zalesky, & Bullmore, 2016). Advances in graph-theoretical analysis provide an opportunity to examine and understand this complex network in human cognition and behavior as well as neurological and psychiatric disorders (Farahani, Karwowski, & Lighthall, 2019; Fornito, 2016; Yu et al., 2015, 2018). Graph theory provides a simple and powerful way to model, estimate and simulate the structure and dynamics of brain networks (Falakshahi, Rokham, et al., 2020; Fornito, 2016; Fornito, Zalesky, & Bullmore, 2016; van den Heuvel & Fornito, 2014; van den Heuvel et al., 2017). A brain network or brain graph is composed of a set of nodes (vertices) as brain regions connected by a set of edges (links) as some measure of structural or functional interaction between those brain regions (Farahani, Karwowski, & Lighthall, 2019; Fornito, 2016).

The combination of graph theory and state-of-the-art noninvasive brain imaging techniques has shown promise in the detection of potential biomarkers of mental disorders (Farahani, Karwowski, & Lighthall, 2019; Yu et al. 2018). Brain function can be localized through functional magnetic resonance imaging (fMRI) that assesses the blood oxygenation level–dependent (BOLD) signal from the brain (Hillman, 2014). There is a vast literature on using graph–based methods to analyze brain functional (network) connectivity patterns applying fMRI data (Farahani, Karwowski, & Lighthall, 2019; Fornito, 2016; Karwowski, Vasheghani Farahani, & Lighthall, 2019; Yu et al., 2018). Different graph metrics were used in these studies to characterize the brain graphs, such as clustering coefficient, modularity, characteristic path length, small-world, and assortativity (Farahani, Karwowski, & Lighthall, 2019; Yu et al., 2018). For instance, clustering coefficient, local efficiency, and global efficiency have been shown to be lower in functional brain graphs in people with schizophrenia (SZ) compared to controls (Fu, Iraji, et al., 2021; Fu, Sui, et al., 2021; Liu et al., 2008). Topological alterations in basal ganglia and limbic systems in patients with sleep-related hyper-motor epilepsy have been reported using resting-state fMRI data (Evangelisti et al., 2018). There are numerous reports on graph-based abnormalities in autism spectrum disorder (ASD) using resting-state fMRI data that have shown lower modularity and clustering coefficient in ASD compared to controls (Di Martino et al., 2013; Farahani, Karwowski, & Lighthall, 2019; Keown et al., 2017; Rudie et al., 2013; van den Heuvel et al., 2017).

However, the focus in neuroimaging data has primarily been on edge and nodal changes or global graphic metrics, with little focus on the information contained in paths connecting brain regions. A path from node X to node Y in a graphical model is a sequence of adjacent edges between nodes X and Y. Paths are of great importance in brain graphs, and details of disruption of these paths in patient group graph may be highly informative for understanding disease mechanisms. Importantly, since different edges in the same path may drop out in different individuals, path-based analysis has the potential to capture information that is typically invisible to approaches that focus only on pairwise relationships.

Our main contribution in this study is to show that a comprehensive path (as opposed to individual edge) analysis may help identify putative path-based biomarkers from neuroimaging data. We compare and analyze paths between the brain graphs of control and patient groups. Interactions in graphs often happen in groups of three or more nodes that cannot be explained clearly in terms of pairwise relationships or single links. Paths analyses allow us to analyze the network beyond pairwise interaction. We provide an algorithm to estimate edges that trigger absent paths (disconnection) and additional paths (abnormal integration) in the patient group with reference to the control group graph. In addition, to examine path differences that are common across groups, we propose the use of a previously reported covariance decomposition method (Jones & West, 2005) in which the covariance between two nodes can be decomposed into a sum of weights associated with each of the paths connecting those two nodes.

We conduct our path analyses under the Gaussian graphical model (GGM) framework to have an interpretable model of brain graphs from resting-state fMRI data. GGM is frequently used to explore networks among multiple random variables and, in this study, represents the interaction between brain components. In a GGM, two variables are conditionally independent given all other variables, if and only if their corresponding off-diagonal element in the precision matrix (inverse of the covariance matrix) is zero, and there is no edge between them. Therefore, the precision matrix in GGM summarizes the conditional dependence of network structure (He et al., 2019).

Previous studies for network analysis, for the most part, focused on Pearson correlation–based matrices to measure the strength of association between nodes in a network. GGM is used to study the partial correlation structure of a set of variables and is a better way to model a complex network because of its interpretation with conditional dependence between two random variables after removing the effect of all other controlling variables. Pearson correlation, on the contrary, ignores the effect of all other variables (He et al., 2019). There are many methods for estimating GGM. In this study, we chose to use the joint graphical lasso, which is a technique for jointly estimating multiple graphical models corresponding to distinct but related classes (Danaher, Wang, & Witten, 2014; Guo, Levina, Michailidis, & Zhu, 2011; Shan & Kim, 2018). The joint estimation borrows strength across classes to boost the estimation of multiple graphical models that share certain characteristics while retaining support for class-specific structures. Additionally, estimating control and patient graph separately using the group-specific data does not consider similarity between the group models. Jointly estimating control and patient group graphs improves identifying the differences between groups by considering the assumption of similarity between the models.

Recent brain imaging studies have demonstrated that brain functional connectivity is dynamic over time (Allen et al., 2014; Bhinge, Long, Calhoun, & Adali, 2019; Damaraju et al., 2014; Iraji et al., 2020). However, prior brain graph–based studies focused predominantly on static functional (network) connectivity with an implicit assumption of stationary brain interactions during the scanning period. Although static functional connectivity has disclosed a great deal of knowledge, evaluating brain connectivity changes across time can reveal additional valuable insight into the inherent dynamic connectivity of the human brain. Therefore, in addition to static, we analyze paths on time-varying GGMs for control and patient groups. To the best of our knowledge, no study has provided a dynamic path–based comparison of time-varying brain graphs in group studies to identify possible brain abnormalities.

We illustrate the utility of the proposed method on resting-state fMRI data of SZ patients. SZ is commonly distinguished as a disorder of brain connectivity (Fornito, Zalesky, Pantelis, & Bullmore, 2012), although its diagnosis is mostly dependent on qualitative symptom–based measures (Rokham et al., 2020). It is characterized by functional disconnectivity (Friston, 1998) or abnormal integration between distant brain regions (Damaraju et al., 2014; Skudlarski et al., 2010). Prior studies have shown functional connectivity abnormalities in SZ in a broad spectrum of systems that frequently involve prefrontal brain regions (Cole, Anticevic, Repovs, & Barch, 2011; Fornito et al., 2013; Skaatun et al., 2017; van den Heuvel & Fornito, 2014). Still, little is known about how this disconnection or abnormal integration manifests. We estimate time-varying and static graphs of the control and SZ group and apply our proposed method for comprehensive path analysis. A preliminary version of this work has been reported in Falakshahi, Vergara, et al. (2020).

The remainder of the paper is structured as follows: Materials and Methods describes the details of our method, including estimating static and time-varying graphs of control and patient groups and our proposed algorithm for path analysis; Results provides details of our path analysis results on resting-state fMRI data of individuals with SZ and control groups; Discussion reviews our results and implications. We provide concluding remarks in Conclusion.

The methods consist of two parts: identifying static and dynamic graphs from resting-state fMRI data and “path analysis,” including identifying edges that are associated with disconnection or abnormal integration in patient group graph and the novel use of a covariance decomposition method to examine the case where there is at least one path between two specific brain components in both control and patient groups. In Figure 1, we illustrate the main steps, and each section will be explained in detail in the following subsections.

In this work, we analyzed eyes-closed resting-state fMRI data collected from 160 controls and 151 age, gender, and mean framewise displacement (FD)–matched individuals with schizophrenia (SZ) (age: p = 0.18; gender: p = 0.39; mean FD: p = 0.97) by the Function Biomedical Informatics Research Network (fBIRN) (Keator et al., 2016). fBIRN demographics are available in Supporting Information Table S1. Written informed consent was obtained from all subjects. Data were collected using 3T scanners with a repetition time (TR) of 2 s, voxel spacing size of 3.44 × 3.44 × 4 mm, a slice gap of 1 mm, and a total of 162 volumes at seven different sites. Data for six of the seven sites were collected on a 3T Siemens Tim Trio System, and for one other site, a 3T General Electric DiscoveryMR750 scanner was used. Further details on the dataset can be found in Fu, Iraji, et al. (2021) and Fu, Sui, et al. (2021).

Data preprocessing was performed using the SPM 12 (https://www.fil.ion.ucl.ac.uk/spm/) toolbox. We performed rigid body motion correction to correct subject head motion, followed by the slice-timing correction to account for timing difference in slice acquisition. The fMRI data were normalized to the EPI template and resampled to 3-mm³ isotropic voxels. Data were spatially smoothed using a Gaussian kernel with a 6-mm full width at half maximum. Subjects with head motion < = 3° and < = 3 mm, and with functional data providing near full brain successful normalization (by comparing the individual mask with the group mask) were selected for further analysis. More details on preprocessing and subject selection can be found in Damaraju et al. (2014), Fu, Iraji, et al. (2021), and Fu, Sui, et al. (2021).

We use independent component analysis (ICA) to determine the nodes of a graph (Yu et al., 2017), which decomposes the whole-brain fMRI into independent spatial components known as ICN. Each ICN is composed of a set of voxels that shows a strong coupling of spontaneous fluctuations in the BOLD signal that can be considered as one functional unit (Iraji et al., 2019). Using fully blind ICA as a method for ICN extraction estimation of functional connectivity measures might result in different components identified across data. This inconsistency of identified components may impede finding replication or comparison. To address this and obtain the same set of ICNs for both control and SZ groups, we used the fully automated Neuromark pipeline, which is a spatially constrained ICA informed by reliable network templates (Du et al., 2020). To generate the Neuromark templates and to identify robust and reproducible ICNs, spatial ICA was applied to two large datasets of typical controls (Human Connectome Project [HCP; 823 subjects after the subject selection] and Genomics Superstruct Project [GSP; 1,005 subjects after the subject selection]). Group ICA with model order 100 was performed on each dataset separately to obtain group-level components. Then independent components from HCP and GSP were matched based on the spatial similarity (correlation > = 0.4), and reproducible ICNs were chosen as the network templates. ICNs identification and functional domain labeling were performed by five neuroscience experts. More details of the Neuromark templates can be found in Du et al. (2020). Next, spatial maps and time courses (TCs) for the fBIRN dataset were obtained by the network templates as prior information within a spatially constrained ICA algorithm. For the spatially constrained ICA algorithm, we used multivariate objective optimization independent component analysis with reference (MOO-ICAR), implemented in the GIFT toolbox from TReNDS (https://trendscenter.org/software/gift) (Iraji et al., 2020), as it shows good performance based on our previous studies (Du & Fan, 2013; Du, Allen, et al., 2016).

To create a static functional network connectivity (sFNC) matrix for each subject, pairwise correlations were calculated using the entire length of the ICNs time course. We then calculated the mean of sFNC matrices across subjects for control and patient groups. Further analysis on sFNC, including group differences using univariate statistical test, estimating group graphs and path analysis described in sections D and E.

Dynamic functional network connectivity (dFNC) investigates the time-varying interactions between brain networks and has been studied extensively (Bhinge, Long, Calhoun, & Adali, 2019; Sakoğlu et al., 2010). dFNC analysis enables us to evaluate how FNCs between ICNs evolve over time. A sliding window approach with square window type was applied to the selected ICNs’ TCs to study the connectivity between selected ICNs in a dynamic manner. ICNs’ TCs were localized by a sliding window size of 20 TR and a window step size of 5 TR to reduce computational demands, which led to 28 windows capturing changes in connectivity over time.

Next, we used the joint graphical lasso (Danaher, Wang, & Witten, 2014) to estimate the inverse of the covariance matrix (precision matrix) for windowed TCs of each subject. We chose to use a joint graphical lasso, which is a technique for jointly estimating multiple graphical models corresponding to distinct but related classes (windows of each subject) (Danaher, Wang, & Witten, 2014; Guo, Levina, Michailidis, & Zhu, 2011; Shan & Kim, 2018). As the simulation study suggested in Danaher, Wang, and Witten (2014), we resample each window data point by using the fast Fourier transform method to augment the sample time from 20 to 400 to robustly estimate covariance matrices and improve the accuracy of estimation. Having the estimated covariance matrix, we then computed the correlation matrix for each window; the detailed formulation is provided as Supporting Information Formula 1.

Recent resting-state studies have shown that fluctuations in the brain networks are not random and reoccur across and within subjects (Damaraju et al., 2014). These dynamic brain states or patterns can be identified by grouping windowed FNC using a k-means clustering algorithmv (Allen et al., 2014; Bhinge, Long, Calhoun, & Adali, 2019; Iraji et al., 2020). Therefore, each cluster contains similar functional connectivity patterns. We identified these functional brain states from computed dFNC matrices across all subjects. First, we flatten the upper diagonal part of each FNC matrix to obtain a feature vector. Then, for each subject, the variance across the feature was computed. Similar to what has been done previously (Allen et al., 2014; Damaraju et al., 2014), initial clustering was applied on a subset of data exhibiting maximal deviation from the mean called subject exemplars. The clustering algorithm was applied to these subject exemplars by using the L1 distance function (Manhattan distance) as the distance measure with random initialization of centroid positions and replicated 50 times to increase the chances of escaping local minima. As the k-means algorithm is very sensitive to the starting point, we applied subject exemplars to select the peaks (strongest contributors) as initialization. Manhattan distance was chosen as it has been suggested that it might be more effective to estimate similarity measures for high-dimensional data (Allen et al., 2014). The centroids from clustering the subject exemplars were used as initial points to cluster all dFNC windows from all subjects (Allen et al., 2014; Damaraju et al., 2014). The optimal number of clusters was determined using the elbow criterion.

We performed the univariate test in static and dynamic functional (network) connectivity to evaluate differences between the control and patient groups. Regarding sFNC, we applied a two-sample t test in a univariate manner on estimated sFNC matrices. To obtain group differences in dFNC, we first computed a subject mean for each state from the subject FNC windows assigned to that state as a representative pattern of connectivity of the subject for that state. Then, to identify significant differences between patient and control groups, we performed t tests per ICN pair of the averaged dFNC for each state. We corrected for multiple comparisons by using the FDR with a threshold of 0.05.

We used GGM for modeling the brain and identifying the graph structure of control and patient groups. The brain graphs gained from GGM are used for encoding relationships among brain components, wherein nodes represent ICNs and edges that demonstrate a partial correlation between nodes. In a GGM, two nodes are conditionally independent, given all other nodes if and only if their corresponding off-diagonal element in the precision matrix is zero. To estimate the precision matrix, we used joint graphical lasso (Danaher, Wang, & Witten, 2014) on ICNs’ time courses.

To estimate static group graphs of control and patient, we concatenated ICNs’ time courses across the subjects of control and patient groups separately and estimated the precision matrix by performing joint graphical lasso estimation on concatenated ICN time courses. To estimate time-varying group graphs of control and patient groups, we extracted a GGM for each state. Clustering the windowed dFNC results in each windowed dFNC being assigned to a specific cluster. Next, the time points of the windows located in one cluster were concatenated to estimate group graphs for each state. Using the concatenated ICN time courses of corresponding windows in each state, we jointly estimated precision matrices for control and patient groups applying the joint estimation presented in Danaher, Wang, and Witten (2014).

Having the precision matrices for each group in static and dynamic states, we computed the partial correlation matrices. The detailed formulation is provided as Supporting Information Formula 2. To obtain the adjacency matrices, the parametric test for statistical significance of the partial correlation was applied for each element of the partial correlation matrix to determine the significant edges. The edge was considered between two ICNs only where the corresponding FDR corrected p value was significant (p < 0.05).

We compared paths between the graphs of control and patient groups for estimated static and time-varying in four possible cases:

1. Disconnection: In this case, two ICNs are reachable from each other, which means there is at least one path between them in the control group graph, but there is no path between them in the patient group graph.

2. Abnormal integration: In contrast to case 1, case 2 occurs when there is no path between two specific ICNs in the control group graph, but there is at least one path between those ICNs’ pair in the patient group graph.

3. Connection in both control and SZ patients: In this case there is at least one path between specific pairs of ICNs in both control and SZ patient group graphs.

4. Disconnection in both control and SZ patients: The last case is when there is no path between two ICNs in both control and patient group graphs.

We considered case 4 as a trivial case and focused on cases 1–3. For cases 1 and 2, we proposed an algorithm to estimate edges associated with disconnection and abnormal integration by mimicking the structure of paths in the control group graph. For case 3, to make the comparison between control and patient graphs, we applied the covariance decomposition in GGM (Jones & West, 2005). We first describe our proposed method for cases 1 and 2 in detail, and then we discuss the proposed method for case 3.

Our proposed method can be applied to estimated patients and control group graphs obtained through GGM to identify which links are associated with mental illness across the fMRI dataset. There are some edges in the control group graph that contribute to creating path(s) between some nodes (ICNs). However, the absence of those edges in the patient group graph leads to the absence of some paths between some nodes. We labeled those edges as disconnectors that are associated with disconnection in the patient group graph. On the contrary, there are some additional edges in the patient group graph that contribute to creating some new paths that do not exist in the healthy group graph. We labeled those edges as connectors that trigger abnormal integration in the patient group graph.

To compare the paths between nodes of estimated graphs, we utilized the concept of the connected component in graph theory. A connected component of an undirected graph is a subgraph in which there is a path between every pair of nodes in that subgraph (Balakrishnan & Ranganathan, 2012). Hence, when two nodes are in the same connected component, it can be concluded that there is at least one path between them.

A precise and succinct description of our method is summarized in Algorithm 1, which receives two graphs as inputs, Graph1 and Graph2 (Line 1 of pseudo-code). The order of the inputs is important. Disconnectors associated with disconnectivity can be obtained by setting the Graph1 and Graph2 inputs as control and patient group graphs, respectively. On the contrary, abnormal integration can be obtained by swapping the inputs and setting Graph1 and Graph2 as patient and control group graphs, respectively.

Algorithm 1 first identifies all connected components in Graph1 and Graph2 (lines 3 and 4 of pseudo-code). Then, for each connected component in Graph1, it checks whether the nodes of that connected component are spread into multiple connected components in Graph2 or not. One way to check this is to see if the nodes set of that connected component in Graph1 intersects with more than one connected component nodes sets in Graph2 (lines 6–12 of pseudo-code). The “Container” set (line 11 of pseudo-code) keeps those connected components in Graph2 that intersects with the current connected component of Graph1 (c1). Then, if the cardinality of the Container set becomes greater than or equal to two (line 12 of pseudo-code), it indicates that the nodes of the current connected component of Graph1 (c1) are observed in at least two connected components in Graph2, and it reveals that the separation between the nodes of c1 occurred in Graph2. Next, to identify edges associated with disconnection, all combinations of two connected components in the Container will be investigated (lines 13–23 of pseudo-code). Then for each 2-combination of connected components in the Container set, the union of their nodes is stored in set “N” (line 14 of pseudo-code). Using node-set N, induced subgraphs of Graph1 and Graph2 will be formed (lines 15–16 of pseudo-code). Next, differences between the edge sets of these two subgraphs will be obtained (set “E”) (line 17 of pseudo-code). As not all the edges in set E cause separation, from the set E, those edges will be considered as edges associated with disconnection if their end points belong to two distinct connected components in Graph2. They also should be part of the current connected component in Graph1 (c1) that we are currently analyzing (line 18–23 of pseudo-code). If set E becomes empty, it indicates that the two connected components in Graph2 have become separated indirectly through another connected component.

Edges associated with abnormal integration (connectors) can be identified similarly, by swapping the order of input variables of the algorithm. In Figure 2, edge (6,8) shown with green is associated with abnormal integration.

To examine case 3, where two distinct brain components are connected and reachable from each other in both control and patient groups, we first need to introduce a theorem proposed in Jones and West (2005).

We applied the covariance decomposition approach for case 3 where paths exist between specific brain components in both control and patient groups to obtain covariance weights for each path between two nodes in both groups. To be able to compare different paths between groups, we used alternative correlation weights, which are constant multiples of the covariance-based weights (Jones & West, 2005). Unlike the covariance-based weights, the correlation-based weights are comparable between paths with different end points.

To identify the differences between control and patient group in case 3, we computed all the paths between each pair of ICNs in both groups and calculated their corresponding correlation weights by using the decomposition method. Lastly, we investigated which paths took unique trajectories in control or patient graphs and studied how their corresponding correlation weights contribute over the whole correlation between the two end points. To this end, for each pair of ICNs that are connected in both control and patient groups, we obtain correlation weights for all simple paths between them via the decomposition method. Next, we identified those brain components that took a distinct path trajectory in comparison with other group and have a significant contribution (more than 50%) to the correlation between paths end points.

We applied the proposed method to resting-state fMRI data collected from control and schizophrenia groups. We selected 53 out of 100 as ICNs for further analysis and categorized them into seven functional domains. The seven labeled functional domains comprise the auditory, subcortical, sensorimotor, visual, cognitive control, cerebellar, and default mode as depicted in Figure 4, and ICN numbers, labels, and domains are presented in Table 1.

To identify significant differences between patient and control groups, we performed two-sample t tests per ICN pair of sFNC matrices. We corrected for multiple comparisons by using the FDR with a threshold at 0.05. Figure 5A shows the group-specific mean of sFNC matrices and the group differences results. The sFNC analysis results show that the schizophrenia group had weaker functional connectivity between sensory domains (AUD, SM, VIS) in comparison with the control group. Similarly, the connectivity between subcortical and cerebellar domains was lower in the patient group. However, the schizophrenia group showed stronger connectivity between sensory domains and subcortical domain. In addition, the connectivity between sensory domains and cerebellar domain was stronger in schizophrenia group.

Within subject, windowed FNC matrices were computed as described earlier. We computed the variance across the windowed FNC matrices for each subject and selected windows corresponding to local maxima that resulted in 7.2 ± 1.3 (mean ± SD) windows per subject with a range of 3 to 11; k-means clustering was applied to subjects exemplars, and the centroids from this clustering were used as initial points to cluster all dFNC windows from all subjects. The optimal number of clusters was determined as five by the elbow criterion, which is within a reasonable range of previous dFNC studies (Allen et al., 2014; Fiorenzato et al., 2019; Fu et al., 2018; Tu et al., 2019). The elbow plot can be seen in Supporting Information Figure S1.

The average of functional connectivity for each state are shown in Figure 6A. In group difference evaluations, the number of significant FDR corrected p values were 16 for State 1, 44 for State 2, 19 for State 3, 37 for State 4, and 30 for State 5. Figure 6A also illustrates the results of the multiple comparison test for the control and SZ groups.

State 2 is the weak-connected state among all states as it shows the weak correlation between all domains. States 1, 4, and 5 have moderate connectivity and State 3 is a strong-connected state. States 5, 4, and 3 have moderate to high functional connectivity among sensory domains (AUD, SM, VIS). The antagonism correlation between SC and sensory domains increases from State 4 to 5 and reaches to its highest value in State 3. Cell-wise differences between each cluster median (FC state) have been illustrated in Supporting Information Figure S2.

Table 2 shows the occupancy measure of a given state as the percentage of windowed FNC that had been labeled with the cluster represented by the given state. About 52% of the windowed dFNCs were clustered as State 2 and State 4, and the occupancy of control group for State 2 and State 4 were 44% and 51%, respectively. In general, the occupancy of different states shows the SZ group were mostly located in weakly connected states (State 1, State 2) in comparison with control, which were observed in strong-connected states (State 3). State 1 was occupied mostly by SZ patients such that 78% of windowed dFNCs of the State 1 belonged to SZ patients, and the occupancy of control group in that state was only 22%. However, 74% of windowed dFNCs in State 3 belonged to the control group. Similarly, the occupancy level of SZ for State 5 was 26% and was lower than the control group. The connectivity between SM and VIS domain in State 3 was stronger than in State 1. In addition, the connectivity between SC and CC, DM, and CB in State 3 was higher than State 1. However, the connectivity between sensory (AUD, SM, VIS) and CB, and between VIS and SB in State 1 was higher than State 3.

We estimated static brain graphs for control and patient groups by using GGM. Figure 5B illustrated the estimated static graphs of control and patient group. The estimated graphs for each state (time-varying graphs) of control and patient groups are illustrated in Figure 6B.

We applied our proposed algorithm to identify missing edges that are associated with disconnection and additional edges that are associated with abnormal integration. In static graph analysis, we identified six missing edges and two additional edges that trigger absent paths and additional paths in the SZ group, both within and between functional domains. We found four missing links associated with blocked paths within visual, subcortical, and cognitive control networks and one missing link between cognitive control and auditory (Insula, Middle temporal gyrus) and one between subcortical and cognitive control (Thalamus, Hippocampus). Two additional edges identified between subcortical, auditory, and visual networks. Figure 5C illustrates differential graphs obtained by our proposed method. Table 3 shows details of missing and additional links (edges) associated with SZ in our static analysis.

The disconnectors and connectors that were identified through dynamic approach are shown in Figure 7, and Table 4 shows details of these edges. In State 1, two missing and three additional edges were identified within the default mode and visual domain. The results show that the connection between the lingual gyrus and the calcarine gyrus in the SZ group was disrupted, and the lingual gyrus and the right middle occipital gyrus and the inferior occipital gyrus connected with new additional edges.

In State 2, two missing links were observed within the cognitive control domain. In addition, one missing link was identified between the insula of the cognitive control domain and the middle temporal gyrus of the auditory domain. In addition, three new connectors were observed in the SZ group within ICNs of the cognitive control, the visual, and the sensorimotor domains.

In State 3, different connections within the subcortical domain were observed in the SZ group. The connection between the caudate and the thalamus and the subthalamus/hypothalamus were disrupted in the SZ group, and new connection between the caudate and the putamen were identified in the SZ group through an additional link. An additional new connection in the SZ group was observed within the sensorimotor domain between the left and the right postcentral gyrus.

In State 4, two disconnectors were identified within the cerebellar and the default mode domains. Five additional links were observed within the cognitive control, the visual, the default mode, and the sensorimotor domains. Also, one additional link was identified between the cognitive and the default mode network domains. The missing links triggered separation within some ICNs of the default mode and the cerebellar domains. However, the additional edges created more connections within the sensorimotor, the visual, and the cognitive control ICNs. Three new additional edges of (superior frontal gyrus, inferior parietal lobule; posterior cingulate cortex, inferior parietal lobule; and precuneus, precuneus) were created for new additional paths between nodes of the cognitive control and the default mode domains in the SZ group.

In State 5, four disruptions within the default mode, the cognitive control, and the auditory domains and one disruption between the cognitive control and the subcortical domains were observed in the SZ group. Three disconnectors were identified that triggered disconnections between nodes of a connected component in the control group that belongs to the default mode, the cognitive control, and the subcortical domains. Two additional integrations also were observed within the sensorimotor and the default mode domains in the SZ group.

We also examined the case that paths exist between two specific brain components in both groups by using the covariance decomposition method. In static analysis of control and SZ group graphs, 67 pairs of ICNs were identified between which paths exist in both control and SZ groups. Ten of these pairs took unique path trajectories in at least in one group and have a contribution of more than 50% of the correlation between path end points obtained via the decomposition method. Table 5 summarize the results of those ICN pairs in static analysis. Results of static analysis showed significant differences within the sensory motor, the subcortical, and the visual domains. The SZ groups are connected through higher number of paths that took unique trajectories, with higher correlation weight in comparison with the control group in the sensory motor and the subcortical domains. However, within the visual domain the control group showed a higher number of paths with higher correlation weights.

In time-varying graphs of control and SZ analysis, 127 pairs of ICNs were identified between which paths exist in both control and SZ groups. States 3 and 5 revealed that pair 4 and 13 took unique path trajectories at least in one group, respectively, which have a contribution of more than 50% of the correlation between path end points obtained via decomposition. Table 6 summarizes the results of those ICN pairs. Within the visual domain in State 3, the SZ group showed a higher number of paths than the control group in all four cases. Also, distinct paths in terms of the sequence of the nodes in the SZ group account for more than 50% of the total correlation weight between path end points. For example, the number of paths between the calcarine gyrus (ICN17) and the right middle occipital gyrus (ICN 21) was six in the SZ group, while the number of paths between these two ICNs was two in the control group. Two paths in the SZ group were common to the control group, and the other four took some unique path trajectories, and decomposition revealed that those distinct paths in the SZ group have a significant contribution (66.39%). As an example, Figure 8 depicts the result of covariance decomposition analysis between the right middle occipital gyrus and the lingual gyrus in State 3. The correlation between these two ICNs in the control and SZ groups is 0.0046 and 0.0772, and the number of paths between them in the control and SZ groups is two and four, respectively. The decomposition approach shows that two new paths in the patient group contribute 93.61% of correlation weight. The other two paths in SZ group that took the same trajectories as the control group contribute 6.39% of correlation weight.

Similar scenarios were observed within the subcortical and the visual domain in State 5. Four cases show that the SZ group has more paths between components of putamen, caudate, and thalamus, which also involve a significant proportion of total correlation between those end points. In addition, two cases that involve the inferior occipital gyrus as one of the path end points show that the distinct paths in the SZ group had a significant effect on the total correlation between the path endpoints.

Five cases in State 5 show the paths from the middle occipital gyrus to other visual domain components have totally different trajectories in control and SZ groups and do not share any common paths between control and SZ groups. Moreover, the path analysis between the cuneus and the right middle occipital gyrus shows that the SZ group has more paths than the control group; however, the unique path of control significantly contributes to the total correlation weights between end points. Similar results were observed between the cuneus and the inferior occipital gyrus, but the distinct paths of both control and SZ groups had a significant contribution to the decomposition.

In this study, we estimated static and dynamic brain graphs from resting-state fMRI data. We applied a data-driven approach for determining the nodes and edges of the graphs. Spatial brain components evaluated by group ICA was used to define data-driven nodes, and for estimating graphs edges, we employed statistical tests and mitigated the bias of relying on ad hoc thresholds for correlation values. After estimation of a partial correlation matrix, we obtain the corresponding p values, and we only considered an edge between two nodes if it passes a significance level (0.05) after multiple comparison correction (FDR). This contribution reduces the need for subjective choices as in other widely used approaches (Karwowski, Vasheghani Farahani, & Lighthall, 2019).

Existing methods for identifying differences between control and SZ patient groups in functional brain imaging literature are limited to correlation networks and pairwise functional connectivity comparison. Most studies focus only on pairwise relationships or “single links.” The path in such networks cannot reveal meaningful interpretation in comparison with the GGM, in which the conditional dependency plays an important role in determining connection and paths between nodes.

We analyzed paths on both static and dynamic manner of control and patient group graphs. We provide an algorithm to estimate edges associated with disconnection (absence of paths) or aberrant integration (additional paths) in dynamic patient group graphs by mimicking the structure of control group graphs. We used the concept of the connected component in graph theory to examine the presence or absence of a path between two nodes of a graph, as there is at least one path between every pair of nodes in a connected component. This results in time complexity reduction of path analysis in comparison with checking the absence or presence of paths for every pair of nodes of a graph. In the first part of the algorithm, we identified connected components in graphs using the NetworkX python package (Hagberg, Swart, & Chult, 2008). Next, through multiple steps in the proposed algorithm, we checked if the nodes of each connected component separated or joined to a new connected component in the second graph. Then, if any separation or abnormal integration occurred, the edges associated with those cases would be returned as an output of the algorithm. This algorithm can be applied to either static or dynamic graphs.

In addition, we studied the cases where two nodes (ICNs) connected in both control and patient groups. We investigated which paths take unique trajectories and computed their corresponding correlation weights by using the covariance decomposition method in a GGM. We applied our method to both static and dynamic time-varying graphs to analyze changes in paths over time.

We illustrated the utility of our proposed method on resting-state fMRI data of SZ patients. Previous studies have been demonstrated that SZ is associated with a disruption of the connections present in the healthy brain (Hunt, Kopell, Traub, & Whittington, 2017; Pearlson, 1997; van den Heuvel & Fornito, 2014; van den Heuvel, Mandl, Stam, Kahn, & Pol, 2010). Since SZ is a psychiatric disorder that can be distinguished by functional disconnectivity or abnormal integration between distant brain components (Damaraju et al., 2014; Friston, 1998; Skudlarski et al., 2010), we analyzed additional and blocked paths in estimated SZ time-varying graphs. However, the proposed method can be applied to any data related to other conditions that can be distinguished by disconnectivity or abnormal integration. Using our proposed algorithm, we found several disconnectors and connectors associated with disconnectivity and abnormal integration (see Table 3 and Table 4).

Our findings show 14 disconnectors and 16 connectors, both within and between functional domains, across all five states in time-varying brain graphs analysis. Static approach reveals six disconnectors and two connectors. Table 7 summarizes the number of identified connector and disconnector across all states in dynamic and static analysis. In static analysis, the disconnection observed within the nodes of the visual, the subcortical, and the cognitive control domains and between nodes of the subcortical, the cognitive control, and the auditory domains. In addition, two abnormal integrations observed between the cognitive control, the visual, and the auditory domains. In dynamic graph analysis within the default mode network, we found four missing links associated with blocked paths in States 1, 4, and 5. Several of the observed disconnectors showed agreement between static and dynamic approach. Four of six disconnections in static analysis, which were observed within and between the subcortical, the cognitive control, and the auditory domain, were also identified in different states of dynamic analysis. In dynamic approach, the disruption identified between the thalamus (SC) and the hippocampus (CC) in State 5; between the insula (CC) and the middle temporal gyrus (AU), and between the supplementary motor area (CC) and the inferior frontal gyrus (CC) in State 2; and between the subthalamus/hypothalamus (SC) and the thalamus (SC) in State 3 were observed in static analysis as well. However, results dynamic approach can capture more detailed information regarding the disconnection and abnormal integration that would be missed in static analysis. Therefore, path analysis on time-varying graphs extends findings from static analysis and provides additional results that are complementary to, and extend, a static analysis.

The default mode network describes a distributed large-scale functional brain network that consists of brain components such as the medial prefrontal cortex, posterior cingulate cortex, and precuneus (Fox, 2007) that is more active during rest or internal cognitive processing compared to external goal-directed cognitive tasks (Raichle et al., 2001). Previous fMRI studies in SZ have consistently reported the disruption and aberrant functioning of the default mode network (Du, Pearlson, et al., 2016; Fornito, Zalesky, Pantelis, & Bullmore, 2012; Mars et al., 2012; Öngür et al., 2010; Pankow et al., 2015). Interestingly, our approach identified three new links associated with additional paths within the default mode network and one new link between the default mode network and cognitive control in precisely the same states as disconnectors were found in dynamic approach, which may be related to a compensatory response in the SZ group that warrants future study.

In addition, disruption in the subcortical domain is consistent with and significantly extends recent reports (Damaraju et al., 2014; Woodward, Karbasforoushan, & Heckers, 2012). In more detail, disconnection between the thalamus and the hippocampus, which were identified by our method in both dynamic and static analysis, has been reported previously in Woodward, Karbasforoushan, and Heckers (2012).

Furthermore, the most significant alteration was observed in cognitive control. Cognitive dysfunction in schizophrenia (McTeague et al., 2017) supports the cognitive control disconnections identified by our method, which were observed both in static and dynamic analysis. In dynamic approach, there were five links associated with blocked paths, of which four were within cognitive control, and one was between the cognitive control and auditory network. Also, three new links were observed associated with additional paths within cognitive control, and one new link was between the default mode network and cognitive control networks. Cognitive control is a cognitive and neural mechanism that contributes to the ability to adapt information processing and regulating behavior according to one’s current goals (Miller, 2000), and according to prior studies, cognitive impairments are the main feature of SZ (Cho, Konecky, & Carter, 2006; Ray, Gold, Silverstein, Barch, & Carter, 2017). Our results thus are consistent with and extend findings from prior studies.

To examine case 3, where two distinct brain components are connected and reachable from each other in both control and patient groups, we applied the covariance decomposition method. The dynamic approach results showed that (only) States 3 and 5 have significant group differences (Jones & West, 2005). Although the paths existed between pairs of brain components, particularly within the subcortical and visual domain in both groups, these paths take unique trajectories, and the decomposition method revealed that those distinct paths have significant weight. The results of both static and dynamic analysis suggest that despite the connection between the nodes in both groups in subcortical domain, schizophrenia have a different type of connection, in terms of unique paths trajectories, and the decomposition method showed that those distinct paths have significant weight between the thalamus and the caudate. In addition, dynamic analysis provided more additional differences, which were not revealed in static analysis in the subcortical domain.

Our work shows that differences between control and patient groups are likely the result of “multilink” disruptions in the paths. The results highlight the importance of studying longer path links, especially in the context of complex mental illness such as schizophrenia. Hence, in comparison with the existing methods, more detailed interpretation regarding the differences between control and SZ patients can be obtained using our proposed method. Future work should focus on expanding the approach and replicating our results in additional datasets.

In this work, we analyzed undirected graphs based on partial correlation. Although the graph-theoretical approach can be used to either assess functional or effective connectivity patterns, most human brain graph investigations have applied undirected networks because of inference constraints of directed graphs (Farahani, Karwowski, & Lighthall, 2019; Karwowski, Vasheghani Farahani, & Lighthall, 2019; Liao, Vasilakos, & He, 2017). Future work could utilize directed network analysis that might provide more information on human brain connectivity.

In sum, we provide an approach using GGM to estimate static and time-varying graphs in control and SZ patient groups on resting-state fMRI data. We comprehensively analyze paths and propose an algorithm to estimate links associated with disconnectivity and abnormal integration in the patient group graphs. Also, in cases in which paths exist between two nodes in two groups, we suggest using a covariance decomposition method. We apply our method to study resting-state fMRI data in SZ versus controls, as SZ is a disorder characterized by disconnectivity. We provided detailed information about how reduced integration or abnormal integration manifests in the SZ group. Several missing links associated with disconnectivity were identified in the SZ group, both within and between functional domains, particularly within the default mode network and cognitive control domains. Also, our proposed algorithm identified additional new edges within these domains associated with abnormal integration, which may be related to a compensatory brain response in the schizophrenia group. In analyzing cases in which paths exist between two specific brain components in both groups, these paths in some cases took unique trajectories, and the decomposition method showed that those distinct paths have significant weight. The proposed path analysis provides a way to characterize individuals by evaluating changes in paths, rather than just focusing on pairwise relationships. As the edges within a path might be distinct in different individuals, a path-based approach can capture important information that is ignored by approaches that focus on pairwise relationships. Our first dynamic path analysis found path-based differences between individuals with schizophrenia and healthy controls, indicating its promise for identifying path-based abnormalities in dynamic connectivity in the healthy and disordered brains. Path analysis on time-varying graphs extends findings from static graphs and provides additional results that are complementary to, and extend, a static analysis.

Data collection was supported by the National Center for Research Resources at the National Institutes of Health (grant numbers: NIH 1 U24 RR021992, NIH 1 U24 RR025736-01).

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00247.

Haleh Falakshahi: Formal analysis; Methodology; Software; Validation; Visualization; Writing – original draft; Writing – review & editing. Hooman Rokham: Conceptualization; Writing – review & editing. Zening Fu: Data curation; Writing – review & editing. Armin Iraji: Writing – review & editing. Daniel H. Mathalon: Resources; Writing – review & editing. Judith M. Ford: Resources; Writing – review & editing. Bryon A. Mueller: Resources; Writing – review & editing. Adrian Preda: Resources; Writing – review & editing. Theo G. M. van Erp: Resources; Writing – review & editing. Jessica A. Turner: Writing – review & editing. Sergey Plis: Writing – review & editing. Vince D. Calhoun: Conceptualization; Funding acquisition; Investigation; Supervision; Validation; Writing – review & editing.

Vince D. Calhoun, National Institutes of Health (https://dx.doi.org/10.13039/100000002), Award ID: R01MH118695. Vince D. Calhoun, National Science Foundation (https://dx.doi.org/10.13039/100000001), Award ID: 2112455.

Brain functional (network) connectivity:

Temporal statistical relationships between brain regions or networks.

Gaussian graphical model (GGM):

An undirected graph of partial correlation that illustrates the conditional independences of multiple random variables. Edges are regarded as the conditional dependencies between the nodes.

Static functional (network) connectivity:

A statistical relationship between brain regions or networks using the entire length of the data.

Intrinsic connectivity network (ICN):

An estimation of a set of voxels obtained using ICA with coherent functional activity over time that can be considered as one functional unit.

Independent component analysis (ICA):

A statistical technique that is used for estimating maximally independent components.

Brain states:

Brain connectivity patterns (states) that reoccur over time and subjects and can be identified by applying k-means clustering algorithm.

Dynamic functional (network) connectivity:

A statistical relationship between brain regions or networks using a portion of the data to study time-resolved changes in connectivity.

Connected component:

A set of nodes in a graph in which there is at least one path between them. That is, it is possible to get from one node to another within a connected component.