A graph comprises of a set of nodes (brain regions) that are interconnected by edges (connectivity weights), and network measures quantify different characteristics of the topology of such graphs. Rubinov and Sporns (2010) illustrate the applicability and interpretation of several complex-network measures in brain imaging. Among them, functional segregation, necessary for optimal specialized processing, informs about dense-connectedness within separate subnetworks. It quantifies whether the regions connected to a given node are connected amongst themselves, thus forming subnetworks wherein majority of the nodes are connected to every other node. For example, during altered consciousness, segregation is reduced, especially in the thalamus (Crone et al., 2013). In the current study, we employed transitivity (global whole-brain-level measure), clustering coefficient and local efficiency (both local node-level measures) to quantify segregation (Rubinov and Sporns, 2010).

In contrast, functional integration captures the ease of interaction between segregated regions (Rubinov and Sporns, 2010). For example, there is elevated segregation in prefrontal and cerebellar subnetworks in attention-deficit hyperactivity disorder (ADHD), but lower integration between these subnetworks (Lin et al., 2014), indicative of characteristic attentional reaction-time deficits observed in this population. In the current study, we employed global efficiency (global measure), shortest path length, and edge betweenness (both local measures at connection-level) to quantify integration (Rubinov and Sporns, 2010). Although traditional connectivity identifies standalone aberrant connections in clinical groups, these measures of integration identify those connections that are not only important by themselves, but are also important for the rest of the connections in the network.

It has been extensively demonstrated that segregation and integration are disrupted in psychiatric disorders [for example, see (Yu et al., 2013; Rocca et al., 2014)]. Most report a narrow, but balanced relationship between them in healthy populations (called metastability) (Hellyer et al., 2015), which is impaired in neurologic and psychiatric disorders (Yu et al., 2013; Rocca et al., 2014). Using resting-state fMRI and our novel framework, we investigated network-level aberrations in soldiers with posttraumatic stress disorder (PTSD) and post-concussion syndrome (PCS) associated with documented mild traumatic brain injury(s) (mTBI). PCS is an outcome of mTBI, in which the individual presents persistent post-concussive symptoms 3 months’ post-injury.

In military populations, there is considerable comorbidity between mTBI and PTSD (Hoge et al., 2008, 2009), often attributed to life-threatening events such as exposure to blast from improvised explosive devices (IEDs), which result in mTBI as well as psychological trauma. With the prevailing clinical approaches focusing on patient reporting, and with substantial overlapping symptoms between PTSD and PCS (Eierud et al., 2014), a better comprehension of the neurobiological-mechanistic basis for PTSD and PCS is imperative for improved diagnosis and treatment outcomes, and for making return-to-duty decisions. Prior fMRI works on comorbid PTSD and mTBI are limited (Spielberg et al., 2015), although its prevalence is considerably high in general society as well as military populations (Veterans, 2015). In the current study, we explored our novel framework involving functional segregation and integration in three groups: soldiers with elevated posttraumatic stress symptoms (PTSD group), PCS + PTSD (comorbid group sustaining both PTSD and PCS), and healthy combat controls. For the sake of disambiguation, we call complex-network modeling as “network-level,” while connectivity modeling is termed “connectivity-level” and activation analysis as “region-level.” Although several studies have identified region-level and connectivity-level aberrations in specific key prefrontal and subcortical regions in mTBI and PTSD (Simmons and Matthews, 2012), a thorough understanding of the aberrations of directional relationships and associated changes in network structure have not emerged from them. We address this limitation in this study.

Graph-theoretic analysis begins from network graphs constructed using pairwise connections, which can be obtained through connectivity modeling. Although functional connectivity (FC) has been the predominant choice so far, we sought to investigate directional networks with causal relationships instead of co-activation (a non-directional entity). It has not been adequately explored, even though it is an equally important mechanism for network-level interactions. Causal connectivity has been discovered even in fMRI timescales (Roebroeck et al., 2005; Abler et al., 2006; David et al., 2008; Deshpande et al., 2011; Deshpande and Hu, 2012; Ryali et al., 2016; Rangaprakash et al., 2018a), indicating that identifying causal networks in addition to co-activation networks is important for a more extensive characterization. Further, PTSD and PCS are typically considered as prefrontal dysregulation disorders (Simmons and Matthews, 2012), meaning that prefrontal causal connectivity is compromised. This provided the impetus for us to further investigate directional connectivity. To our surprise, there have been hardly any fMRI studies investigating effective connectivity (EC) in either PTSD or PCS or the comorbid condition.

