Predicting Dimensional Antidepressant Response to Repetitive Transcranial Magnetic Stimulation using Pretreatment Resting-state Functional Connectivity

Repetitive transcranial magnetic stimulation (rTMS) is an effective treatment for depression and has been shown to modulate resting-state functional connectivity (RSFC) of depression-relevant neural circuits. To date, however, few studies have investigated whether individual treatment-related symptom changes are predictable from pretreatment RSFC. We use machine learning to predict dimensional changes in depressive symptoms using pretreatment patterns of RSFC. We hypothesized that changes in dimensional depressive symptoms would be predicted more accurately than scale total scores. Patients with depression (n=26) underwent pretreatment RSFC MRI. Depressive symptoms were assessed with the 17-item Hamilton Depression Rating Scale (HDRS-17). Random forest regression (RFR) models were trained and tested to predict treatment-related symptom changes captured by the HDRS-17, HDRS-6 and three previously identified HDRS subscales: core mood/anhedonia (CMA), somatic disturbances, and insomnia. Changes along the CMA, HDRS-17, and HDRS-6 were predicted significantly above chance, with 9%, 2%, and 2% of out-of-sample outcome variance explained, respectively (all p<0.01). CMA changes were predicted more accurately than the HDRS-17 (p<0.05). Higher baseline global connectivity (GC) of default mode network (DMN) subregions and the somatomotor network (SMN) predicted poorer symptom reduction, while higher GC of the right dorsal attention (DAN) frontoparietal control (FPCN), and visual networks (VN) predicted reduced CMA symptoms. HDRS-17 and HDRS-6 changes were predicted with similar GC patterns. These results suggest that RSFC spanning the DMN, SMN, DAN, FPCN, and VN subregions predict dimensional changes with greater accuracy than syndromal changes following rTMS. These findings highlight the need to assess more granular clinical dimensions in therapeutic studies, particularly device neuromodulation studies, and echo earlier studies supporting that dimensional outcomes improve model accuracy.

Depression has long been characterized as a network-based disorder involving disrupted connectivity and metabolic activity of the subgenual anterior cingulate (sgACC), hippocampus and DLPFC [8][9][10]. While modulation of sgACC activity is widely associated with successful antidepressant response [11], rTMS is only capable of directly stimulating super cial cortical regions, though these act as windows that allow transsynaptic distal and indirect modulation of connected nodes. Targeting of the left DLPFC is therefore motivated by its relative accessibility as a more super cial target within this putative depression network, and its downstream connections with nodes relevant to mood disorders: successful rTMS response is thought to be mediated, in part, by increased inhibition of the sgACC through the modulation of the DLPFC [2,8,12]. Despite the hypothesized network-level mechanisms of antidepressant response to rTMS, few studies have investigated whether rTMS treatment outcomes are predictable using pretreatment MRI-based resting-state functional connectivity (RSFC) measures. Among these, a recent study investigating RSFC predictors of rTMS in 47 patients reported that regional blood oxygen level dependent (BOLD) signal power and connectivity of the default mode network (DMN) could predict treatment response with 85-95% accuracy [13]. Additional correlative ≥ studies have reported that the pretreatment volume of the left superior frontal and caudal middle frontal cortices [4] and thickness of the left rostral ACC [14] are predictive of rTMS response. Another study using accelerated intermittent theta burst stimulation (aiTBS) in 50 patients to treat depression identi ed a correlation between pretreatment caudal cortical thickness measures ACC and subsequent treatment response [15]. These studies suggest that structural or functional properties of networks and regions commonly affected by depression might be predictive of rTMS response.
The symptomatic heterogeneity of depression is a well-known obstacle to the identi cation of treatmentpredictive biomarkers. Concretely, a DSM-based diagnosis of major depressive disorder (MDD) requires ve out of nine symptoms allowing for 227 potential symptom constellations to arrive at the same diagnosis (i.e., equi nality), and even the possibility that 2 patients diagnosed with MDD do not have overlapping symptoms. This varied range of potential symptom constellations may be underpinned by dysfunctions in unique and/or overlapping neural circuits. For example, a review of 59 fMRI studies on patients with current or remitted depression reported that two separable aspects of anhedonia, reward learning and reward liking, are characterized by reduced frontostriatal sensitivity to positive feedback and nucleus accumbens/ventral striatal hypoactivation [16]. This issue is compounded in biomarker studies that use syndromal severity scale total scores as primary outcomes, as total scores rarely re ect nuanced patterns of individual symptomatology. Additionally, there is no speci c brain region or network functionally relevant for all 9 symptom dimensions captured in the DSM nosology, and hence correlations of biological and clinical measures will necessarily be plagued by low signal-to-noise scenarios. Strategies to contend with this problem have generally used datadriven approaches to either cluster patients into more homogenous biotypes on the basis of neuroimaging or genetic features [17,18] or to identify homogenous symptom dimensions to model as primary outcomes rather than total scores. Using this latter strategy, we have previously shown that neuroimaging biomarkers of antidepressant outcomes in electroconvulsive therapy (ECT) and serial ketamine infusion (SKI) can be improved by modeling dimensional components of depressive symptoms derived by applying exploratory factor analysis (EFA) to the 17-item Hamilton Depression Rating Scale (HDRS). Speci cally, the EFA model identi ed dimensions of core mood and anhedonia (CMA), somatic disturbances (SoD), and insomnia [19,20]. Here, we used machine learning to predict changes in depressive symptoms captured by the HDRS-17 and HDRS-6, as well as the CMA, SoD, insomnia subscales of the HDRS-17, and, as a secondary outcome, the Snaith-Hamilton Pleasure Scale (SHAPS) following rTMS using pretreatment RSFC measures as predictive features. We anticipated that the change in dimensional symptoms of the HDRS would be predicted more accurately than changes in the HDRS-17 total score, as we observed in ECT and SKI studies. We expected that pretreatment connectivity of the DMN and frontoparietal control network (FPCN) would be particularly predictive of treatment outcomes across scales and subscales and that outcomes would be predicted by overlapping but unique patterns of RSFC.

