ABSTRACT
Network diffusion models are a common and powerful way to study the propagation of information through a complex system and they offer straightforward approaches for studying multimodal brain network data. We developed an analytic framework to identify brain subnetworks with perturbed information diffusion capacity using the structural basis that best maps to resting state functional connectivity and applied it towards a heterogenous dataset of internalizing psychopathologies (IPs), a set of psychiatric conditions, including depression and anxiety, in which similar brain network deficits are found across the swath of the disorders, but a unifying neuropathological substrate for transdiagnostic symptom expression is currently unknown. This research provides preliminary evidence of a transdiagnostic brain subnetwork deficit characterized by information diffusion impairment of the right area 8BM, a key brain region involved in organizing a broad spectrum of cognitive tasks, that may underlie previously reported dysfunction of multiple brain circuits in the IPs. We also demonstrate that models of neuromodulation involving targeting this brain region normalize IP diffusion dynamics towards those of healthy controls. These analyses provide a framework for multimodal methods that identify both brain subnetworks with disrupted information diffusion and potential targets of these subnetworks for therapeutic neuromodulatory intervention based on previously well-characterized methodology.
INTRODUCTION
The dynamics arising from the interaction of individual elements of a complex system are commonly investigated by representing such a system as a network or graph [13]. This has been applied extensively towards studying neural connectivity, where brain regions (represented by nodes) and the links between them (represented by edges) define a brain network, or connectome [71,72]. As in any network, the underlying structure of connectivity constrains the functions and processes that take place within it [78]. Indeed, the mechanisms of propagation of and neural response to brain pathologies have been closely tied to network topology [31].
A number of studies have found that the basis formed by the anatomic structure of the brain and the patterns of observed spatiotemporal neural activity are intimately related [1, 2, 5, 23]. Many of these investigations use a network diffusion model to study the propagation of neural impulses, or information (given by functional connectivity), throughout the structure formed by the white matter tracts of the brain. These so-called network diffusion-based approaches all utilize the properties of the structural graph Laplacian, which encodes how a diffusive process spreads throughout a network over time [20]. The basis of diffusion given by the graph Laplacian, or its eigenmodes, can thus be used to model the flow of neural information [82], and has been shown to play important roles in healthy brain network organization [78]. More specifically, such an approach was used successfully to demonstrate that rs-fMRI functional connectivity is predictable from DTI structural eigenmodes [1, 2]. Furthermore, resting state networks were found to closely match spatial patterns of structural eigenmodes [5]. Of note, many other approaches for investigating this interplay have been studied, including models based on epidemic-spreading [73], threshold models [55] and neural mass models [40] (for a review, see [6]).
Network diffusion-based methods have also been applied towards investigating the patho-physiology of neurological diseases. For example, Raj and colleagues [65] show that white matter structural eigenmodes closely resemble known patterns of dementia and correlate strongly with regional atrophy. A closely related methodology to that of studying structural eigenmodes is using the graph Laplacian-derived heat kernel of brain networks, a matrix which encodes how much information is transferred between every pair of network nodes after a given time, or diffusion depth. Heat kernel methods have been used to characterize perturbations in brain network information transfer in autism [68] and to predict future adverse motor function resulting from premature birth using white matter structural connectomes [19].
In the context of psychiatric conditions, the application of methods that integrate multimodal perspectives in their analytical approaches has great potential to elucidate the distributed perturbation of neurocircuitry that underlie the complex cognitive and behavioral disruption found in patients suffering from these disorders. This may hold particularly true for the internalizing psychopathologies (IPs), including mood (e.g., MDD, dysthymia) and anxiety (e.g., panic disorder (PD), social anxiety disorder (SAD), generalized anxiety disorder (GAD)) disorders. Many previous neuroimaging studies have been conducted on IPs, however, most of which use unimodal data analysis. Importantly, findings about potential neuropathological substrates are more often overlapping between the IPs than distinct to a specific IP. In addition, IPs, as traditionally categorized, are often comorbid with one another and present heterogeneously with a spectrum of related symptoms and disruptions to emotion regulation and negative valence system (NVS) processes [17, 44, 59]. Furthermore, the available first line treatments for the IPs, either cognitive behavioral therapy (CBT) or selective serotonin reuptake inhibitors (SSRI), are equally effective across the swath of the disorders [26]. Many studies have concluded that similar structural and functional brain networks involving regions commonly implicated in the expression of fear, anxiety, negative affect and other NVS features are dysfunctional in these disorders [29, 38, 50]. In addition, previous findings also implicate the disruption of similar canonical resting state networks (RSNs), such as the default mode network (DMN), in both depression and anxiety using both structural and functional neuroimaging [45, 52, 60, 70, 76]. To address this pattern of findings, the National Institute of Mental Health’s Research Domain Criteria (RDoC) initiative was developed [22, 41] in order to reorient the study of psychiatric disorders away from traditional diagnostic categories and towards that of data-driven approaches such as identifying transdiagnostic cognitive domain disruptions, biomarkers of treatment response and targets for neuromodulatory intervention.
In line with the aims set forth by this initiative, we study a dataset from an RDoC clinical trial consisting of DTI and rs-fMRI neuroimaging scans of a treatment naive, heterogeneous IP patient (PT) cohort and age and sex matched healthy controls (HC). PTs were then randomized to either 12 weeks of SSRI or CBT and completed Inventory of Depression and Anxiety Symptoms (IDAS-II) [79] self-reports pre- (Pre) and post-treatment (Post) to assess transdiagnostic dimensions of symptom burden. Several previous analyses have been conducted on this dataset but they are unimodal and have a hypothesis-driven focus on a priori brain region, often studied in association with specific task-based measures of NVS subsystems [14, 15, 32, 36, 48, 64, 77].
Because IPs share similarly dysfunctional brain networks and respond to similar treatments, we hypothesize that there exist pathophysiological features common to all IPs that result in impaired emotion regulation and the subsequent heterogeneous expression of IP symptoms. In this paper we ask (1) can the presence of such unifying perturbations can be identified by incorporating both structural and functional connectivity data, (2) do these perturbations predict response to CBT and/or SSRI treatment and (3) can we identify candidate brain region as targets of neuromodulatory intervention to normalize connectivity dynamics to those found in healthy controls? To answer these questions, we use a multimodal data-driven analysis based on network diffusion models. In this approach, the diffusion basis of the structural connectome, given by the eigenmodes (eigenvector-eigenvalue pairs) of the normalized graph Laplacian, is central to the methodology. The eigenmode basis of a network encodes its information diffusion properties, and the eigenmodes of a graph, scaled exponentially by the diffusion depth or time scale allowed for diffusion to occur, represent an embedding of nodes such that their pairwise Euclidean distance in the embedding is inversely proportional to the ability for information to diffuse between nodes at the given diffusion depth [19, 82, 83]. Next, the mapping between each subject’s structural and functional network that minimizes the error between the empirical and estimated functional connectivity is computed by optimizing over diffusion depth parameters as described by Abdelnour and colleagues [1]. We then define a novel representation of structural connectivity, the structural diffusion distance (SDD) connectome, where the edge weights between each brain network node are given by the pairwise Euclidean distance from the diffusion embedding at the optimal time scale that best maps to empirical rs-fMRI activity. We chose such a model for multimodal incorporation as it has both been well characterized in the study of structural to functional mappings in brain networks [1] and allows for the computation of the diffusion-based network embedding given by the structural eigenbasis that most accurately underlies empirical functional connectivity. Thus, we achieve correspondence to functional connectivity in SDD connectomes by computing the embedding at the optimal diffusion depth, βt, for each subject (see Methods section for details). Furthermore, a network diffusion model easily allows for the simulation of information flowing through the network. Simple mathematical properties of the graph Laplacian-derived heat kernel corresponding to the network embedding can be exploited to identify network nodes to supplement with additional ‘heat’ supplement in order to normalize deficient heat transfer throughout a PT subnetwork towards the desired diffusion dynamics of a target (mean HC) subnetwork. Our study applies this analytic framework to provide evidence for transdiagnostic IP subnetwork disruption and neuromodulation treatment targets, opening up application of multimodal network diffusion-based methods towards other brain disorders.
RESULTS
Data analyzed in this study are taken from a previously conducted Research Domain Criteria (RDoC) clinical trial on predictors of IP treatment outcomes (ClinicalTrials.gov identifier: NCT01903447) in which PT were randomized to either 12 weeks cognitive behavioral therapy (CBT) or selective serotonin reuptake inhibitor (SSRI) treatment. In addition, the patient population has heterogeneous clinical presentations and IP comorbidity in order to study common pathological features without regard to traditional diagnostic categories. Baseline rs-fMRI and DTI scans as well as IDAS-II Panic and Depression scores were assessed from both HC (n = 22) and PT (n = 65) subjects. The subset of PT (n = 50 total; n = 28 CBT cohort; n = 22 SSRI cohort) completed the IDAS-II self-reports following treatment are used for identifying correlates of treatment response. Of note, there were no group differences in age or sex between HC and PT. As expected, baseline IDAS-II Panic and Depression scores were significantly different between groups at baseline (table 1).
Figure 1 provides a brief visual overview of the analytic framework for computing the diffusion-based embeddings and subsequent structural diffusion distance (SDD) connectomes used in this study, discussed in detail in the Methods section. We first compute the structure to function mapping as described in [1] for each subject. The fit of the empirical rs-fMRI connectome to the rs-fMRI connectome estimated from the DTI structural connectome Laplacian eigenmodes is then determined by computing the Pearson r statistic between these adjacency matrices. The mean fit values for HC were slightly higher than PT (r = 0.26 ± 0.03 and r = 0.25 ± 0.03, respectively, table 1) but the difference between their mean values was not statistically significance via t-test (p = 0.06), indicating that the relation between resting state functional connectivity and the structural diffusion basis is likely preserved at the global scale. Finally, the distance matrix connectomes used for subnetwork discovery are computed from the structural diffusion basis at the diffusion depth, βt, that optimally maps structural to functional connectivity.
Identification of Diffusion Impaired Subnetwork
After SDD connectomes were computed, the Network-Based Statistic (NBS) algorithm [86] was used to discover subnetworks with altered diffusion properties in PT relative to HC. NBS was applied with one-sided t-tests used for edge comparison, over a range thresholds t = {3.0, 3.5, 4.0, 4.5}. For a t-test contrast (left-sided t-test) of HC < PT, the range of subnetwork sizes varied greatly as a function of t-statistic threshold, resulting in subnetworks from 314 to 2 nodes with thresholds t = 3 and t = 4.5, respectively. With the threshold of t = 4 a single significant subnetwork (SN1) consisting of 48 nodes and 73 edges was identified (Figure 2A) and was chosen for subsequent analyses because of the biological interpretability of its size, and for adequate estimated statistical power, 1 − β = 0.81, with a t-statistic threshold of t = 4. Note that because NBS was conducted on SDD connectomes, edges are pairwise distance between nodes in the diffusion-based network embedding. This indicates that SN1 represents a subnetwork whose nodes are significantly further apart in this diffusion space, i.e., diffusion of information occurs more slowly between the nodes of SN1 along the edges of SN1 in PT relative to HC. For NBS with a t-test contrast defined as HC > PT, no significant subnetworks were found.
To determine the brain regions that are most central to the diffusion impairment in SN1, we calculated the mean strength, strSDD, of each node using the SDD values of edges in SN1 (see Methods section for details). The nodes with the greatest strSDD values are from brain regions found in the bilateral anterior cingulate (ACC) and medial prefrontal (mPFC), right inferior parietal (IPC) and left insular (INS) and frontal opercular (FOP) cortices (Figure 2A). Specifically, the brain region with the greatest average strSDD value in both HC and PT is the right area 8BM, which has recently been reported to be a core region of a ‘multiple demand’ (MD) subnetwork that is active during a broad spectrum of tasks and may play an important role in general cognitive control [4]. Additionally, we found that the hubness of area 8BM was not present in structural connectomes, determined by computing graph theoretical metrics, indicating that the centrality of this brain region is unique to the SDD representations of brain networks (Figure S1).
Baseline Heat Kernel Correlates of Symptom Severity and Treatment Response
To investigate the association of information diffusion between brain regions of SN1 and clinical scales, we compute heat kernel matrices which encode the amount of possible heat transferred between each node after a specified time scale, given by each subject’s diffusion depth, βt. Each value, i, j, of a heat kernel matrix then describes pairwise diffusion of information from the ith to the jth brain network node at the time scale best associated with functional connectivity. We took a hypothesis-driven approach by focusing our correlation analyses on the heat transfer between the bilateral ACC and all other SN1 cortical regions, as the nodes in these cortical regions are central to the diffusion impairment of SN1. The heat kernel values (HKVs) are averaged by cortical region and are then used for computing Spearman correlations with IDAS-II Depression and Panic subscales, which map to distress and fear domains, respectively, based on a previously used approach to broadly divide and assess symptom domains of IPs (Figure 3) [63, 64]. In addition, we use only PT (n = 65) for these analyses, as mean IDAS-II subscale values expectedly differ highly significantly between HC and PT, and correlations would likely be due to group differences. Finally, we report at most the top three most significant correlations that survive FDR correction (n = 44 comparisons) in this section of the manuscript, but all significant correlations are available in Table S1.
First, we determine the HKV correlates of baseline IDAS-II subscales (Figure 3A). Heat transfer between the left ACC/mPFC and INS/FOP (ρ = −0.37, p = 0.003, q = 0.04) and the right ACC/mPFC and left inferior frontal cortex (IFC) (ρ = −0.36, p = 0.004, q = 0.04) was negatively correlated with IDAS-II Depression scales, indicating that less information diffusion between these brain regions is associated with higher IDAS-II Depression values. A positive correlation was found with heat transfer between the left ACC/mPFC and posterior cingulate cortex (PCC) (ρ = 0.39, p = 0.001, q = 0.04). No significant correlations were found using IDAS-II Panic.
To determine whether baseline HKVs could predict response to treatment, we determined correlations between heat transfer and percent improvement in IDAS-II subscales in the subset of PT (n = 50) that completes post-treatment IDAS-II self-reports (calculated as , where pre and post are the IDAS-II subscales before and after 12 weeks of treatment, respectively). We first grouped both SSRI and CBT therapy cohorts of PT together to study common correlates of treatment response. Significant negative correlations with IDAS-II Depression and heat transfer between the left ACC/mPFC and IFC (ρ = −0.45, p = 0.001, q = 0.04), left auditory association cortex (ρ = −0.42, p = 0.002, q = 0.04) and INS/FOP (ρ = −0.41, p = 0.003, q = 0.04) were found (Figure 3B). Similar to baseline symptom correlations, no significant associations were found using IDAS-II Panic.
Next, we segregated PT by cohort to investigate treatment specific HKV predictors of therapeutic response. Interestingly, in the SSRI cohort (n = 22), significant negative correlations were found only with IDAS-II Panic improvement and heat transfer between the left ACC/mPFC and IFC (ρ = −0.60, p = 0.003, q = 0.15), INS/FOP (ρ = −0.54, p = 0.009, q = 0.17) and right NS/FOP (ρ = −0.53, p = 0.012, q = 0.17) (Figure 3C). In the CBT cohort (n = 28), we observed significant negative associations with IDAS-II Depression percent improvement and heat transfer between the right ACC/mPFC and superior parietal (SPC) (ρ = −0.51, p = 0.006, q = 0.09), inferior parietal (IPC) (ρ = −0.51, p = 0.006, q = 0.09) and left auditory association (ρ = −0.50, p = 0.008, q = 0.09) cortices (Figure 3D). We also found significant negative correlations with IDAS-II Panic percent improvement and heat transfer between the right ACC/mPFC and SPC (ρ = −0.55, p = 0.003, q = 0.12) and the left ACC/mPFC and right lateral temporal (LTC) cortex (ρ = −0.50, p = 0.007, q = 0.14) (Figure 3E). Interestingly, the brain regions involved significant correlations of treatment response in all PT and the SSRI cohort are all Frontal cortical regions that have been tied to emotion regulation in MDD [39, 66, 67], while the regions found with the CBT cohort are members of canonical RSNs. For example, the IPC and LTC are members of the default mode network (DMN), while the SPC is part of the dorsal affective network (DAN) [12, 27]. Please note that we increased the threshold for significance after FDR correction from q < 0.05 as for baseline and percentage improvement of IDAS-II subscales in all PT to q < 0.2 given the decreased sample size and decreased power of these correlations. As such, these results should be understood as exploratory and preliminary.
Identification of Subnetwork Targets for Modulation
Once we had characterized the disrupted information diffusion of SN1 in IP PTs, we next sought to determine brain subnetwork regions that may serve as neuromodulatory targets to normalize aberrant subnetwork information diffusion properties towards those found in HC. To this end, we start by computing the difference between the mean HC heat kernel matrix and each kth PT’s heat kernel matrix to obtain a residual heat kernel matrix. A favorable modulation strategy would then add supplemental diffusion activity to a PT’s heat kernel with a pattern that closely resembles the deficits encoded in a PT’s residual matrix. To identify such a strategy, the error of the difference between this residual matrix and the product of a heat value, , and a heat kernel modifier matrix, which encodes the effect of supplemental heat at a the ithnode on the network, is minimized over and repeated for all possible target nodes in order to identify an optimal brain region and heat value for stimulation. Once a modulation strategy is determined, the product of the modifier matrix and the heat value are then added to a PT’s heat kernel to model the effects of the intervention (see Methods section for details). For each potential target brain region, we then obtain a heat and error value corresponding to the optimal effectiveness of a brain region as a modulation target. As the error values for all target node tested were empirically observed to be uniform in distribution with small regional variance (Figure 4C), we defined optimality as the brain region requiring the smallest amount of supplemental heat. Furthermore, this allows us to identify brain region neuromodulatory targets that potentially more efficiently disseminate stimulation with fewer off-target effects from a greater modulatory energy.
We first conducted this analysis on the mean cortical regions of SN1 as target nodes and found the right ACC/mPFC, which was identified as a hub of SN1 by mean strSDD, to be the optimal region to target for heat supplementation, with the lowest mean heat across all PT (Figure 4A and C). In addition, we determined subject-specific optimal brain regions and the right ACC/mPFC was also found to be the modal optimal target (34 out of 65 PT) (Figure 4E). Next, we applied the same target node identification procedure to the individual brain regions of the right ACC/mPFC in SN1. We found that the right 8BM was the optimal nodal modulatory target (Figure 4B and D). Note that the right 8BM was an SN1 hub with the greatest strSDD of all individual subnetwork nodes.
Correction of Subnetwork Diffusion Impairment
As a final step in our analysis, we sought to confirm whether the heat modulation approach identified as above successfully normalizes subnetwork information diffusion in PT to that of HC statistically. To carry this out, we performed left sided t-tests (contrast: HC > PT) at each heat kernel value (HKV) between the brain regions of SN1 at baseline in HC versus PT. Only this direction of t-test was assessed at baseline as this modulation model can only add heat, and is not presently able to correct hyperconnected HKVs in PT. We then modified PT heat kernels at the optimal target nodes and heat values as described above and repeated the edgewise t-test analysis with left and right sided t-tests in order identify significantly increased HKVs in PT relative to HC as a result of heat supplementation. An HKV was defined as significant if it survived FDR correction (q < 0.05) with comparisons, where N is the number of brain regions in the subnetwork.
We conducted this analysis first on mean cortical regions using the right ACC/mPFC (N = 23) as a target for modulation and identified 24 significant baseline HKV differences (Table S2). Following heat supplementation, 19 of the original 24 HKV differences were no longer significant (Figure 4A, Table S2). Three HKVs remained significant in the same direction as baseline: connections between the right dorsal visual stream and the right auditory association area, INS/FOP and LTC. On the other hand, 2 HKVs gained new significance (PT > HC) following modulation; connections between the right ACC/mPFC (target) and the left ventral visual stream and the right IPC.
We next conducted the same analysis on individual brain regions in SN1 (N = 48), using the right area 8BM as the modulatory target (Figure 4B, Table S3). Of the 31 HKVs differences at baseline, 26 were corrected following heat modulation. Five HKVs remained significant following modulation, which predominantly follow the regional pattern as described above. All 6 of the newly significant (PT > HC) HKVs involved connections with the right area 8BM (target node).
DISCUSSION
This research provides preliminary evidence of a transdiagnostic subnetwork deficit, that resembles the cingulo-opercular network, characterized by information diffusion impairment of the bilateral ACC and mPFC. Central to this impairment is more specifically the right area 8BM, a key brain region involved in organizing a broad spectrum of cognitive tasks, which may underlie previously reported dysfunction of multiple brain circuits in the IPs. In addition, this is also the first report of using network diffusion models to study psychiatric disorders. We also demonstrate that models of neuromodulation involving targeting this brain region normalize PT diffusion dynamics towards those of healthy controls. These analyses provide a framework for multimodal methods that identify diffusion disrupted subnetworks and potential targets for neuromodulatory intervention based on previously well-characterized methodology.
From our analyses we discovered a subnetwork (SN1) with increased diffusion embedding distance, i.e., decreased information diffusion capacity at the time scale required for optimally capturing functional connectivity patterns, in HC relative to IP PT at baseline. Hub regions (as defined by HC or PT mean strSDD, indicating regions with the greatest total diffusion impairment) of SN1 include the bilateral ACC/mPFC and INS/FOP, which bears resemblance to the cortical aspect of the cingulo-opercular network (CON) (Figure 2). The CON, which includes the dorsal ACC, anterior INS, PFC, hypothalamus, thalamus and amygdala has been shown to contribute to many brain functions, including the processing of pain and negative affect, as well as the maintaining general cognitive control during goal-oriented behaviors [24]. Dysfunction of the CON has also been demonstrated in IPs, including MDD [37] and anxiety disorders [75]. The dorsal ACC, a core hub of the CON, has been shown to be critical in both MDD [81] and anxiety disorders [75]. In particular, investigations of emotion regulation in anxiety disorders have demonstrated that increased ACC activity [16] and coupling to anterior INS activity [46, 49] is present in HC relative to PT, indicating a potential regulatory role for the ACC during cognitive control of emotion processing within the CON.
SN1 also includes brain regions that are part of other functional brain networks, such as the LTC and IPC (a hub of SN1), part of the DMN, and the SPC, part of the DAN, in addition to the IFC, which is important for the coordination of executive function and emotion processing [66, 69]. Of note, the IPC is also a member of the frontal parietal network (FPN), an executive function network that antagonistically deactivates the DMN [7]. Increased activity of the DMN, a network defined by functional connectivity at rest associated with both rumination and worry, relative to other RSNs has been implicated in MDD [27] and anxiety disorders [45]. Studies of MDD have also indicated the widespread disruption in the salience and central executive networks [35, 43]. Indeed, accumulating evidence suggests that dysfunctional coordination between multiple brain networks may better explain the pathophysiology of IPs.
Our findings from interrogating interregional information transfer within SN1 indicate that connections between the bilateral ACC/mPFC and brain regions involved with canonical brain networks are associated with baseline symptom severity. IDAS-II Depression scores correlated positively with heat transfer between the left ACC/mPFC and the PCC. Reduced connectivity between the PCC, a DMN hub, and DLPFC [51] and other frontal regions [84] has been previously reported in MDD. In the macaque, retrograde tracing studies have revealed afferent connections to the 8BM subregions of the ACC/mPFC from the PCC [28]. Taken together, our results may indicate that depressive symptoms in the IPs at least partly result from the hyperconnectivity of the DMN to the MDN and CON, which may skew cognitive resources away from executive control of emotion processing. On the other hand, negative correlations were found with IDAS-II Depression and heat transfer between ACC/mPFC and the IFC and INS/FOP. These results are similar to those from previous studies, where the ACC activity was found to negatively correlate with depression symptom scales [81] and with anxiety symptom scales in the context of its potentially antagonistic role with the INS in the CON [46, 49]. Overall, these baseline associations provide evidence for ACC/mPFC activity and connectivity to regions important for emotion regulation as critical to transdiagnostic depression symptoms in the IPs.
Our results also indicate that internodal information transfer of the bilateral ACC/mPFC in SN1 is associated with therapeutic response in both treatment modality agnostic and specific manners. Interestingly, predictors of general and SSRI specific treatment response involved similar brain regions to those associated with baseline symptom severity. In all PT, negative correlations were found with IDAS-II Depression symptom improvement and the heat transfer between the left ACC/mPFC and the IFC, INS/FOP and auditory association cortex. As with all PT at baseline, no associations were found with IDAS-II Panic. In the SSRI cohort, however, significant correlations were only found using IDAS-II Panic improvement, yielding negative associations with connections between similar regions as with all PT: heat transfer between the left ACC/mPFC and the IFC and bilateral INS/FOP. These findings indicate that lower baseline information transfer between the ACC/mPFC and CON or IFC are predictive of general and SSRI specific responsiveness to treatment, and these connections may be the substrate of treatment action. Interestingly, associations of baseline internodal information transfer and treatment response in the CBT cohort were found with unique brain regions compared to all PT and the SSRI cohort. Significant associations were also found using both IDAS-II Depression and Panic scales. Depression improvement was negatively associated with information transfer between the right ACC/mPFC and SPC and IPC, regions of the dorsal affective network and DMN. Using IDAS-II Panic improvement, we found significant negative correlations with information transfer between the right ACC/mPFC and SPC the left ACC/mPFC and right LTC, another DMN subnetwork region. These results indicate the presence of CBT specific neural substrates of treatment prediction, and again, that lower baseline information transfer between these regions is indicative of greater treatment efficacy. In addition, these findings are in line with a previous study of CBT for social anxiety disorder, where increased baseline activation of the ACC and LTC [47] was predictive of treatment response.
Central to the diffusion impairment of SN1 is a subregion of the right ACC/mPFC cortical region, area 8BM, the caudal aspect of the dorsomedial PFC (dmPFC), which borders the dorsal aspect of the ACC. Retrograde tracing studied in the macaque monkey have revealed both afferent and efferent neuronal connections between 8BM and the ACC, indicating likely functional synergy of these areas [28, 57]. Importantly, 8BM is brain region that has recently been found to be a key region of the multiple-demand subnetwork (MDN), a fronto-parietal system that is co-activated during a broad range of cognitively demanding tasks [4]. As such, the MDN is likely critical for the organization and recruitment of multiple task-specific subsystems in the brain towards current cognitive needs. Additionally, our network-diffusion model identified the most effective brain regions as putative targets for neuromodulatory stimulation on a subject-specific level. Examining the results of this investigation provide evidence that the right ACC/mPFC is an optimal target for neuromodulation to correct SN1 information diffusion dynamics of PT towards those observed in HC, and, further, that area 8BM is the most efficient subregion of the ACC/mPFC for such an intervention.
These preliminary findings have important clinical implications for the therapeutic application of neuromodulatory interventions for IPs, such as repetitive transcranial magnetic stimulation (rTMS), a focal, non-invasive brain stimulation method. rTMS, which is typically delivered to the dorsolateral PFC (dlPFC), is a first line treatment for SSRI-refractory MDD [33] but has also been shown to significantly reduce anxiety [18, 25] and nicotine dependence symptoms [3]. Although rTMS is generally effective, MDD remission following treatment is 30-40% [33], indicating a need for developing better therapeutic paradigms, such as personalized treatment protocols. In line with these goals, a study that used MRI guidance to enhance targeting of specific regions of the dlPFC, all rTMS-resistant MDD PTs responded to this new approach [58]. Another study used fMRI guidance to individually targeted cortical emotion regulation systems to improve rTMS efficacy [53]. Our results offer an additional strategy for identifying brain network targets of rTMS treatment for correcting subject-specific structural deficits that may underlie functional expression of IP symptoms.
Limitations and Future Directions
A significant question about the implications of our findings is whether the addition and propagation of heat on a subnetwork sufficiently models the effects of rTMS stimulation. It has been shown that the propagation of direct cortical electrode [74] and rTMS stimulation and subsequent functional activity are best predicted by structural connectome topology [11, 56], providing encouraging evidence that network diffusion-based analyses such as those proposed in this study are an appropriate model for the effects of rTMS neuromodulation. Nevertheless, further studies are necessary to verify the utility of the proposed model. In addition, while this is the first report of using a multimodal information diffusion model towards a transdiagnostic sample of IPs, the findings presented in this paper pertaining to the pattern of diffusion impairment of SN1 should be replicated in a larger sample of IP PT and HC to determine the generalizability of our results. This is especially true for the analyses obtained by segregating PT into SSRI and CBT cohorts, as the sample sizes of these groups are very small, and, accordingly, these findings must be interpreted as preliminary and exploratory. Further, availability of post-treatment brain imaging would be required to make reliable conclusions about network diffusion-based substrates of treatment response, and, as such, should be included in future studies.
Although our structural to functional mapping quality, as measured by Pearson correlation between the predicted and empirical functional connectomes, is in the range of values obtained in a replication study of related mappings, other network-diffusion related models have been developed that result in better predictive mappings [23]. Future studies of these proposed methods should investigate how the mapping parameters from these related methods can be incorporated into the currently proposed model for better structural-functional fusion. Another area for further development in our proposed network diffusion methods is that neuromodulation target identification and simulation can only model the effects of the addition of heat to subnetwork’s diffusion dynamics. Therefore, if a brain circuit is found to be pathologically hyperconnected, other strategies will have to be developed in order model either direct or indirect inhibition of brain activity.
Future implications of this research are broad. We present here an analytical framework assembled from well-characterized neuroimaging and graph theoretical methods that can be used as is to study multimodal brain networks for other brain disorders. Aside from the application to other datasets and to larger similar datasets for replication, this model offers the ability to determine the diffusion embedding from functional connectivity other than that observed during resting-state. For example, functional connectomes can be generated from tasks that are known to be implicated in a certain disorder, and as a result the diffusion embedding will be formed from the structural basis that best maps to the functional connectivity observed during this task. This allows for the researcher to investigate the structural connectivity potentially most pertinent to forming the functional brain activity observed during specific tasks, and therefore permits a more granular interrogation of complex features of brain disorders and states. Perhaps the most promising and immediate application of our proposed methodology is that of identifying brain regions best suited for therapeutic for neuromodulatory intervention, such as with rTMS. As the proposed methods provide both subject-specific regional targets and magnitudes of stimulation that best modify subnetwork dynamics towards those of a desired (HC) network, future work is in line with improving the outcomes of rTMS intervention by personalizing treatment features.
Conclusion
There are two major outcomes of this study. First, we report impaired information diffusion between area 8BM and other SN1 regions, many of which have been previously implicated in both the pathology of multiple IPs and the function of the CON, DMN and DAN. Second, we found that hubs of SN1, found to be critical for organizing brain function for a variety of cognitive and emotional processes, are optimal targets for modelled neuromodulatory intervention. Taken together, our results may indicate the presence of a concerted disruption of multiple brain networks pertinent to cognitive control of emotion regulation in IPs. This dysregulation of connectivity could result in a loss of “top-down” executive control of emotion processing via connections between the MDN and other task-specific (and task-negative) brain networks. Such a perturbation in brain network dynamics involving the integration of multiple complex subsystems via a hub of the MDN may represent the underlying pathophysiological brain network features that are common to all IPs and give rise to the heterogeneous expression of transdiagnostic symptoms.
Data Availability
Data and software specific to this manuscript are available at the lead author's github page. All other resources are detailed in the methods section of the manuscript.
METHODS
Resource table
Clinical trial and research participants
Subjects were recruited from the greater Chicago area through advertisements and through University of Illinois at Chicago (UIC) outpatient clinics and counselling centers as part of a larger Research Domain Criteria (RDoC) [22] investigation on predictors of IP treatment outcomes (ClinicalTrials.gov identifier: NCT01903447). A heterogeneous study population was recruited in order to obtain a sample with a broad range of symptom severity and functioning. Details regarding inclusion/exclusion criteria, participant recruitment, clinical characteristics and treatment have been previously described [36]. In brief, this study was approved by the UIC Institutional Review Board, and written informed consent was obtained for each participant. The inclusion criteria for subjects were age between 18 and 65 years, and the need for randomization to 12 weeks of treatment with SSRI or CBT, as determined by a consensus panel consisting of at least three trained clinicians or study staff. Subjects were excluded from the study if they have a history of current or past manic/hypomanic episodes or psychotic symptoms, active suicidal ideation, presence of contraindications or history of SSRI resistance (no response to >2 SSRIs despite adequate duration and dose), psychopathology not appropriate for the treatment algorithm, or current cognitive dysfunction or impairment. The SCID-5 [30] was used to determine current and lifetime Axis I diagnoses. The study was a parallel group randomized control trial with 1:1 allocation ratio to either 12 weeks of CBT or SSRI. For the SSRI cohort, PTs were administered one of 5 drugs (sertraline, fluoxetine, paroxetine, escitalopram or citalopram) with a flexible dosing schedule with a goal of obtaining target dose by 8 weeks. SSRI PTs met at 0, 2, 4, 8 and 12 weeks with their study psychiatrist for medication management. For the CBT cohort, PTs received 12 once-weekly 60 min sessions led by a PhD-level clinical psychologist. CBT procedures were based on the PT’s principal diagnosis and predominant symptoms [15]. Each participant was scanned at enrollment and IP subject scans were acquired before treatment was administered.
At the time of enrollment (Pre) and after 12 weeks of treatment (Post), severity of IP symptoms was assessed in all subjects using the Inventory of Depression and Anxiety Symptoms (IDAS-II) [79]. Subjects responded to each of the 99 items in this inventory using a 5-point Likert-type scale ranging from 1 (not at all) to 5 (extremely), yielding 17 empirically derived and symptom-specific scales (Suicidality, Lassitude, Insomnia, Appetite Loss, Appetite Gain, Ill-Temper, Well-Being, Panic, Social Anxiety, Traumatic Intrusions, Traumatic Avoidance, Mania, Euphoria, Claustrophobia, Checking, Ordering, Cleaning, General Depression and Dysphoria). As we are interested in transdiagnostic NVS construct disruptions in IPs, we use the IDAS-II Panic and Depression subscales, as these scales have been shown to map well to ‘fear’ and ‘distress’ dimensions, respectively, which is a previously used approach for broadly dividing and assessing these symptom domains [63, 64, 79].
Image acquisition and processing
Structural MRI
A diffusion tensor imaging (DTI) scheme was used, and a total of 32 diffusion sampling directions with 4 b0 images were acquired on a 3 Tesla GE Discovery MR750 System (Milwaukee, WI) at the UIC Center for Magnetic Resonance Research. The b-value was 1000 s/mm2. The in-plane resolution was 0.9375 mm. The slice thickness was 2.5 mm. Two sets of DTI images were acquired, with opposite phase encoding directions. For preprocessing, the two sets of DTI images were realigned and corrected for motion artifacts and eddy current induced distortions using the FMRIB Software Library (FSL) [42]. White matter tract reconstruction was done using DSI Studio software [85]. A deterministic fiber tracking algorithm was used with whole brain seeding with a total of 10000000 seeds, an angular of 70 degrees, step size of 1 mm and anisotropy threshold of 0.1. The fiber trajectories were smoothed by averaging the propagation direction with 10% of the previous direction. Tracks with length shorter than 10 or longer than 300 mm were discarded. Cortical areas were defined by the MMP1.0 parcellation comprising of 360 regions of interest (ROIs), 180 per hemisphere, from Glasser et al. [34]. A 360 x 360 structural connectivity matrix or connectome was created for each subject where each edge is given by the fiber count of reconstructed white matter tracts, normalized by median tract length, between ROIs.
Functional MRI
Whole-brain blood-oxygen-level dependent (BOLD) functional images were acquired using a T2* weighted gradient-echo echo-planar imaging (EPI; 2s TR, 25ms TE, 82° flip, 64 × 64 matrix, 200 mm FOV, 3 mm slice thickness, 0 mm gap, with 44 axial slices) sequence optimized to reduced susceptibility artifact, along with a high-resolution T1 structural scan. During this scan, subjects were asked to view a fixation cross on a blank background for 5 minutes. Subjects were instructed to keep their eyes open and focused on the cross, and to try not to think of anything in particular for the duration of the scan. Functional connectivity was measured using the resting-state fMRI toolbox, CONN [80], using standard preprocessing and denoising pipelines as described in [62]. The principal components of the white matter, CSF signal, realignment parameters, and scrubbing were regressed out of the signal using the CompCor method in CONN [8]. BOLD signal data was passed through a band-pass filter of.008 to.09 Hz. Using 360 regions of interest (ROIs) defined by the MMP1.0 multimodal parcellation [34], functional connectivity measures were derived using pairwise Pearson correlations on the BOLD time series data between each ROI. A 360 x 360 functional connectivity matrix or connectome was created using these Pearson r statistics for each subject and used for the structure-to-function mapping.
Network notation
A brain network may be represented as a graph with nodes being grey matter regions and edges being the connection between these regions. For structural connectomes, edge weights are assigned based on the number of fiber counts between nodes, normalized by the median fiber length. For functional connectomes, edge weights are the Pearson correlation r of the time series of BOLD signals between nodes. Formally, a graph is defined as G = (V, E) where V is the set of nodes of size N and E is the set of edges linking nodes in V. In addition, w : E → ℝ, is a weight function that assigns weights to the E according to the modality of imaging from which the brain graph is constructed. G can then be encoded as an adjacency matrix, A ∈ ℝNxN, where each entry Ai,j corresponds to the connection weight between nodes vi and vj. The diagonal strength matrix is then defined as , which can be used to define the graph Laplacian matrix, L = D − A, and, the normalized graph Laplacian is defined as . The eigendecomposition of the graph Laplacian is given by ℒ = U ΛU T, where U is a matrix with columns as eigenvectors and Λ is the diagonal matrix of eigenvalues. The spectral properties of the normalized graph Laplacian have been extensively studied [21].
Diffusion in networks
Diffusion in networks over time can be determined analytically using the heat equation given by where ℋ(t) is the heat kernel and fundamental solution to the heat equation, and t is time. The heat kernal can be understood as describing the flow of information between all nodes of a network over time, taking into account all possible pathways of information flow into and out of all nodes in the network. Throughout this manuscript, heat or information with be used interchangeably in this context. The heat kernel is then given by which can then be used to solve for the heat distribution on network nodes, h(t) after time = t, given an initial condition h(0):
The product of an element of the heat kernel and initial condition vector, ℋ(t)i,j * h(0)i, then represents the amount of heat at the jth node at time = t that has diffused from the ith node. The heat kernel can be computed as a sum of the outer product of eigenvectors of the graph Laplacian, scaled by exponentiating the corresponding eigenvalues with time:
Structure to function mapping
In this study, we use the methods proposed by Abdelnour and colleagues [2] [1] to identify the graph diffusion-based mapping of functional connectivity from the structural basis. Briefly, the observed functional activity in rs-fMRI, i.e., the transfer of information from brain region (node) i to region j as measured by pairwise correlation of temporal BOLD signals, can modeled by first order diffusion-like dynamics given by: where β is a diffusion constant. This extends to the entire network then as the heat equation, . The mapping from the structural diffusion dynamics, yielding the estimated functional connectome, 𝒞est,is then given by [2] and has been updated as [1] where a and b are additional model parameters and I is the identity matrix. Note that the summation allows for the exclusion of eigenvector-eigenvalue pairs. As in [1], we leave out the first eigenvector with corresponding zero valued eigenvalue, as it largely represents uniform background connectivity (captured by b parameter above) and is typically regressed out of rs-fMRI data. The model is then fit by minimizing the normalized predictive error with respect to the model parameters given by: where ℱ is the empirical rs-fMRI connectivity matrix. This model is fit for each individual to obtain subject specific βt parameters, used in the diffusion-based embedding discussed below.
Diffusion-based network embedding
The embedding of nodes based on the diffusion properties of a graph has been studied extensively for dimensionality reduction and clustering of multi-dimensional data [9, 54, 61]. These approaches center around the embedding of network vertices via the eigenvectors of the graph Laplacian. Each element of an eigenvector corresponds to the coordinate of the corresponding node such that nodes that are closer together by geodesic distance on the underlying graph or manifold have more similar coordinate values, i.e., are closer together in Euclidean distance, in the embedding. A subset of these eigenvectors of size k = 1, 2, …, N can then be used to embed nodes in ℝk. In the context of this manuscript, then, brain nodes that are able to more efficiently pass information between one another are embedded closer together via the Laplacian eigenmodes of a brain network [78].
To compute the temporally dependent diffusion-based embedding for a network, we follow the methods discussed by Xiao and colleagues, where node coordinates are obtained from the Young-Householder decomposition of the heat kernel [82, 83]. The embedding matrix, Y = (y1| y2|…| YN), with columns as embedding coordinate vectors for network nodes can determined with the heat kernel by ℋ(t) = YT Y. The Euclidean distance between nodes i and j is then where the distance in each dimension is scaled by exponentiating the product of the corresponding eigenvalue and diffusion time, t. A pairwise Euclidean distance matrix, 𝒟, of embedded nodes is computed for each subject with the time parameter βt obtained via the structure to function mapping as described above. This provides a newly defined structural connectome where edges are distances in the embedding space spanned by the eigenvectors of the graph Laplacian, and the scale of distance in each embedding dimension is dependent on the time parameter, βt, for each subject:
Subnetwork identification
The structural diffusion distance (SDD) connectome edges then represent the distance between brain regions in diffusion embedding at the diffusion depth that best represents empirical resting-state functional connectivity. Many previous findings from studies of IPs implicate perturbations in resting state functional connectivity, especially within the DMN, which is defined by activity at rest. As such, SDD connectomes then used for further analysis to identify structural subnetworks with aberrant diffusion characteristics pertinent to rs-fMRI dynamics using the Network-based Statistic (NBS) algorithm [86]. Briefly, NBS is carried out by computing a t-statistic at each network edge, applying a pre-determined threshold to the resulting t-statistics and determining the connected components formed by supra-threshold edges. To determine the significance of each identified subnetwork, the size of each component is compared to a null distribution of maximum sized connected components obtained by shuffling group labels for t-statistic calculation. NBS was carried out using SDD connectomes in HC versus PT at baseline, using a left sided t-test with t-statistic thresholds of t = {3.0, 3.5, 4.0, 4.5}. In doing this, a subnetwork, S ⊂ G, where G is the full brain network graph, with significantly greater embedding distances, or significantly decreased diffusion capacity, in PT relative to HC was.
NBS subnetwork hubs
Hub node identification for brain regions within the significant NBS subnetwork was performed using edges from both the SDD and standard fiber count structural connectomes, and hub or centrality metric values were then averaged by group. For SDD connectomes, strength, or weighted degree, strSDD(n), for the ith node, ni, in S is defined as the sum of all edge weights (distances) that are within the identified subnetwork, S, i.e.,
A high strength value would then indicate that a node has relatively less diffusion between other nodes in the subnetwork, thus identifying nodes that may have the greatest diffusion impairment. For structural connectomes, standard measures of nodal strength, betweenness centrality, local efficiency and clustering coefficient [13, 72] were computed on the whole brain network graph and subgraph, S ⊂ G, for each node in S. For each of the above scenarios, metrics were also determined by cortical region by averaging region of interest (ROI) level metrics by their cortex assignment as defined by the MMP1.0 360-node parcellation [34].
Heat kernel edge-based analyses
In order to study more granular characteristics of subnetwork diffusion, heat kernel values between nodes within the significant NBS subnetwork were investigated using t-tests in HC vs PT at baseline. Correlates of baseline symptom severity (using PT only) of and treatment response (defined as were found by computing non-parametric Spearmen rho statistics between symptom scales and heat kernel values. As above, where cortex-averaged metrics were computed, heat kernel edge-based analyses were conducted cortex-averaged heat kernels in order to simplify interpretation and reduce the number of mass univariate tests performed. To further focus this analysis, statistics other than baseline t-statistics were only computed on heat kernel values corresponding to pairwise links between SSD hubs (as determined above) and all other subnetwork nodes. P-values for all statistics were corrected for multiple comparisons using False discovery rate correction (FDR) [10].
Subnetwork targets for supplemental heat
To identify potential brain regions for neuromodulatory treatment, we model a brain subnet-work receiving supplemental heat by simply adding a heat kernel modifying matrix, whose rows are made by repeating rows of the original heat kernel, to the original heat kernel. As discussed previously, the product of an entry of the heat kernel matrix, ℋ(t)[i, j], and the heat at node i at time t = 0, h(0)i, indicates the amount of heat transferred from node i to node j at time t. The product of the ith row of the heat kernel matrix, ℋ(t)[i, :], and h(0)i then yields the distribution of heat values for each jth node in the network. If the ith node of a network is supplied with supplemental heat independent of the initial distribution of heat on the network, h(0), the network heat distribution at time t is given by where is a vector of all zeros except for the ith whose value is the amount of supplemental heat at node i. Using this model, we can then identify both a node and supplemental heat value that may best ‘normalize’ a heat kernel representing impaired diffusion processes to a reference heat kernel. In this study, we use the mean HC heat kernel matrix as the reference and identify patient-specific nodes and supplemental heat values for optimal correction of the diffusion dynamics encoded in the heat kernel, described as follows.
Given the mean HC heat kernel, ℋHC, and the kth patient’s heat kernel, , we compute the residual heat kernel matrix as . For the ith node in the subnetwork, we construct the heat kernel modifier matrix as as the outer product of the all ones vector and the ith row of the kth patient’s heat kernel matrix. In order to find the optimal supplementary heat value, for a given node and subject, we minimize the following objective function.
Which can be easily solved analytically by where tr(·) denotes the matrix trace.
It is important to note, that the above model is also constrained such that only positive values in the residual matrix are used for optimization, as this simple network diffusion model only allows for the addition of heat.
This process is repeated for each unique node-patient pair, and the mean or patient-specific optimal node and corresponding supplemental heat value can be determined. We have experimentally observed that the minimum error between residual and heat kernel modifier matrices has very little variation between brain regions, and we define optimal region as the node with the lowest value supplemental heat value, such that neuromodulatory efficacy could be achieved with lower amounts of stimulation and thus potentially yield fewer undesired effects.
Visualizing effects of heat supplementation
To assess whether the significant group differences at baseline between HC and PT heat kernel values are corrected with heat supplementation, we add the product of the optimal supplemental heat value, , at the ith node, and the heat kernel modifier matrix, to each kth patient’s heat kernel. We then compute a t-statistic for each heat kernel entry, as discussed above.
ACKNOWLEDGEMENTS
This study was supported by funding from the National Institute of Mental Health of the National Institutes of Health (NIMH-NIH) grants R01MH101497 (to KLP) and the NIMH grant 5T32MH067631-14 (support for PJT).
Footnotes
This version of the manuscript contains revisions for clarity, grammatical errors and improved figures to better represent results.
References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].↵
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].↵
- [63].↵
- [64].↵
- [65].↵
- [66].↵
- [67].↵
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].↵
- [82].↵
- [83].↵
- [84].↵
- [85].↵
- [86].↵