Granger causality (GC), an exploratory technique, was employed to quantify EC (Deshpande et al., 2010b). It is the most prevalent technique for deriving causal relationships in natural systems (Kirchgässner et al., 2012). Both recent simulations (Ryali et al., 2011; Wen et al., 2013) and experimental results, including optogenetics and electrophysiology (David et al., 2008; Katwal et al., 2013; Ryali et al., 2016; Wang et al., 2016), demonstrate that GC is reliable for drawing inferences regarding directional relationships between brain regions when used after deconvolving the hemodynamic response function (HRF) from fMRI data (as done in the current study). Several recent fMRI works have also employed this technique (Deshpande et al., 2013; Sathian et al., 2013; Grant et al., 2014; Lacey et al., 2014; Wheelock et al., 2014; Feng et al., 2015; Grant et al., 2015; Hutcheson et al., 2015; Bellucci et al., 2016).

Most studies investigating fMRI connectivity assume connectivity as stationary over time, although static connectivity does not capture dynamic variations of connectivity. While an fMRI scan endures for several minutes, mental processes occur within a few milliseconds to a few seconds’ time, implying that connectivity varies over the timescales of fMRI scans, and that those variations contain biologically relevant information (Hutchison et al., 2013), which are different from that contained in static connectivity (Jia et al., 2014). Recent works have found connectivity dynamics to be a unique and important marker of brain functioning (Hansen et al., 2015; Jin et al., 2017). Therefore, the current study utilized both static EC (SEC) and dynamic EC (DEC). Brain networks were constructed from strength (SEC) and temporal variability (DEC) of directional connectivity, using which we obtained strength and variability of segregation/integration measures, respectively. Such a characterization of dynamic network properties is one of the important novel contributions of this work. While dynamic connectivity has prevailed in neuroimaging for some time (Hutchison et al., 2013), for the first time we introduce dynamic modeling of segregation and integration in a novel framework.

Lower variability of connectivity over time is associated with both psychiatric and neurologic conditions (Garrett et al., 2013; Jia et al., 2014; Miller et al., 2016; Rashid et al., 2016; Rangaprakash et al., 2017a, 2018a), often corresponding to a lack of cognitive flexibility. Compromised behavioral performance is linked with reduced temporal variance of connectivity in both clinical and non-clinical populations (Sakoğlu et al., 2010; Jia et al., 2014; Rangaprakash et al., 2017a, 2018a). Such reduction is linked to impaired ability in dynamically adjusting to changing conditions (thoughts, behaviors, etc.). A healthy biological system is flexible in response to continual momentary changes within the internal and external milieu of the organism. In those terms, temporally “frozen” connectivity and/or complex-network properties point to compromised brain health. Such a characterization has been done in recent connectivity studies (Jia et al., 2014; Rangaprakash et al., 2017a). Higher variability of connectivity is also considered a marker of greater mental flexibility (Zhang et al., 2016).

In this work, we extend these concepts to the reduced temporal variability (or rigidity) in network properties instead of individual connection strengths. We hypothesized that PTSD and mTBI are characterized by altered strength and lower temporal variability of segregation and integration in directional brain networks. We associated the connections exhibiting suppressed network properties with deflation, given that reduced engagement of certain prefrontal-subcortical and prefrontal-cortical pairwise connectivities may be considered as an outcome of impaired regulation from prefrontal regions (Gross, 2014). Similarly, we associated the connections exhibiting elevated network properties with inflation, or pathologically enhanced network-level engagement, given that pairwise hyper-connectivity is seen as an outcome of neurological disruption (Hillary et al., 2015), and has been noticed in PTSD (Hayes et al., 2012; Simmons and Matthews, 2012; Cisler et al., 2014). Within this framework, we sought to identify such networks properties which were (i) affected by PTSD but not mTBI (we call this hypothesis-1, see Figure 1A), and (ii) affected by PTSD as well as comorbid PTSD and mTBI (we call this hypothesis-2; see Figure 1B). Such dichotomy would enable us to identify both common (hypothesis-1) and distinguishing (hypothesis-2) network features between PTSD and mTBI, given the high comorbidity and overlapping symptomatology between them (Spielberg et al., 2015). Notably, we tested the hypothesis on whole-brain data, in a data-driven manner without imposition of any priors, using resting-state fMRI, which is not task dependent. With the network properties that fit our hypothesis, we assessed their association with relevant behaviors (neurocognitive functioning, and symptom severity in PTSD and PCS).