Materials and Methods
Participants Twenty-six patients experiencing a depressive episode (mean age = 41.1 14.1 years, 50% male) diagnosed by the Structured Clinical Interview for DSM-IV were assessed between February, 2015 and November, 2020. Patient ± clinical and demographic characteristics are outlined in Table 1. All patients provided written informed consent and the study procedure was approved by the local ethics committee of the Massachusetts General Hospital.

Image Acquisition and Processing
Images were collected on a 3-Tesla Siemens Skyra MRI scanner (Siemens Healthcare, Malvern, PA, USA). Structural data was acquired with an anatomical T1-weighted multi-echo magnetization prepared rapid gradientecho sequence with parameters: Repetition time (TR) = 2530ms, echo time (TE) = 1.69ms, inversion time (TI) = 1100ms, ip angle = 7.0°, number of excitations = 1, slice thickness = 1mm, eld of view (FoV) = 256mm, in-plane resolution = 1.0 x 1.0mm, and a matrix of 256x256. Resting-state BOLD data was collected with a whole-brain echoplanar imaging sequence with the following parameters: TR = 3000ms, TE = 30ms, ip angle = 85°, slice thickness = 3.0mm, in-plane resolution = 3.0x3.0mm, FoV = 216mm. Subjects were instructed to keep eyes open and try to stay focused passively on a xation cross displayed behind them during the resting-state scan (white cross on black background).
FMRI data were preprocessed using fMRIprep software package v. 20.2.1 [21]. Anatomical images were corrected for intensity bias with Advanced Normalization Tools N4BiasFieldCorrection (ANTs v.

2.3.3) (22).
Images were segmented across tissue types in both ANTs and Freesurfer, and the previously estimated brain mask was re ned by reconciling these two tissue segmentations using Mindboggle [23]. The T1 images were warped to MNI space using nonlinear registration with antsRegistration (ANTs 2.3.3).
Functional images were skull-stripped and aligned to the anatomical reference image using boundary-based registration with six degrees of freedom in Freesurfer. Functional images were then slice-time corrected (AFNI v. 201660207) [24] and motion-corrected (FSL v. 5.0.9) [25]. These images were then normalized to MNI space, using the warp previously computed for transforming the T1 image to MNI space. Fmriprep's implementation of anatomical CompCor[26] was used to estimate principle components for signal originating in white matter and CSF, respectively. The data were then blurred using a 5mm smoothing kernel and submitted to confound regression using nilearn.glm. rst_level. This regression included regressors for six motion parameters and their derivatives, and enough aCompCor principle components to capture 50% of the variance of signal from white matter and CSF, respectively (a maximum of ve components). Regressors were also added to scrub the rst 3 TRs and high motion TRs. High motion TRs were de ned as those having a framewise displacement over .3 mm.
Runs with greater than 50% of TRs meeting this threshold were excluded from analysis. Sessions with fewer the 5 minutes of data after censoring were excluded.

