Individual interregional perfusion between the left dorsolateral prefrontal cortex stimulation targets and the subgenual anterior cortex predicts response and remission to aiTBS treatment in medication-resistant depression: The influence of behavioral inhibition

BACKGROUND
Accelerated intermittent Theta Burst Stimulation (aiTBS) has been put forward as an effective treatment to alleviate depressive symptoms. Baseline functional connectivity (FC) patterns between the left dorsolateral prefrontal cortex (DLPFC) and the subgenual anterior cortex (sgACC) have gained a lot of attention as a potential biomarker for response. However, arterial spin labeling (ASL) - measuring regional cerebral blood flow - may allow a more straightforward physiological interpretation of such interregional functional connections.


OBJECTIVES
We investigated whether baseline covariance perfusion connectivity between the individually stimulated left DLPFC targets and sgACC could predict meaningful clinical outcome. Considering that individual characteristics may influence efficacy prediction, all patients were also assessed with the Behavioral Inhibition System/Behavioral Activation System (BIS/BAS) scale.


METHODS
After baseline ASL scanning, forty-one medication-resistant depressed patients received twenty sessions of neuronavigated left DLPFC aiTBS in an accelerated sham-controlled crossover fashion, where all stimulation sessions were spread over four days (Trial registration: http://clinicaltrials.gov/show/NCT01832805).


RESULTS
Stronger individual baseline interregional covariance perfusion connectivity patterns predicted response and/or remission. Furthermore, responders and remitters with higher BIS scores displayed stronger baseline interregional perfusion connections.


CONCLUSIONS
Targeting the left DLPFC with aiTBS based on personal structural imaging data only may not be the most optimal method to enhance meaningful antidepressant responses. Individual baseline interregional perfusion connectivity could be an important added brain imaging method for individual optimization of more valid stimulation targets within the left DLPFC. Additional therapies dealing with behavioral inhibition may be warranted.


Introduction
Current European and North American clinical guidelines put repetitive Transcranial Magnetic Stimulation (rTMS) forward as an approved treatment option for patients with major depression [1]. Given the relatively long duration of current rTMS treatment courses, with daily single sessions spread over four to six weeks, accelerated (a)rTMS protocols have been introduced recently. A reduction in the number of clinical visits and a faster time to respond can only be beneficial for the patient's comfort [2]. Encouragingly, classical daily rTMS and arTMS protocols show similar clinical efficacy [3], indicating that the major advantage of accelerated stimulation is the shorter treatment duration and a likely faster onset of clinical response. In addition, delayed clinical responses are possible, and depending on specific clinical symptoms, patients may follow distinct treatment response trajectories [4,5]. Notwithstanding the potential time gain of such accelerated rTMS protocols, unfortunately, not all patients will benefit. Recently, a new form of rTMS has been introduced -theta-burst stimulation (TBS) -thought to produce similar outcome effects in comparison to the standard rTMS protocols. Because the stimulation pulses are provided in bursts, the duration of a single TBS session is significantly shortened. In addition, it is also assumed that accelerated TBS treatment protocols may show similar clinical efficacy when compared to accelerated rTMS in treatment-resistant depression [2].
Because individual functional and structural characteristics may influence efficacy, to personalize rTMS treatment parameters, applying brain imaging as a predictor tool gains momentum (for a review see Ref. [6]). Indeed, current research predominantly focuses on the predictive value of individual functional connectivity (FC) patterns between the stimulated area (mostly the left dorsolateral prefrontal cortex (DLPFC) and (regions within) the subgenual anterior cingulate cortex (sgACC). This in addition to a more accurate structural localization of the optimal target zone based on individual anatomical information. In a landmark paper of Fox et al. [7], these authors found that the antidepressant efficacy of daily rTMS treatment was associated with baseline resting state (rs)FC anticorrelation between the more anteriorly located left DLPFC rTMS target sites and the sgACC seed area. Several other studies, applying daily rTMS [8], as well as accelerated stimulation [9,10], showed -but not always -significant similar (anti)correlations in those depressed patients who benefited most from this kind of treatment. Methodological differences in rsFC assessment and analysis may explain these discrepancies to some extent. Arterial spin labelling (ASL) may provide a more straightforward way to determine the predictive value of the individually targeted left DLPFC and functionally connected sgACC in response to rTMS treatment. Considered as a reliable physiological marker of neuronal activity [11e13], ASL is a noninvasive fMRI technique which does not use a radioactive agent but rather arterial water as an endogenous tracer providing stable absolute quantification of regional cerebral blood flow (rCBF) [14]. ASL is quantitative, stable over time, and less variable across subjects [15], which makes it a specific, useful noninvasive method to measure rCBF in clinical studies [16]. Earlier research, applying covariance analytical approaches on structural anatomical brain data, has shown that baseline cortical morphological abnormality patterns between different regions of interest (ROIs) can increase clinical outcome prediction accuracy on an individual basis [17,18]. Given the more stable rCBF quantifications measured with ASL, the application of covariance analysis between the individual targeted areas within the left DLPFC and the sgACC may be in particular suitable to predict treatment efficacy.
Furthermore, former research has suggested that the optimal left DLPFC stimulation targets may comprise symptom-specific clusters, defined as "dysphoric" and "anxiosomatic" [19]. Even more, advanced statistics and machine learning approaches have been used to finetune symptom clustering and to define subtypes of MDD, to increase antidepressant efficacy in rTMS treatment (for a recent review, see Ref. [6]). Based on the clustering of anxiety and anhedonia dimensions and associated resting state fMRI connectivity patterns, Drysdale and colleagues [20] identified two biotypes to be more responsive to rTMS treatment when applied over the dorsomedial prefrontal cortex (DMPFC). This study suggests that clinical response patterns to rTMS treatment could be associated with certain depression subtypes, which differ in behavioral inhibition (anxiety) and/or behavioral (des)activation (reward) patterns. In this regard, introduced by Carver and White [21], the Behavioral Inhibition System (BIS) and Behavioral Approach/Activation System (BAS) -respectively assumed to measure avoidance/anxiety and reward/impulsivity -could be an important asset to such brain imaging FC measurements. Furthermore, both BIS/BAS constructs not only have been implicated in the pathophysiology of mood disorders, but they have also been evaluated as potential outcome predictors [22]. However, in clinical (accelerated) rTMS research, no studies have examined such associations with BIS/BAS. Given the assumption that rTMS treatment outcomes in depression can depend on behavioral characteristics, it would be of high interest to verify whether traits related to inhibition and reward would influence these connectivity patterns.
Consequently, we hypothesized that clinical response and/or remission can be predicted by baseline covariance perfusion connectivity between the individually stimulated left DLPFC targets and the sgACC. Given that both BIS/BAS (sub)scales have been linked to clinical outcomes, we also explored whether the predictive value of this interregional covariance connectivity was influenced by individual scores on this self-report measure.

Participants
This randomized double-blind sham-controlled crossover study was approved by the Ghent University Hospital ethics committee. The study was registered in the Clinical Trials.gov database (http:// clinicaltrials.gov/show/NCT01832805) and is part of a larger project investigating the influence of aiTBS on various neurocognitive and neurobiological markers. Out of fifty enrolled patients, both ASL scan data and individual MNI coordinates for the left DLPFC were available for forty-one patients (thirty-one females; ten males; mean age ¼ 40.20 years, SD ¼ 11.83 years; mean baseline 17-item Hamilton Depression Rating Scale (HDRS [23]; Mean ¼ 20.98, SD ¼ 5.59). That relatively more females than males were included has to be considered as coincidence. The full behavioral data were published in Duprat et al. [4].
Inclusion criteria for the study were 1) a major depressive episode selected with the Mini-International Neuropsychiatric Interview (MINI [24], 2) right-handedness, 3) at least stage I treatment resistance according to the Rush et al. [25] criteria, 4) and at least two weeks free from psychotropic agents. Only habitual benzodiazepine agents were allowed. The maximum allowed dose of benzodiazepines was the equivalent of 40 mg diazepam. Any changes during the stimulation sessions would have resulted in drop-out from the study. No patients dropped out according to this criterion. See also Duprat et al. [4].
Severity of depressive symptoms was assessed with the 17-item HDRS by a certified psychiatrist, not related to the study and blinded to the actual treatment of the patient. All patients gave written informed consent.