Statistical group separation is the analysis framework for our hypothesis. However, statistical separation does not automatically attribute them with predictive diagnostic ability (Deshpande et al., 2010a) at the individual-subject level. Machine-learning classifiers have been successfully utilized on fMRI data for such diagnostic prediction in disorders like major depression (Deshpande et al., 2009), PTSD (Liu et al., 2015), Parkinson’s (Marquand et al., 2013), dementia (Chen et al., 2011), ADHD (Deshpande et al., 2015), prenatal-cocaine-exposure syndrome (Deshpande et al., 2010a), autism (Deshpande et al., 2013; Libero et al., 2015), and many others. However, to the best of our knowledge, there have been no works utilizing complex-network properties in PTSD/mTBI classification. Given the unique high-level information contained in network properties, we expected network measures to possess predictive ability. Neuropsychiatric conditions such as PCS and PTSD are currently diagnosed through clinical observation and self-report, hence classification using neuroimaging-based network signatures can be useful in obtaining more accurate diagnoses in these highly comorbid conditions. Hence, we employed a machine learning technique to identify highly predictive features by recursively eliminating unimportant complex-network features in a data-driven way. In addition, we sought to find an overlap between connections with network properties satisfying our primary hypotheses (Figure 1), and those possessing high predictive ability. As our secondary hypothesis, we hypothesized that such network properties would predict the diagnostic membership of a new subject better than available non-imaging measures (neurocognitive, behavioral and self-report measures), thus underscoring their relevance to the underlying neuropathology of mTBI and PTSD. We place special emphasis on network properties having all the desirable qualities assessed in this work: high statistical separation, behavioral relevance and high predictive ability. Our study illustrates the utility of our methodological framework using the PTSD/mTBI cohort as an example.

Active-duty soldiers (aged between 18 and 50 years) were recruited from Fort Benning, GA, United States and Fort Rucker, AL, United States to participate voluntarily in the study. The study was conducted in accordance with the Declaration of Helsinki, and the procedures were approved by Auburn University’s Institutional Review Board (IRB) as well as the Headquarters United States Army Medical Research and Materiel Command, IRB (HQ USAMRDC IRB). Written informed consent was obtained from all participants.

Eighty-seven male, active duty U.S. Army soldiers were enrolled in the study, which included 17 with PTSD, 42 with comorbid PCS and PTSD (PCS + PTSD), and 28 combat controls (all groups were matched in age, education and race), all having combat experience in Iraq (Operation Iraqi Freedom, OIF) and/or Afghanistan (Operation Enduring Freedom, OEF). Participants were grouped based on symptom severity in PTSD using the “PTSD Checklist-5” (PCL5) score, post-concussive symptoms using the “Neurobehavioral Symptom Inventory” (NSI) score, clinician referral and medical history. (i) Participants with post-concussive symptoms, clinician referral, history of medically documented mTBI, and scores ≥ 38 on the PCL5 and ≥ 26 on the NSI were grouped as the comorbid PCS + PTSD group. (ii) Participants with no history of mTBI in the last 5 years, a score ≥ 38 on PCL5, and < 26 on NSI and clinician referral were grouped as PTSD. (iii) Participants with score < 38 on PCL5 and < 26 on NSI, no mTBI within the last 5 years, no DSM-IV-TR or DSM-V diagnosis of a psychiatric disorder (based on medical records), and no history of moderate-to-severe TBI were grouped as combat controls. All participants were assessed by a licensed medical practitioner, and reported being deployed to a combat environment. Those with psychotic, mood or substance dependency disorders were excluded.

A battery of psychological health measures were administered to the participants prior to their MRI scan. The battery consisted of the Brief Traumatic Brain Injury Screen [BTBIS; (Schwab et al., 2007)], PCL-5 (Blevins et al., 2015), NSI (Cicerone and Kalmar, 1995), Life Events Checklist (LEC; (Gray et al., 2004)), Combat Exposure Scale [CES; (Guyker et al., 2013)], Childhood Environment [CE; (King et al., 2003)], Zung Depression Scale [ZDS; (Zung et al., 1965)], Zung Anxiety Scale [ZAS; (Zung, 1971)], Alcohol Use Dependency Identification Test [AUDIT; (Saunders et al., 1993)], and the Epworth Sleepiness Scale [ESS; (Johns, 1991)]. In Supplementary Section “Psychological Health Measures,” (SI-2.1) we present, in more detail, the measures that were most relevant for this study [PCL5, NSI, and CNS vital signs (CNS-VS]). Neurocognitive composite index (NCI) was derived from CNS-VS domain scores (Gualtieri and Johnson, 2006) as an aggregate measure of neurocognitive functioning.

For procedures done prior to the fMRI scans during the scheduled appointment, see Supplementary Section “Procedures” (SI-2.2).

Standard resting-state fMRI data pre-processing steps were done including realignment, normalization to MNI space, detrending and regressing out nuisance covariates such as six head-motion parameters, white matter signal and cerebrospinal fluid signal, and band-pass filtering (0.01–0.1 Hz). The largest permitted head motion was half the voxel-size (1.5 mm); no significant group differences were observed in participant head-motion (p > 0.05) (also see Table 1). Pre-processing was performed using Data Processing Assistant for Resting-State fMRI (DPARSF v1.7) (Chao-Gan and Yu-Feng, 2010), which is based on Statistical Parametric Mapping (SPM8) (Friston et al., 2007) and Resting-State fMRI Data Analysis Toolkit (Song et al., 2011).