Predictive Features
Pretreatment rs-fMRI scans were parcellated using multiple atlases. We included 200 cortical regions using the Schaefer atlas [27], 24 regional parcellations from the Kelly insula atlas [28], 21 subcortical parcellations from the Harvard-Oxford Atlas [29], 16 nuclei from the Pauli subcortical atlas [30], 10 parcellations from the Pauli amygdala nuclei atlas [31], the bilateral bed nucleus of the stria terminalis, and the grey matter of the TMS target for 274 regional measures. Supplementary Table 2 lists each predictor. Global connectivity (GC) values for each ROI were computed subjectwise using the graph theoretic node degree as previously reported in [20]. NiLearn [32] scripts were used to compute correlations of time series data between each ROI. Subjectwise correlation matrices were thresholded at Z |0.4| to create a binary adjacency matrix. The graph theoretic node degree for each parcellation was computed as the number of other regions with which the given region was correlated above the threshold of Z |0.4|, serving as a proxy of GC for each ROI.

Outcome Measures
Primary outcomes included changes along the depression severity HDRS-17 total score [33] and the HDRS-6 subscale [34]. We also included a previously identi ed three-factor solution [19,35]

Predictive Modeling
Random forest regression (RFR) models were trained to predict change in each outcome. Models were trained and tested using 10-repeated 10-fold cross validation with a nested grid search for parameter optimization. The parameter grid tuned the number of features selected, k, using mutual information regression with k={10, 50, 100, p}, where p is the number of all predictors; number of regression trees n_tree={50, 100, 500, 1000}. Model performance was evaluated using the sums of squares formulation of ; i.e., 1 -, where is the predicted outcome of the i-th subject and is the average outcome across all subjects in the testing fold [37].
The signi cance of model performance was tested using permutation tests with B = 100 permutations. Model signi cance was adjusted for multiple comparisons across primary outcomes using a Bonferroni correction,

Post-hoc Analysis
Components of the DMN and dorsal attention network (DAN) were predictive of CMA, HDRS-17, and HDRS-6 outcomes and the activity of both networks have been associated with ruminative symptoms [39][40][41]. Therefore, as post-hoc analyses, we conducted Pearson correlation tests to determine whether the pretreatment Rumination Response Scale (RRS) [42] total score and brooding and re ection subscales were associated with (1) changes in the CMA, HDRS-17, and HDRS-6 scales, and (2) with the top six most predictive GC measures across those outcomes.

Cohort Characteristics
With treatment response de ned as 50% reduction in the HDRS-17 total score, 42% of the cohort were treatment responders. Female patients experienced greater symptom reduction than males (t = 2.10, p < 0.05). Patient age was not signi cantly associated with the degree of HDRS-17 symptom reduction. The degree of symptom reduction across all outcomes is reported in Table 1.

Model Performance
Change along the CMA dimension was predicted signi cantly above chance ( =0.09, p < 0.01). The HDRS-6 subscale ( =0.02, p < 0.01) and HDRS-17 total scores =0.02, p < 0.01) were also predicted signi cantly. The mean for the CMA outcome was signi cantly higher than the HDRS-17 (p < 0.05). The remaining primary and secondary outcomes (SoD, insomnia, and SHAPS) were not predicted signi cantly (all p > 0.05). Model performance is outlined in Table 2 and distributions are illustrated in Fig. 1.    Figure 2(a) illustrates regional parcellations predictive of change for the CMA outcome. Table 3 outlines characteristics of regional predictors for each signi cantly predicted outcome.

Post-hoc Evaluation of Ruminative Symptoms
Pearson correlation tests identi ed no associations between the pretreatment RRS total score, Brooding, or Re ection subscales and changes in the CMA, HDRS-17, or HDRS-6 scales or GC measures predictive of changes in those outcomes (all p > 0.05).