Behavioral Inhibition System/Behavioral Activation System scales (BIS/BAS)
The BIS/BAS Scale is a 20-item self-report measure (with four additional filler questions) that assesses different types of reinforcement of two biological motivational systems -activation and inhibition -respectively regulating approach and withdrawal behaviors in response to environmental cues [21]. Participants are asked to respond to every item using a one to four scale, ranging from one (quite untrue of you) to four (quite true of you). The scale comprises four subscales: one BIS scale and three BAS subscales (Drive, Reward Responsiveness, and Fun Seeking). Whereas the unidimensional BIS scale is thought to measure the degree of negative affect experienced in response to threatening cues, the BAS Reward Responsiveness scale focuses on affective responding, the BAS Drive scale exclusively assesses behavioral responding, and the BAS Fun Seeking scale assesses for both affective and behavioral responding [22].

Experimental setup
After a washout period, all patients were at least two weeks antidepressant (AD) free. After the baseline ASL (ASL baseline ) scan, patients were randomized to receive 20 sessions of active or sham aiTBS treatment. Following a crossover pattern the next week, strictly following the same treatment schedule but with a change of stimulation. Depression severity symptoms were assessed at baseline and each time three days after the last run of stimulation. To detect delayed aiTBS clinical effects, with a delay of exactly two weeks, the 17-item HDRS was assessed for the last time. See also Supplemental Fig. S1 for more details.
Accelerated intermittent TBS stimulation was applied using a Magstim Rapid2 Plus1 magnetic stimulator (Magstim Company Limited, Minneapolis, USA) with a 70 mm Double Air Film figureof-eight shaped cooled coil. The sham coil is identical in all aspects to its active variant, but without stimulation output. For every individual, we visually located the left DLPFC on the 3D surface rendering of the brain based on the known gyral morphology and marked the center part of the mid prefrontal gyrus as the left DLPFC target, comprising Brodmann area 9/46 (See also Figs. 2 and 3). To accurately target this part of the left DLPFC, we used Brainsight neuronavigation (Brainsight™, Rogue Resolutions, Inc). The individual native space coordinates of the targeted DLPFC areas were saved in the Brainsight module and were normalized into personal MNI coordinates with the forward deformation field produced by SPM12 structural segmentation.
The 20 iTBS sessions were spread over four days at five sessions per day, amounting to a total of 32.400 stimuli. Patients received 1620 pulses per session in 54 triplet bursts with a train duration of 2 s and an intertrain interval of 8 s (from start to end and including the 2-s train duration). Between the two sessions, there was a pause of approximately 15 min. Single-pulse TMS in combination with motor evoked potentials (MEP) -measured with surface electromyography (EMG) registration -determined the optimal right abductor pollicis brevis muscle (APB) response at rest (rMT). Positive MEP responses of at least 50 mV (peak-to-peak amplitude) had to be produced in at least five of 10 consecutive trials before the selected part of the cortex was accepted as the motor cortex related to the contralateral APB muscle. A stimulation intensity of 110% of the subject's rMT was used during the entire aiTBS experiment. All patients reached the full dosage of aiTBS treatment.
Multi-delay pulsed arterial spin labeled (pASL) images with a 3D GRASE readout were obtained with the following parameters: TR ¼ 3.4s, TE ¼ 14.46 ms, labelling duration ¼ 1400 ms, postlabelling delay changing from 300 to 3000 ms in steps of 300 ms, resulting in twelve pairs of slice-selective (SS) and nonselective (NS) images, scan duration ¼ 5.26 min. Patients were asked to stay awake with their eyes closed. All confirmed they had stayed awake during scanning.
To compare among different subjects, the individual left DLPFC localization stimulation coordinates, representing the patientspecific stimulation site, were saved in the BrainSight neuronavigation system and all converted into Montreal Neurological Institute (MNI) coordinate space, using SPM12 (WellcomeTrust Center for Neuroimaging, London, UK). Based on the left DLPFC template, following the method of Fox and colleagues [26], we defined the individual DLPFC stimulation targets (MNI coordinates) as spheres with a 15-mm radius as shown in Supplemental Material Fig. S2.
Similar to Cole et al. [10], and our former accelerated HF-rTMS treatment study combined with 18 FDG-PET [27], we used the WFU_pickatlas to define the sgACC region of interest. Given that the sgACC seed selection by Fox et al. [7] and our own sgACC seed used in our former FC studies (e.g., Ref. [9]) were located within the right sgACC, for this ASL study, we only focused on the right sgACC anatomical region. According to the Rosen et al. [28] paper, using the Yeo et al. [29] network atlas, we also verified whether the individual left DLPFC targeted spots were related to different functional networks; for response/remission and nonresponse.