Deconvolution was then carried out on voxel-level time series, because confounds arising from spatial and inter-subject variability of the hemodynamic response function (HRF) (Handwerker et al., 2004; Rangaprakash et al., 2017c) could lead to a scenario wherein two fMRI time series show high effective connectivity but the underlying neural signals are not highly connected, and vice versa (refer to Figure 2 for an illustration) (Rangaprakash et al., 2018b, c). Such phenomena have been specifically found in the case of PTSD and mTBI with functional connectivity (Rangaprakash et al., 2017c). Additionally, causal connections could potentially switch directions in case the underlying HRFs possess different times-to-peak. In this respect, it has been demonstrated that deconvolution results in improved estimation of effective connectivity (David et al., 2008; Ryali et al., 2012, 2016). The viewpoint of cellular neuroscience on BOLD fMRI presented in a recent paper (Hall et al., 2016) discussed several caveats in the interpretation of fMRI results, wherein careful consideration is warranted based on the underlying cellular mechanisms. Neurovascular dynamics or HRF variability is one such primary issue, about which they comment as follows: “advances in cellular neuroscience demonstrating differences in this neurovascular relationship in different brain regions, conditions or pathologies are often not accounted for when interpreting BOLD.” They advise employing computational modeling (e.g., deconvolution) to mitigate the issue. We employed a popular blind deconvolution algorithm (Wu et al., 2013). Many recent papers have employed it [see for example (Amico et al., 2014; Lamichhane et al., 2014; Boly et al., 2015; Rangaprakash et al., 2017a, b)]. The deconvolution is blind since both the HRF and underlying latent neural time series are estimated only from the recorded fMRI data. Resting-state fMRI data is modeled as event-related using point processes with randomly occurring events; then, voxel-specific HRFs are estimated using Wiener deconvolution. This technique is date-driven; hence, we do not encounter overfitting issues that often plague model-based approaches.

Since whole-brain fMRI data has high dimensionality, 125 functionally homogeneous brain regions spread out across the cerebral cortex and encompassing it completely, determined using spectral clustering [known as the Craddock-200 atlas (Craddock et al., 2012)], were taken and mean deconvolved time series were obtained from them. All further analyses (carried out on the Matlab^® platform) utilized these 125 time series from every participant.

A precursor to obtaining network-level characterization is to first get the connectivity network itself, which is, computing SEC and DEC matrices from pre-processed fMRI data. Whole-brain SEC and DEC were computed using GC (Deshpande et al., 2010b), which is an exploratory technique used to quantitatively measure directional relationships between brain regions. While SEC uses a multivariate autoregressive (MVAR) model, DEC employs a dynamic MVAR model evaluated in a Kalman filter framework using variable parameter regression (Wheelock et al., 2014).

The concept of Granger causality (GC) is that, if future values of time series “Y2” can, in a mathematical sense, be predicted by the past values of time series “Y1,” then a causal influence is inferred from time series Y1 to time series Y2 (Granger, 1969). GC’s MVAR model predicts one time series from the other quantitatively, as described briefly next. Given k different time series Z(t) = [z₁(t), z₂(t), … z_k(t)], wherein k corresponds to 125 ROIs of this work, the MVAR model of order p is given by:

Here, E(t) is the model error, while M(0) … M(p) are model coefficients. Like in earlier studies (Deshpande et al., 2010b), this formulation included a zero-lag term with coefficient M(0) which would eliminate the contribution of zero-lag cross-correlation between the time series. Since M(0) represents co-variance between time series and not used in GC computation, the effect of zero-lag cross-correlation gets ignored in GC. Given that diagonal elements of M(0) are set to zero, we only model the instantaneous cross-correlation, and not the auto-correlation between the time series.

The coefficients were estimated using multivariate least-squares estimation. It computes the set of optimal coefficients with model error being minimized in a least-squares sense. The model order p can either be chosen by utilizing a mathematical principle such as the Bayesian Information Criterion (BIC) (Roebroeck et al., 2005) or based upon the needs of the application being considered. In neuroimaging, causal relationships corresponding to neural delays less than or equal to the TR are of interest (Deshpande et al., 2013), hence we used a first order model. Given that fMRI’s temporal resolution is relatively low, a first-order model has been shown to capture the most relevant directional information (Deshpande and Hu, 2012).

The degree to which the past Z(t-p) is able to predict the present Z(t) is given by the coefficient matrix M(p). The sum of all such coefficients would then correspond to the degree to which the past values put together can predict the present. As in prior works (Kaminski et al., 2001), GC was formally derived, based on the model coefficients, as:

(2) G ⁢ C i ⁢ j = ∑ n = 1 p m i ⁢ j ⁢ ( n )