Discussion
We investigated whether dimensional changes in depressive symptoms following rTMS were predictable using multivariate patterns of pretreatment RSFC. Changes along the CMA dimension and HDRS-17 and HDRS-6 scores were predicted signi cantly and our ndings are convergent with earlier studies in ECT and SKI which suggested that change in the CMA dimension is predicted more robustly than other symptom clusters. and HDRS-6 total scores were predicted by two GC patterns that overlapped with the CMA dimension: The right primary motor cortex (spanning BA4) and the right fusiform gyrus (spanning BA37). The GC of the right supramarginal gyrus (BA40) within the FPCN was a common predictor for both changes in the HDRS-17 and HDRS-6 outcomes. Beyond these, important predictors were unique to each outcome; however, unique parcellations within overlapping large-scale resting-state networks exhibited similar relationships with symptom changes across primary outcomes. For example, increased GC of the DMN and SMN predicted poorer symptom changes across all primary outcomes while increased GC of components within the DAN and VN generally predicted more symptom reduction.
Dimensional clinical measures are predicted with greater accuracy than syndromal severity outcomes A central aim in precision psychiatry is the identi cation of biomarkers predictive of individual outcomes following antidepressant treatment. Identi cation of such biomarkers is obscured by the symptomatic heterogeneity of depression and the use of scale total scores as model outcomes. Critically, distinct constellations of symptoms may be related to separable patterns of aberrant neural connectivity and indicate differing likelihoods of responding to a treatment. Thus, total score outcomes are prone to simply agglomerate multiple possibly separable syndromal dimensions and fail to re ect potentially more nuanced and clinically relevant aspects of an individual's symptomatology. The framework we have introduced in previous studies on ECT and SKI [19,20] and applied here with rTMS, is one means by which to contend with symptom heterogeneity and these ndings collectively support that modeling homogenous symptom dimensions may improve prediction of treatment outcomes.

Functional neuroanatomy and clinical implications
Here, we will outline the potential clinical signi cance of the more informative predictors and their related circuitry in predicting antidepressant outcomes following rTMS.

Default Mode Network (DMN)
The DMN is often conceptualized as a network of three subsystems [43]: a midline core encompassing the PCC and anterior mPFC, a dmPFC subsystem which also includes temporal regions [44,45], and a medial temporal lobe system. Activity within the DMN has been related to self-referential thought facilitated by hippocampocortical connectivity [46]. Dysfunctional DMN connectivity has been associated with maladaptive ruminative symptoms [39,40] and emotion dysregulation [47] which, in turn, is predictive of treatment resistance [48]. Echoing this, we observed that elevated pretreatment connectivity of the dmPFC and pCun/PCC, key nodes in two DMN subsystems, predicted treatment resistance.
The dmPFC has been implicated in roles of emotional and behavioral regulation [49][50][51][52]. Given its involvement in these domains, the dmPFC has also been used as a potential TMS target site for a variety of disorders including depression [49,53], obsessive-compulsive disorder [54], and borderline personality disorder [55].
Anterior components of the DMN (aDMN) are hyperconnected in depressed patients [56-58] and this hyperconnectivity may be predictive of antidepressant response following rTMS [59]. An earlier study of ours reported that the aDMN was predictive of CMA improvement following SKI treatment [20]. In the SKI study, however, increased connectivity of the aDMN was associated with greater therapeutic response in CMA symptoms. This discrepancy may relate to differences in treatment mechanisms, in line with recent studies that identi ed differential biological mechanisms of established interventional antidepressants, i.e. ECT and TMS [60]. Alternatively, slight variations in boundaries of aDMN components may partially account for this as a different atlas was used in our previous study.  [70] reported that the PCC's connectivity is elevated in depression while Blum et al.
reported reduced connectivity between the PCC and the caudate in unmedicated patients [71]. The previously mentioned SKI study from our group using a similar predictive framework identi ed that increased pretreatment functional connectivity of the PCC (BA v23ab) predicted less reduction of ruminative symptoms captured by the RRS Re ection subscale [42] which parallels our ndings of increased PCC GC predicting refractory symptoms.

Somatomotor Network (SMN)
Higher pretreatment SMN connectivity also predicted refractory CMA symptoms. A large study of 848 patients with MDD and 794 unaffected controls reported reduced within network connectivity of the SMN and disrupted SMN-DAN connectivity in depressed participants [72]. Another study in adolescents with depression reported similar patterns of hypoconnectivity of the SMN relative to controls [73]. The RSFC of the SMN has also been identi ed in predictive studies. For example, in a cohort of 163 patients with depression receiving escitalopram, sertraline or venlafaxine-XR, greater connectivity between the SMN-DMN as well as the SMN with cinguloopercular network and DAN was predictive of remission after eight weeks of treatment [74]. A study by Leaver et al. also reported that pretreatment RSFC of components of the SMN were predictive of antidepressant response to ECT [75].