Behavioral data
All behavioral data needed for our hypothesis were analyzed using SPSS 26 (IBM, Statistical Package for the Social Sciences, Chicago). The significance level was set at p < 0.05, two-tailed, for all analyses. Where necessary, we applied the Greenhouse-Geisser correction to ensure the assumption of sphericity. For this ASL study, we focused on a meaningful clinical outcome, defined as response (50% reduction of the baseline HDRS score), and/or remission (having a score on the HDRS 7) at any given time point, in the text further abbreviated as RR (Response and/or Remission). Patients with no response and/or remission at any of these time points were considered as non-responders, abbreviated further as NR. Given that (accelerated) rTMS treatment paradigms might be prone to placebo responses [2], those patients having had sham aiTBS and yielded RR in the first week, but not any further, were also considered NR. See also Supplemental Fig. 1.

Imaging data
The pASL and anatomical images were preprocessed and analyzed using FSL (FMRIB, Oxford, UK) and SPM12. The individual structural images were segmented into grey matter, white matter, and cerebrospinal fluid with fsl_anat. All NS and SS images were first corrected for motion and registered to the anatomical image. Then, twelve perfusion-weighted images were generated by surround subtraction, i.e., the differences between the paired SS and NS images. The perfusion-weighted images were submitted for CBF estimation using the Oxford ASL (oxford_asl) tool in FSL. The partial volume effects in the generated CBF maps were corrected by a regression algorithm in the PETPVE12 toolbox (https://github.com/ GGonEsc/petpve12). Global mean normalization was applied to the CBF data. Finally, the generated CBF maps were spatially normalized into MNI space and smoothed with a 6 mm full-width halfmaximum Gaussian kernel.

Individual interregional covariance connectivity
Baseline CBF map was used to calculate the interregional covariance connectivity estimated with a previously proposed measure of similarity [17]. The similarity (s) between the individually targeted left DLPFC and the sgACC was computed as: where m k and s k denote the mean and standard deviation of the regional CBF, respectively.

Statistical analysis
To compare the covariance connectivity strength between two groups (RR and NR), we applied a nonparametric permutation test to compute the statistical significance of the between-group difference, with age, gender, and physical distance as covariates. For the latter, the physical connection distance was estimated as the mean Euclidean distance between the individual stimulus site and each voxel in the (right) sgACC anatomical area. The physical connection distance has been shown to be a central parameter for efficient information propagation [30]. Furthermore, we calculated the actual T-value and the T-values of 10000 randomly mixed data splits. The p-value of the actual difference is then estimated from this null distribution.
To assess the relationship between the individual interregional covariance connectivity strength, Spearman partial correlation analysis, with age and gender as covariates, was further performed for each BIS/BAS subscale separately. To correct for multiple comparisons, the significance of p < 0.05 was Bonferroni corrected for four BIS/BAS subscales and two conditions: RR and NR.

Behavioral results
Twenty medication-resistant depressed (MRD) patients received active iTBS treatment in the first week and sham stimulation in the second week. Twenty-one patients followed the reverse order. There were no significant order differences in gender (c 2 (40) From the forty-one included MRD patients, at least at one of the predefined clinical assessment points twenty-two patients experienced a response, and fourteen of them reached remission. This means that according to our definition of RR, in total twenty-three MRD patients met this criterion at least at one of these assessments. Only one patient responded to the first week of sham aiTBS, but after active treatment she went further in remission, so she was also considered RR. For a detailed overview, see Supplemental Table 1  For the BIS/BAS scales, three data sets were incomplete or went missing. Consequently, these analyses were performed for thirtyeight MRD patients. Baseline BIS/BAS scores, corrected for age and gender, were not significantly different between RR and NR (BIS, t(34) ¼ -1.87, p ¼ 0.07; BAS drive, t(34) ¼ 1. 35 Fig. 3.
On the assumption that the left DLPFC target locations would differ between RR and NR -in relation to the connected functional networks [20] -according to the Yeo et al. [21] atlas nearly all our left DLPFC targets were located within the default mode network (DMN), regardless of clinical outcome. For more details see Figs. 2 and 3, and Supplemental Table 3.

Influence of BIS/BAS on individual baseline interregional covariance connectivity
Spearman (partial) correlation analysis showed that only in the RR group the BIS scores were positively associated with the individual interregional covariance connectivity (r ¼ 0.6058, p ¼ 0.006). This correlation was significantly stronger in RR when compared to the NR group (z ¼ 2.6651, nonparametric p ¼ 0.0116). No significant correlation between any of the BAS subscale scores and the individual interregional covariance connectivity was observed (p's > 0.05). See also Fig. 4 and Supplemental Table 2.

Discussion
By using baseline perfusion (pASL), we retrospectively examined whether the functional connections between the individually targeted left DLPFC coordinates and the (right) sgACC could predict clinical efficacy in a well-defined ADfree cohort of MRD patients following aiTBS treatment. We found that stronger baseline interregional covariance perfusion patterns were predictive for MRD patients experiencing a meaningful clinical response (RR) at least at one of the predefined assessment moments. Although the individually determined physical connection distances were not significantly different between RR and NR, the individual stimulation targets for RR were more widely spread around the structurally predefined left DLPFC spot, whereas for NR the targets were more closely concentrated around this region. Importantly, the Euclidean distance from the individual stimulus site to their mean coordinates was significantly different between groups. This implies that the closer to the predefined structural stimulation area, the weaker the interregional covariance perfusion and the less likely MRD patients will respond. This also suggests that when targeting the left DLPFC based on individual structural cortical data only, this may result in suboptimal clinical efficacy.
Here, it has to be noted that our group mean targeted MNI coordinates are close to the mean targeted area when using the 'EEG F3' method (See Ref. [31], and Supplemental Fig. 4), a region found to be less clinically effective than the more anteriorly located left DLPFC targets as proposed by Fox and colleagues [7]. Notwithstanding that it has been postulated that a more correct structural targeting of the DLPFC may be necessary to increase response rates (see also [6]), our current findings indicate that when only including individual anatomical cortical information this might not be sufficient to significantly increase clinical outcome, at least for aiTBS in MRD. This observation also agrees with the few clinical studies in depression available unable to find clear-cut clinical superiority when comparing neuronavigated to non-neuronavigated rTMS treatment [32]. Of note, by definition (see formula (1)), our interregional covariance perfusion patterns always show a positive connectivity. Because stronger anticorrelations also represent stronger connectivity, in essence, this does not really disagree with the expected anticorrelation patterns in relation to clinical efficacy, 1 Of note, by adding the temporal signal-to-noise (SNR) ratio as covariate into the analysis, this did not change the outcome result, t ¼ 2.121, nonparametric pvalue ¼ 0.022. as often is reported in resting state FC research [7,8,10]. However, because such temporal dynamics were not part of the pASL sequence, we cannot further interpret the directions of the interregional covariance perfusion patterns.
Considering the influence of the BIS/BAS scales on our outcome measurements, first, none of them were significantly different between RR and NR. This can partly be explained by the inclusion of a well-defined sample of unipolar MRD patients excluding psychosis, substance dependency, acute suicidality, and nonsignificant baseline depression severity symptom scores. Furthermore, the presence of comorbid psychiatric illnesses, such as generalized anxiety disorder or social phobia were not significantly different between RR and NR (see Supplemental Table 1). When analyzing the effects of interregional covariance connectivity, only those MRD patients responsive to aiTBS treatment showed significantly stronger connectivity correlations with the individual BIS scores. This suggests that when targeting the individual left DLPFC targets, in addition to a stronger connectivity strength with the (right) sgACC, that individuals more sensitive to threat and punishment could be more susceptible to this kind of treatment. On the other hand, as it has been documented that in depressed patients higher BIS scores were positively associated between maladaptive emotion regulation patterns -including rumination, catastrophizing, and self-blame [33] -our findings may also indicate that when optimally stimulated, the effects of aiTBS may result in a more adaptive regulation of depressive emotional processes. This assumption does not disagree with our former observations where stronger sgACC FC in parts of the medial orbitofrontal cortex were associated with prompt attenuation of negative thinking (hopelessness) [9]. Nonetheless, the contribution of BIS/BAS scales as an influential mediator on the examined interregional covariance connectivity should interpret carefully, as the original underlying key neuroanatomical regions and (serotonergic) pathways for BIS comprise the hippocampus and amygdala [34], and not directly the two brain regions examined in our study. On the other hand, our connectivity observations may also imply an underlying vascular mechanism which directly drives the neural activity underlying the aiTBS treatment response based on reward responsiveness, which is considered related to the dopaminergic system [35,36]. Of note, our neuronavigated left DLPFC stimulation target closely corresponds to the medial and posterior located 'anxiosomatic' cluster target as defined by Siddiqi et al. [19], which, besides insomnia and sexual dysfunction, comprise feelings of failure and indecisiveness, related to behavioral inhibition. Following the reasoning of Rosen and colleagues [28], it is possible that in RR the DMN was more successfully targeted via the individual left DLPFC in MRD patients with a more 'anxiosomatic' profile. It has been suggested that a hyperconnectivity between the DMN and the sgACC may mediate interactions with the 'affective' salience network (SN) via "affective-laden behavioral withdrawal", integrating self-referential DMN-mediated processes with behavioral manifestations of depression [37]. Furthermore, the depressed state is associated with reduced excitatory connectivity mainly within the DMN, and between the default mode and SN [38]. Proper DMN functioning depends on a good controlling SN, modulating cognitive control and self-regulation, appropriate behavior, and adequate emotional responses, which becomes dysfunctional when clinically depressed [39]. Of interest, it has been documented that rTMS treatment can modulate this SN in regaining cognitive control when being effective [40,41]. However, given that we did not include neuronal network analysis, we can only speculate that our observed interregional perfusion covariance connectivity findings in relation to behavioral inhibition may also be attributed to a more preserved integrity of the SN.
Although our study has a couple of major advantages -such as a well-defined antidepressant-free treatment-resistant depressed patient sample -besides the modest sample size, some limitations should be discussed. All interpretations should be limited to accelerated theta burst stimulation and should not be extrapolated to classical daily rTMS paradigms. Because we aimed at meaningful categorical clinical outcomes (RR), we did not use HDRS scores as a continuous variable. Although we did not observe order effects, the final clinical evaluation after active aTBS was not the same for those receiving sham first. Therefore, our crossover design might not have been the optimal setup. Furthermore, despite taking great care in blinding the patients (including the use of a sham coil that was allegedly identical to the one used for the active stimulation, producing scalp sensations and sound when stimulated), sham undeniably differs from the active condition which may become especially obvious in crossover designs. Because our latest clinical assessment was performed only two weeks after the last stimulation session, delayed RR could have remained undetected. Because temporal dynamics were not considered in our pASL sequence, future work could consider pseudo-continuous ASL as an alternative to the BOLD-fMRI to analyze FC with perfusion time series, which enable us to compare with the anticorrelation patterns as reported in BOLD-fMRI studies. Concerning the sgACC, we focused solely on the right anatomical region, and left or other right sgACC subregions were not examined. This is not a trivial limitation, given that laterality differences in sgACC structure and function have been reported and given that the stimulation localization is mostly over the left DLPFC [42]. It has to be noted that the BIS/BAS scales were only used to evaluate its modulation impact on individual interregional covariance connectivity and not to differentiate between depression subtypes, and of course, other personality scales could also have provided extra information. Lastly, although benzodiazepine intake did not differ between RR and NR, recent studies have demonstrated that benzodiazepines are a negative prognostic factor for rTMS clinical outcome (e.g., Refs. [5,43]).
In conclusion, targeting the left DLPFC based on individual structural data only may not be sufficient to optimize clinical efficacy. The combination of individual baseline perfusion interregional covariance connectivity measurements could be necessary to improve outcome rates. However, to truly challenge the hypothesis that when targeting the (stronger) functional perfusion connections between the left DLPFC and the (right) sgACC will yield higher response and remission rates, larger prospective aiTBS studies are needed. Given that individual baseline left DLPFC-sgACC functional connections were influenced by behavioral inhibition, our brain perfusion findings may be especially relevant to combine aiTBS with psychotherapeutic interventions focusing on maladaptive emotion regulation patterns such as depressive rumination [44].

Declaration of competing interest
None.