Wherein m_ij are the elements of matrix M and GC_ij refers to the SEC value from ROI i to ROI j. Notably, a single coefficient matrix was computed for the entire duration of the experiment, that is, coefficients are not varying with time. A deeper theoretical rendering of GC can be found here (Deshpande et al., 2010b). GC-based methods have been experimentally validated for fMRI EC analysis (David et al., 2008; Katwal et al., 2013), and they have been extensively utilized for fMRI EC modeling in recent times [see for example (Deshpande et al., 2011, 2015)].

Next, DEC was computed by employing time-varying dynamic Granger causality (DGC), evaluated using a Kalman filter framework. Dynamic multivariate vector autoregressive (dMVAR) model was employed for estimating DGC (Grant et al., 2014; Wheelock et al., 2014). DEC is the underlying time-varying physiological process, while DGC is the mathematical measure that quantifies it. This technique has also been used in several recent studies (Deshpande et al., 2013; Wheelock et al., 2014; Hutcheson et al., 2015). Unlike GC formulation, dMVAR model coefficients M’(p,t) are a function of time, hence the model is “dynamic.”

A Kalman filter framework which uses variable parameter regression (Büchel and Friston, 1998) was used to estimate dynamic model coefficients, which involved imposing a forgetting factor (which was chosen as 1 in our case). DGC was thus computed as:

(4) D ⁢ G ⁢ C i ⁢ j ⁢ ( t ) = ∑ n = 1 p m i ⁢ j ⁢ ( n , t )

Where m_ij are the elements of matrix M and DGC_ij(t) is the value of DEC from ROI i to ROI j at a given time point t. Like in GC, zero-lag cross-correlation effects were compensated here also. Further, a forgetting factor of 1 was used to make the system well-conditioned so that the coefficients may be estimated accurately.

A 125 × 125 whole-brain SEC matrix was obtained for every participant by computing GC between all combinations of connections between the 125 ROIs. With DEC, the dynamic MVAR model coefficients are a function of time, hence, with our fMRI data having 1000 time points, we obtained a 125 × 125 × 1000 DEC matrix per participant. SEC and DEC matrices were used for further complex-network analysis. To illustrate the concepts of SEC and DEC in the context of neuroimaging, we show a simple illustration using a pair of fMRI time series from our experimental data (see Figure 3).

We first describe the network measures of segregation and integration, and then explain how they were used in the context of this work. As noted earlier, given the complexity of our hypothesis, we dealt with weighted directed networks in this work. Functional segregation was quantified using transitivity (global measure, one value for whole brain per participant), clustering coefficient and local efficiency (both local measures, one value per node/region per participant). Functional integration was quantified using global efficiency (global measure), shortest path length and edge betweenness (both local measures, one value per connection per participant). We obtained source codes for these measures from the Brain Connectivity Toolbox (April 2014 release) (Rubinov and Sporns, 2010), and implemented the entire pipeline in the Matlab^® platform through custom codes. A detailed account of these measures can be found in Rubinov and Sporns (2010). For the benefit of readers, we have explained each of these measures in detail using a simplified example network in Supplementary Section “Complex-Network Analysis” (SI-2.3).

Briefly, transitivity is a global measure of overall efficiency of local processing in the brain. Clustering coefficient (CC) gives a transitivity-type characterization for every node. Local efficiency (EffLoc) is closely related to CC, wherein nodes with powerful neighbors that are involved in several shortest paths have higher EffLoc, indicating that the node is important in the sub-network for specialized processing. While CC and EffLoc usually give similar (but not same) results, their interpretations are different. In this work, along with transitivity as the global measure, we employed both CC and EffLoc as local measures, which are the two popularly used local measures of segregation. We took an overlap (intersection) of the final significant group differences for the two measures, so that the affected nodes had differences in both the measures, thus providing more conservative results with a broader interpretation.

Global efficiency (EffGlob) is a global measure indicating the aggregate ease of communication in the entire network. Shortest path length (SPL) is a measure of how easy it is to reach one node from the other, and is analogous to meta-connectivity. Edge betweenness (EB) measures the number of all shortest paths in the entire network that contain the given connection. Like SPL, EB is an important network measure because it characterizes the importance of a connectivity path not only through its pairwise connectivity value but also through the significance of the connectivity path for other connectivity paths present in the network. If a connectivity path matters a lot for other paths, i.e., for communication between various other nodes, then the given path would have high integration ability (i.e., SPL and EB). Such a characterization can be obtained only through complex-network modeling since traditional pairwise connectivity informs us only about the strength of interaction between just two brain regions. In this work, we employed EffGlob as the global measure and both SPL and EB as local measures of integration. As with segregation, we took an overlap (intersection) of the final significant group differences in the two local measures, so that the affected paths had differences in both measures, thus providing more conservative, but potentially more reliable, results with wider interpretation.