Dorsal Attention Network (DAN)
The DAN is involved in external orientation of attention and has also been implicated in the pathophysiology of depression. A large meta-analysis involving 556 patients with depression and 518 controls noted widespread hypoconnectivity between the DAN and the FPCN [41] in patients with depression. A study using both task-based and resting-state functional connectivity identi ed altered low-frequency oscillations of the DAN in unmedicated patients with depression [76]. Disrupted DAN connectivity has been hypothesized to be a signature of attentional bias towards internally-oriented attention (e.g., rumination) at the expense of externally-oriented attention [41]. In contrast with the DMN, we observed that increased pretreatment GC of regions within the DAN predicted more reduced symptoms across the CMA, HDRS-17, and HDRS-6 scales. Taken together with our DMN ndings, we hypothesized that connectivity patterns that were predictive for this set of outcomes may be related to the severity of ruminative symptoms which might also be related to subsequent symptom changes; that is, symptoms of rumination might mediate associations between regional connectivity and symptom reduction.
Exploring this, however, we observed no evidence for these associations.

Visual Network (VN)
Higher pretreatment connectivity of the right fusiform area, a component of the visual network, was predictive of more reduced symptoms following rTMS. Perceptual processing de cits have previously been reported in MDD. For example, eye tracking studies report that patients with depression disproportionately gaze at dysphoric stimuli relative to controls [77] and perceptual processing de cits for non-aversive stimuli have been reported in depressed patients [78]. Several studies have further reported speci c functional connectivity abnormalities of the fusiform area in depression. A study of 62 patients with depression and 61 controls identi ed disrupted connectivity between the fusiform area and sensorimotor regions of the pre-and post-central gyrus [ . In line with this and our own ndings, a recent fMRI study reported that decreased GC of the FPCN with the rest of the brain was associated with more severe depressive symptoms in the general population [87]. Taken together with previous investigations, our ndings may re ect that patients with more intact cognitive control are more responsive to rTMS.

Limitations
There are several limitations to consider. Notably, this is a small cohort and the number of predictive features is much larger than the number of subjects which leaves models more prone to over tting. To combat this, however, we used a conservative 10-repeated 10-fold cross validation approach wherein model performance is evaluated on held-out data and performance is aggregated across repeated cross validations. Furthermore, TMS parameters were not uniform across participants; however, this re ects the natural variation seen in clinical practice.

Conclusions
We used a data-driven approach to identify candidate biomarkers of dimensional symptom changes following rTMS. Importantly, these ndings converge with two earlier studies in electroconvulsive therapy and serial ketamine infusion wherein the degree of change in CMA symptoms was predicted more accurately than broader changes of the HDRS-17 and HDRS-6, highlighting the importance of dimensional measures in therapeutic studies and biomarker development research. Notably, the patterns of pretreatment connectivity informing symptom change are in networks and regions affected by the pathophysiology of depression: anterior and posterior components of the DMN, components of the DAN, VN, and the SMN. We observed that increased connectivity of the DMN and SMN robustly predicted more treatment resistance while increased connectivity of the DAN, FPCN, and VN largely predicted higher therapeutic e cacy. This work may help us to characterize patients who are well-or poorly indicated for rTMS and further highlights the importance of adopting dimensional approaches to modeling outcomes as a viable strategy to accelerate the discovery of biomarkers for treatment strati cation. Figure 2 Partial dependence plots for signi cantly predicted outcomes illustrating the expected change in symptoms (yaxis) across the observed range of the top six imaging predictors (x-axis). We also highlight the Schaefer atlas parcellations with predictive global connectivity measures below each plot. Section (a) highlights partial dependence plots for the core mood/anhedonia (CMA) dimension. From left to right RH_SomMot_19 (right dorsal primary motor cortex); LH_SomMot_15 (left dorsomedial primary somatosensory cortex), RH_Default_PFCdPFCm_4 (right anterior medial PFC), LH_Default_pCunPCC_2 (precuneus/posterior cingulate cortex), RH_Vis_3 (right fusiform area), RH_DorsAttn_Post_2 (right extrastriate cortex). Section (b) partial dependence plots for 17-item Hamilton Depression Rating Scale (HDRS). From left to right, predictors are