Next, we describe how these six network measures were used in the context of this work. SEC and DEC connectivities were used separately to construct static and time-varying networks with brain regions as nodes and connectivity strengths between them as the weighted directed edges of the network graphs. Absolute value of connectivity was used to construct the network graphs. With SEC, a single network was constructed for the entire duration of time in the data, giving a “connectivity strength” network, which was used to obtain each one of the six complex-network measures for every run of every participant. With DEC, we considered each time point of the DEC time series as a snapshot of the network at that time instant, and then constructed a graph using the nodes and edge values from that snapshot. We computed network measures for that snapshot, and repeated the procedure independently for the rest of DEC time series to obtain a time series of values for each network measure. Then, for each network measure, we computed the variance of the network measure time series to obtain a single value for the entire duration of the data. This gave us a network with paths whose weights corresponded to the temporal variability of complex-network measures. This was obtained for every measure and for every run of every participant, similar to SEC.

Statistically significant differences in these strength and variability networks were obtained, in accordance with our hypothesis (p < 0.05, FDR corrected). We corrected for 31250 comparisons: 125 comparisons of segregation (125 ROIs), 15500 comparisons of integration (125 × 125 − 125), each for both static connectivity and variance of dynamic connectivity networks. Differences were controlled for age, race, education, and head-motion [using mean frame-wise displacement, as defined by Power et al. (2012)]. That is, we found significant group differences with both SEC and DEC derived complex-network measures separately for these three pairwise comparisons (thus giving a total of six comparisons per network measure): Control vs PTSD, Control vs. PCS + PTSD, PTSD vs. PCS + PTSD. We then identified the common network measures among four of these comparisons (hypothesis-1) which excluded PTSD vs PCS + PTSD comparison, and we also identified common network measures among all the six comparisons (intersection, hypothesis-2), all of which also fit our hypothesis, that is, conformed to the increasing/decreasing trend as we moved from Control to PTSD to PCS + PTSD.

It is notable that we have taken a conservative approach in this work. We opted to look for common differences in pairwise statistical comparisons, rather than performing a single three-way statistical comparison, which is less conservative. We obtained common differences in static as well as dynamic network measures, and we also constrained the differences to conform to a trend as per our hypothesis. Additionally, we computed two local measures in segregation as well as integration, and considered only common differences in them, which added another level of constraints on our findings. In addition to these, we notably discarded any paths which had significant network-level differences in local measures of integration (i.e., SPL and EB), but not significant pairwise effective connectivity differences themselves. That is, we included only those paths which had significantly different SEC and variance of DEC in accordance with the trend set out in our hypothesis (p < 0.05, FDR corrected, controlled for age, race, education and head-motion), in addition to having significantly different local measures of integration (i.e., SPL and EB). This was done to ensure that, irrespective of network-level disturbance, the significant connections that emerged in this work would have also cleared whole-brain multiple-comparisons-corrected statistical threshold with traditional static and dynamic effective connectivity like in most other studies. This reassured that our results conformed to multiple layers of validation, verification and statistical standards, and that evidence of network disruption were obtained via multiple analysis approaches, in addition to providing novel insights through network characterization.

In an effort to assess the behavioral relevance of complex-network measures, we first obtained the association of the strength and variability of complex-network measures (only those which fit our hypothesis) with symptom severity in PCS (NSI score) and PTSD (PCL5 score), as well as neurocognitive functioning (NCI score and subtests). Neurocognitive functioning (e.g., executive functioning, cognitive flexibility) is often impaired in psychiatric disorders such as PTSD and PCS (Simmons and Matthews, 2012), hence identifying such network properties associated with it would be important. We report significant associations between complex-network properties and behavioral/clinical measures.

In order to obtain additional insight into how network properties of the ensemble of identified connections mapped on to the ensemble of behaviors, we performed partial least squares regression (PLSR) analysis (Krishnan et al., 2011), which we employed to predict neurocognitive functioning (NCI and subtests) and symptom severity (PCL5, NSI) from strength and variability of network measures obtained from our prior analysis. We present the percentage variance in behaviors explained by the complex-network measures.

For predicting the diagnostic membership of a novel subject based on a novel measurement using the measure, success in hypothesis testing is neither necessary nor sufficient. A mechanism to quantify the predictive ability of the features is not available with the hypothesis testing framework, requiring us to acknowledge what a technique like hypothesis testing can do, and cannot do. Statistically significant network properties necessarily need not have high predictive ability, and vice versa. Hence, those network properties that are both statistically significant (in accordance with our hypothesis) and are top-classifiers (high predictive ability) carry superior importance and relevance. Therefore, we used machine learning techniques to identify such network properties (features) which can accurately classify individuals between controls, PTSD, and PCS + PTSD. A Recursive Cluster Elimination based Support Vector Machine (RCE-SVM) classifier (Deshpande et al., 2010a) was used to classify the participants based on whole-brain network properties (both strength and variability). Notably, findings from prior complex-network analysis were not used to bias the machine learning analysis as whole-brain data was used. A detailed account of this technique can be found in Rangaprakash et al. (2017a), and we have explained it thoroughly in Supplementary Section “Machine Learning Classification Analysis” (SI-2.4) to benefit the reader.

Briefly, RCE-SVM iteratively eliminates features to minimize prediction error. The training data is clustered, and upon SVM classification the clusters are scored using testing data. Low-scoring clusters are eliminated (RCE step) and the procedure is repeated until only the top-predictive features remain. In this work, we made the following parameter choices. The training set consisted of 80% of the participants, while the testing set consisted of the remaining 20%. We began the algorithm with forty clusters in the first RCE step. Based on performance, the bottom 20% of the clusters were eliminated in every subsequent RCE step. Two clusters containing the top-predictive features remained in the final RCE step. With a hundred random iterations, sixfold cross validation was performed in every iteration, resulting in a total of 600 iterations over the complete execution.

To be conservative, we obtained the worst-case classification accuracy by evaluating the lowest accuracy value gathered from test data among all 600 iterations (sixfolds × 100 repetitions). Statistical significance of the accuracies was computed through estimating p-values using a binomial null distribution B(η,ρ), with ρ being the probability of accurate classification and η being the number of participants like in previous studies (Pereira et al., 2009). Only accuracies with p < 0.05 (Bonferroni corrected) were taken as statistically significant.

We repeated this procedure and performed classification independently with 32 available non-imaging measures as input features instead of network measures. The 32 measures were: (i) psychological health measures: Perceived Stress Scale, Epworth Sleepiness Scale, Pittsburgh Sleep Quality Index, Zung Depression Scale, and Zung Anxiety Scale; (ii) behavioral measures: all CNS-VS measures including the NCI score; (iii) exposure/injury descriptives: CES, lifetime concussions, and Life Events Checklist. Worst-case accuracies and top-classifying features were obtained, with them being compared with the results obtained by using complex-network measures.

Network properties having statistical significance in accordance with our hypothesis, having behavioral relevance as well as having high predictive ability were attributed distinctive importance in this study. Using such features, we finally performed support vector regression (SVR) to predict PTSD and PCS symptom severity, in order to assess those features dimensionally. Similar to the classification analysis, we performed sixfold cross-validated linear SVR over one million iterations. Specifically, in each iteration, the regression model was developed using 5/6th of the randomly chosen participants. The model used features described above as inputs and learned the underlying function which maps onto the PCL5 and NSI scores. Subsequently, the model was used to predict PCL5 and NSI scores in the remaining 1/6th participants. Our machine learning classification and regression analyses involved no hyperparameter optimization. We report correlation (R²) between predicted and measured symptom severity scores.

Figure 4 summarizes the processing pipeline of all the methods.

The demographics (for the three groups) are presented in Table 1. There were no significant group differences in age, p = 0.70, or education, p = 0.15. The results indicated a difference in the frequency of reported psychotropic use between the groups, τb = 0.24, p = 0.01, with the highest percentage of medicated participants being in the comorbid group. The number of reported lifetime mTBIs also had significant group differences specifically between control group and the PCS + PTSD group [F(2,171) = 5.81, p = 0.004], but not the control versus PTSD groups or PTSD versus PCS + PTSD groups, p > 0.05.

The results revealed significant differences between the three groups in posttraumatic symptoms (PCL5), [F(2,81) = 101.65, p < 0.001], post-concussive symptoms (NSI), [F(2,78) = 49.79, p < 0.001], and CES, F(2,79) = 40.69, p < 0.001. All p-values remained significant after corrections for multiple comparisons. As observed in Table 1, the PCS + PTSD group had the highest scores out of the three groups on these respective measures.

The results indicated that, after corrections for multiple comparisons, the control group displayed significantly higher scores than the PCS + PTSD group on all neurocognitive measures, p < 0.05, except for reaction time and verbal memory, p > 0.05. The PCS + PTSD group also had significantly lower scores in executive functioning, cognitive flexibility, and the NCI compared to the PTSD group, p < 0.05. The findings suggest that both the PTSD and PCS + PTSD groups display lower scores than controls, but also, the comorbid group had greater impairments than the PTSD group (see Table 1).

We used SEC and DEC connectivity matrices to compute six complex-network measures (two global and four local measures). With global measures (Figure 5), we found significantly lower strength and variability of both segregation and integration in PTSD and PCS + PTSD compared to controls. Our finding indicates that both specialized processing and efficient communication are compromised in the disorders at the whole brain level. However, no significant differences were found between PTSD and PCS + PTSD groups, indicating that PTSD might contribute to global aberrations whereas the effect of mTBI might be more localized.

Strength and temporal variability of functional integration values of four paths, which were significantly different between all three groups (the yellow connections in Figure 8), as well as the strength and temporal variability of functional segregation of MFG and Insula (Figure 6A) showed significant associations (p < 0.05 Bonferroni corrected) with neurocognitive functioning (NCI) and severity of both PTSD symptoms (PCL5-score) and post-concussive symptoms (NSI-score), thus highlighting their relevance to the underlying pathophysiology (see Table 2). Notably the associations followed the expected trend: increase in symptom severity and decrease in behavioral performance corresponded to higher strength of integration in inflated paths and lower in deflated paths, and lower variability (i.e., rigidity) in integration in all paths (similarly with segregation). However, those connections which were not different between PTSD and PCS + PTSD (green paths in Figure 8), as well as other nodes in Figure 6A and global complex-network measures had no significant associations with symptoms and neurocognitive performance.

Since multitude of network paths and nodes had relevant associations with multiple measures of symptoms and neurocognitive performance (which we will now refer to as neurobehavioral indices), it would be interesting to measure how much variance in the neurobehavioral indices could be explained by those set of network measures. This was accomplished using PLSR (Krishnan et al., 2011), which finds the combined ability of the strength and variability of functional integration of the four connections and functional segregation of two nodes to predict neurobehavior. We found that the strength of network measures could explain 48.95% variance in the neurobehavioral indices, while the temporal variability of network measures could explain 57.17% variance. When both were combined, they could explain 61.74% variance in the neurobehavioral indices. A significantly large association between these network measures and neurobehavior (R² = 0.56, R = 0.75, p = 3.5 × 10^–32) was observed in the latent space (see Figure 9A for linear fit). The latent space consists of categorical variables that represent all network measures and all neurobehavioral indices included into the model, so that their relationship in the latent space could be considered the effective association of all the included network measures with all the neurobehaviors. As such, the latent space variables contain more “information” in them than the individual variables themselves, consequently explaining more variance than individual measures. For this reason, our finding of higher R²-value must not be surprising (Vul et al., 2009), and this is fundamental to the multivariate nature of PLS, as elaborated by Krishnan et al. (2011). Our finding reiterates that the strength and variability of functional integration of the four paths and that of segregation of the two nodes identified in this work are behaviorally relevant.

One a side note, head motion [mean frame-wise displacement (Power et al., 2012)] was not significantly correlated with behavioral measures in latent space (R = 0.049, p = 0.52), complex network measures in latent space (R = 0.056, p = 0.46), PCL5 (R = 0.047, p = 0.54) or NSI (R = 0.015, p = 0.84) symptom severity scores. This enhanced our confidence in the results.

Top predictors are those that, among all network measures, possess the highest ability in predicting the diagnostic membership of a novel subject. To identify the top-predictors, we performed RCE-SVM classification (Deshpande et al., 2010a). Classification was done with two different paradigms: (i) classification using the 32 non-imaging measures (NIMs), and (ii) classification using strength and temporal variability of network measures taken from the entire brain (all data, nothing left out). Results showed that classification using network measures provided significantly better accuracy (approximately 10% more, p < 0.05 Bonferroni-corrected) than classification using NIMs (Figure 9B). This result indicates that network measures have superior predictive ability in identifying individuals with PCS and PTSD as compared to NIMs.

Table 3 shows the worst-case accuracies and top predictive features (for average accuracy, please see Supplementary Section “Supplemental Machine Learning Classification Results” (SI-3.4) and Supplementary Figure S3). Also of considerable interest are the top-predictors that resulted in highest classification accuracy. For classification using network measures, strength, and temporal variability of functional integration of the following four paths were the top predictive features: L_MFG → L_Insula, L_Insula → L_Amygdala, L_Amygdala → R_Hippocampus and R_Hippocampus → L_Precuneus). Coincidentally these four paths also showed statistically significant differences in static as well as time-varying network properties (the yellow paths in Figure 8, which were significantly different between all three groups). Also, coincidentally, these were the same four paths whose network measures had significant associations with neurocognitive functioning and symptom severity. To expand upon this, our findings revealed that, in addition to behavioral relevance and statistical separation, these paths also possessed the highest predictive ability, all obtained in a data-driven fashion from whole-brain complex-network data.

Finally, using these network properties (the yellow paths in Figure 8), we performed SVR to predict PCL5 (and NSI) scores. Over one million iterations, we found that predicted and measured PCL5 scores were significantly correlated (R² = 0.72 ± 0.05, R = 0.85 ± 0.03, P = 3 × 10^–7 ± 10^–6), as also were predicted and measured NSI scores (R² = 0.66 ± 0.04, R = 0.81 ± 0.03, P = 7 × 10^–7 ± 10^–6). With such high predictive ability (i.e., 66–72% variance explained in symptom severity), these network properties assume considerable importance in the context of PTSD and PCS + PTSD pathology. Figure 4 summarizes the processing pipeline of our entire work, along with corresponding results.