Focal application of accelerated iTBS results in global changes in graph measures

Abstract Graph analysis was used to study the effects of accelerated intermittent theta burst stimulation (aiTBS) on the brain's network topology in medication‐resistant depressed patients. Anatomical and resting‐state functional MRI (rs‐fMRI) was recorded at baseline and after sham and verum stimulation. Depression severity was assessed using the Hamilton Depression Rating Scale (HDRS). Using various graph measures, the different effects of sham and verum aiTBS were calculated. It was also investigated whether changes in graph measures were correlated to clinical responses. Furthermore, by correlating baseline graph measures with the changes in HDRS in terms of percentage, the potential of graph measures as biomarker was studied. Although no differences were observed between the effects of verum and sham stimulation on whole‐brain graph measures and changes in graph measures did not correlate with clinical response, the baseline values of clustering coefficient and global efficiency showed to be predictive of the clinical response to verum aiTBS. Nodal effects were found throughout the whole brain. The distribution of these effects could not be linked to the strength of the functional connectivity between the stimulation site and the node. This study showed that the effects of aiTBS on graph measures distribute beyond the actual stimulation site. However, additional research into the complex interactions between different areas in the brain is necessary to understand the effects of aiTBS in more detail.

frequency rTMS, delivering pulses at a frequency higher than 5 Hz, is currently FDA approved as treatment for patients with medication resistant major depressive disorder (MDD), which is approximately one-third of all MDD patients. Left prefrontal high-frequency rTMS has shown to be an effective and safe treatment in adult MDD patients documented as medication resistant (Pascual-Leone, Rubio, Pallardó, & Catalá, 1996;George et al. 2010;George, Taylor, & Baron Short, 2013;Padberg & George, 2009;Baeken et al., 2013;Loo, McFarquhar, & Mitchell, 2008). The rationale to stimulate these parts of the cortex is based on earlier studies showing clear involvement of the prefrontal cortex (PFC) in the pathophysiology of MDD (Koenigs & Grafman, 2009). More specifically, the ventromedial PFC (VMPFC) shows hyperactivity, whereas the dorsolateral PFC (DLPFC) shows hypoactivity, as demonstrated by multiple imaging studies (Mulders, van Eijndhoven, Schene, Beckmann, & Tendolkar, 2015). Reversing these effectsdecreasing the activity of the VMPFC or increasing the activity of the DLPFC-has been proposed as a possible mechanism by which rTMS treatment can achieve response and remission from depressive symptoms (George, 2010;Seminowicz et al., 2004).
Standard rTMS guidelines to treat depression follow mostly a daily pattern, with applied frequencies from 1 to 20 Hz, repeated for 4-6 weeks (Perera et al., 2016). With such protocols, clinical effectiveness remains however rather modest. To improve clinical outcome, new treatment parameters are currently under investigation.
One new approach is accelerated rTMS, where a similar amount of stimulation sessions is concentrated over a couple of days instead of the more conventional daily sessions, spread over multiple weeks.
Another line of research focuses on theta burst stimulation (TBS) (Huang, Edwards, Rounis, Bhatia, & Rothwell, 2005), where a particular set of parameter deliverables applies bursts of 3 stimuli at 50 Hz and is repeated every 200 ms (5 Hz, theta range). TBS has shown comparable clinical efficacy compared to rTMS but stimuli are delivered during a shorter period and usually with a lower intensity (Blumberger, Vila-rodriguez, Knyahnytska, et al., 2018). Intermittent TBS (iTBS), the administration of 2 s of TBS alternated with 8 s rest, has been investigated for treatment of MDD (Bakker et al., 2015;Chistyakov, Rubicsek, Kaplan, Zaaroor, & Klein, 2010;Li et al., 2014), based on the excitatory character of the standard iTBS protocol (600 stimuli at 80% active motor threshold) (Huang et al., 2005).
To maximize clinical efficacy within a shorter time period, an intensive accelerated iTBS (aiTBS) protocol, consisting of multiple iTBS sessions per day, was recently tested as possible treatment for depression in our group. Duprat et al. (2016) showed a rapid significant decrease in depression severity symptoms assessed with the 17-item Hamilton Depression Rating Scale (HDRS) (Hamilton, 1967) after 4 stimulation days. Although clinical effects were found both after sham and verum aiTBS, the most meaningful clinical outcomes regarding response and remission were observed 2 weeks after the aiTBS protocol, during follow-up. While only 28% of the patients showed a 50% reduction of their initial HDRS score at the end of the stimulation procedure, response rates mounted up to 38% 2 weeks later, indicating delayed clinical effects. Furthermore, 30% of the responders were considered in clinical remission.
How aiTBS has the potential to improve depression symptoms over such a limited period in medication resistant MDD patients remains to be elucidated. Because it is known that the effects of stimulation are propagated through the brain via anatomical and functional connections (Amico et al., 2017;Fox et al., 2014), the effect of aiTBS might occur on a network level. In this study, the effect of this aiTBS protocol on the brain's network topology is investigated by means of graph analysis derived from resting-state functional MRI (rs-fMRI) data of a group of MDD patients. Graph analysis is a mathematical concept to quantify networks, for example, brain networks, according to various neurobiologically meaningful properties such as integration and segregation (Bortoletto, Veniero, Thut, & Miniussi, 2015;Rubinov & Sporns, 2010). Combining rs-fMRI datasets before and after a brain stimulation protocol with graph analysis allows one to map the network changes throughout the whole brain induced by TMS, instead of just looking at single connections at a time, as it is done in many functional connectivity studies.
Previous studies have investigated the brain's network topology in patients with MDD. Graph analyses were performed based on cortical thickness (Mak, Colloby, Thomas, & O'Brien, 2016), voxel based morphometry measures (Lim, Jung, & Aizenstein, 2013), structural connectivity using diffusion MRI data (Ajilore et al., 2014;Chen et al., 2016;Korgaonkar, Fornito, Williams, & Grieve, 2014;Singh et al., 2013), or functional connections using rs-fMRI datasets (Bohr et al., 2013;Li et al., 2016). The reported differences in graph measures between healthy volunteers and MDD patients were ambiguous. On one hand, some studies did not find differences, and on the other hand, increases in clustering coefficient, local efficiency, and path lengths were reported.
To study the effects of stimulation on network level, only few studies have been performed combining brain stimulation and graph theory: for example, Shafi et al. (2014) and Deng et al. (2015) used resting EEG data to examine the effects of continuous TBS and rTMS respectively. Shafi et al. (2014) showed frequency band dependent effects of stimulation on clustering coefficient and local efficiency: the beta band showed increases in clustering coefficient after cTBS, whereas alpha band showed decreases in clustering coefficient along with increased path length. Deng et al. (2015) showed reduced smallworldness in the beta frequency band after stimulation. Vecchio et al. (2018) performed source localization on EEG data recorded before and after transcranial direct current stimulation (tDCS) and showed that anodal tDCS over the motor cortex reduces small-worldness. Park et al. (2014), Polanía et al. (2011), andCocchi et al. (2015) studied the effects of various stimulation techniques using task fMRI and rs-fMRI data. Park observed a correlation between the motor performance change and the increase and decrease in global and local efficiency respectively, induced by 10 Hz rTMS (Park et al., 2014).

Cocchi showed different effects of continuous versus inhibitory TBS
represented by modularity, out-degree participation index, and within-module degree (Cocchi et al., 2015). Polanía et al. (2011) combined anodal tDCS over the motor cortex with rs-fMRI derived graph measures and found increases in path length in the somatomotor areas after stimulation.
Besides the effect of aiTBS on graph measures, it will clinically be relevant to investigate if graph measures can be used as biomarkers to predict the outcome of this stimulation protocol. Previously, it has been shown that rs-fMRI connectivities can be used for this purpose. Drysdale et al. (2016) derived four depression subtypes that seem to respond differently to rTMS treatment. And Fox et al. (2012 demonstrated that the clinical effects of rTMS are linked to the functional anti-correlation between the subgenual anterior cingulate cortex (sgACC) and the stimulation spot in the left DLPFC. This anticorrelation between the sgACC and parts of the left superior medial prefrontal cortex was also suggested to have predictive value for the outcome of accelerated rTMS in a cohort of MDD patients (Baeken et al., 2014), although in another accelerated iTBS this was not found to be that straightforward (Baeken, Duprat, Wu, De Raedt, & van Heeringen, 2017a). Nevertheless, Downar et al. (2014) showed in a cohort of MDD patients that the graph measure betweenness centrality can be used to distinguish responders from nonresponders to rTMS to the dorsomedial prefrontal cortex.
Specifically, this is the first study using graph analysis to investigate the clinical effects of the relatively new aiTBS treatment protocol. Graph analysis was performed on the whole-brain level, using the clustering coefficient, global efficiency, small-worldness, and modularity, and on the nodal level, using the degree, and the betweenness centrality as graph measures. Due to the presumably excitatory character of iTBS, we hypothesized that aiTBS would increase all four whole-brain graph measures. On nodal level, we expected to find mostly increases in degree and betweenness centrality in nodes related to the pathophysiology of MDD. Furthermore, we expected that changes in graph measures would be linked to the clinical response. We also hypothesized that changes in functional connectivity, expressed by graph measures, would not only occur in the stimulated area (the left DLPFC), but will also be present in functionally connected regions.

| METHODS
This study (http://clinicaltrials.gov/show/NCT01832805) was approved by the local Ghent University Hospital ethics committee and is in accordance with the declaration of Helsinki (2004). All patients gave written informed consent.

| Inclusion criteria
Fifty right-handed MDD patients were included in this study. MDD was diagnosed using the structured Mini-International Neuropsychiatric Interview (MINI; Sheehan et al., 1998). All patients were at least stage I treatment resistant according to the Rush criteria (Rush, Thase, & Dube, 2003). They had a minimum of one unsuccessful treatment trial with selective serotonin reuptake inhibitors/serotonin and norepinephrine reuptake inhibitors (SSRI/SNRI). Medication was tapered off before the aiTBS treatment period, so all were medication-free for at least 2 weeks before the start of the first stimulation session. More extensive information about the patients and clinical outcome can be found in Duprat et al. (2016).
During the resting-state measurement, patients were asked to stay awake with their eyes closed. On Days 2-5 and Days 9-12, verum or sham aiTBS was applied depending on the randomization order. A Magstim Rapid 2 Plus 1 magnetic stimulator (Magstim Company Limited, Wales, UK) connected to a verum or sham figure-of-eight shaped coil (Magstim 70 mm double air film [sham] coil) was used to apply the verum and sham stimulation respectively. On the 8th day (T2) and on the 15th day (T3), so 3 days after the stimulation, the imaging protocol was repeated. At the same days when imaging was performed (T1, T2, and T3) and additionally 2 weeks after the last stimulation (T4), depression severity symptoms were assessed using the 17-item HDRS questionnaire (Hamilton, 1967).
Before the first stimulation session, the resting motor threshold (rMT) was determined based on motor evoked potentials (MEPs) induced in the right abductor pollicis brevis (APB) after applying single pulses to the hotspot. During four consecutive days, five daily sessions of iTBS were applied at 110% rMT to the left DLPFC: the center part of the midprefrontal gyrus (Brodmann area 9/46) based on structural MRI of each individual (Peleman et al., 2010). Positioning of the coil was maintained with the BrainSight neuronavigation system (Brainsight™, Rogue Resolutions, Inc). One iTBS session consisted of 54 trains of 10 bursts of 3 stimuli. Two seconds of stimulation were given in an 8 second cycling period. This adds up to 1,620 stimuli per session with a total number of 32,400 stimuli during the four-day treatment. There were breaks of~15 min between the stimulation sessions. During the stimulation, patients were blindfolded, wore earplugs, and were kept unaware of the type of stimulation (sham or verum) they received.

| Graph analysis
Functional connectivity analyses were performed using the rs-fMRI data from T1 and T2. In this first week of the study design, patients received either sham or verum aiTBS depending on the order of randomization. The second part of the study protocol, the period between T2 and T3 after the cross-over, was not used to be able to study the pure effects of sham and verum aiTBS. The duration of the after-effect of 4 days aiTBS is not yet known and as there was only a weekend between the stimulation weeks, effects of verum and sham might be crossed over into the second week.

Data were preprocessed with MATLAB 2015b (The Mathworks
Inc., Natrick, MA) and SPM12 (Wellcome Trust Centre for Neuroimaging, London, UK) according to standard steps. After realignment, volumes with excessive motion, quantified as >0.3 mm framewise displacement, were discarded for further analysis. Complete datasets were excluded if more than 100 volumes had to be removed (see Appendix A, Figure A1). Six motion regressors, and additionally a white-matter and cerebrospinal fluid regressor were used to correct the data using SPM's REST toolbox (Song et al., 2011). The latter two regressors were defined as the mean of the time-series within the eroded white-matter and cerebrospinal fluid masks, respectively. Temporal bandpass filtering was applied with cutoff frequencies of 0.1 and 0.01 Hz. The brain datasets were parcellated using the parcellation scheme from Drysdale et al. (2016), using the 264 parcels, further referred to as nodes, from Power et al. (2011) and additionally 13 subcortical gray matter structures (see Appendix A, Table A1 for additional information). For all nodes, the mean time-series was computed by averaging all the voxel time-series belonging to that node. The temporal signalto-noise ratio (tSNR) criterion was used to remove nodes with unreliable time-series from further analyses. Nodes were discarded if tSNR <40 in more than 10% of the datasets (Liston et al., 2014). Furthermore, if more than 10% of the nodes within one dataset had tSNR <40, the dataset (both T1 and T2) were removed from further analysis (Appendix Figure A2).
For every patient and for both time points (T1 and T2), a connectivity matrix was calculated as the Pearson correlation between all the node time-series, herewith rejecting the first 10 volumes to ensure scanner stability. The connections in this connectivity matrix are further referred to as edges. All edges are scaled to be in the range between zero and one (Schwarz & McGonigle, 2011) in a three-step process. First, the range of the connectivity matrix was defined by subtracting the minimum value from the maximum value. Second, all edge-values were divided by the range. Last, the minimum value of the new matrix was added which results in a scaled matrix between zero and one. This method was repeated for every subject and for every time-point separately.
The MATLAB-based Brain Connectivity Toolbox (Rubinov & Sporns, 2010) and the Graph Analysis Toolbox (Hosseini, Hadi, Hoeft, & Kesler, 2012) were used to calculate graph measures that quantify the brain's network organization (Bullmore & Sporns, 2009Rubinov & Sporns, 2010). On whole-brain level, four weighted graph parameters were calculated from every connectivity matrix: clustering coefficient, global efficiency, small-worldness, and modularity. Here, high clustering coefficients are associated with high local efficiency regarding information transfer and robustness (Bullmore & Sporns, 2009). The modularity measure represents the way in which a network can be subdivided into modules: groups of nodes with a high number of within-group links and a low number of between-group connections (Girvan & Newman, 2002;Newman, 2004). Functional integration can be described by path length and efficiencies. High functional connectivity values can be translated to short path lengths and high efficiencies. The path length is the average of the shortest routes of information flow between pairs of nodes. Global efficiency can be calculated by inverting the path lengths. Moreover, the smallworldness was calculated. Small-world networks are assumed to be efficient, both locally and globally (Rubinov & Sporns, 2010). To calculate the small-worldness, the clustering coefficient and path length were normalized by dividing them by their equivalents derived from random networks. Random networks were obtained using 20 randomization steps, leaving the degree of the connectivity matrix unchanged.
On the nodal level, two graph measures were calculated: the betweenness centrality and the degree. The betweenness centrality represents the fraction of shortest paths that pass through a certain node. Degree is a measure of interaction and can be calculated as the summation of all functional connections per node.
In general, graph measures are known to depend on the number of nodes and the average degree within a network (Wijk, Van, Stam, & Daffertshofer, 2010). Therefore, to obtain robust measures, every graph measure was calculated for a range of densities. The lowest density was set to 28% to prevent disconnected networks. The full density range comprises densities between 28 and 50% (in steps of 2%). Above 50%, connections are thought not be physiologically meaningful (Hosseini, Hoeft, & Kesler, 2012;Kaiser & Hilgetag, 2006).
The area under the curve was calculated over this whole density range to obtain one robust, representative value for the graph measure per patient, per time-point, and in case of the nodal analysis also per node.

| Statistical analysis
In this study, functional connectivity, represented by various graph measures, was compared between T1 and T2 ( Figure 1). Here, ΔGM is the change in graph measure (GM T2 − GM T1 ), and referred to as the effect size. Because of non-normality of the graph parameters (see Appendix B), nonparametric permutation tests using 1,000 permutations were performed to investigate the difference between sham and verum stimulation on graph measures (ΔGM sham vs ΔGM verum ).
Significance level was set to p < .05 for the whole-brain analysis.
On the nodal level, additional multiple comparison correction was applied via the Holm-Bonferroni method, using the number of nodes for correction, but all findings with p < .05 were reported. Post-hoc t tests were used to investigate the direction of the effects.

| Spatial distribution
To study the assumption that the effect of aiTBS distributes via func-

| Biomarker investigation
To investigate the predictive value of graph parameters on the clinical response to aiTBS, the baseline graph measures were correlated with the change in HDRS in terms of percentage (T2 with respect to T1 in the subgroup of patients receiving verum stimulation). Here, this means the lower the scores on HDRS changes in terms of percentage, the better the clinical response. Only significant correlations (p < .05) were reported.

| RESULTS
Given five drop-out patients (due to a different diagnosis retrospectively, clinical improvement before the stimulation, or incomplete or wrongly timed MRI datasets), exclusion of seven patients (due to excessive motion in the MRI dataset at either T1 or T2), exclusion of three subjects based on the tSNR criterion, and three subjects did not have connected graphs within the density range, data from 32 patients were used for analysis. Of these patients, 14 received sham stimulation between T1 and T2 (arm A in Figure 1) and 18 received verum stimulation (arm B in Figure 1). Patient details and results on the clinical outcome of this stimulation protocol can be found in Duprat et al. (2016). Based on the tSNR criteria, 19 nodes (represented in red in Figure 2) were removed. Detailed information about the excluded nodes can be found in Appendix A, Table A2.

| Whole-brain network topology changes
On the whole-brain level, stimulation caused a significant effect on clustering coefficient and global efficiency (p values <.01, <.01, .072, and .607 for clustering coefficient, global efficiency, modularity, and small-worldness, respectively) (Appendix B). However, the effects did not differ between the subgroups receiving sham and verum stimulation. An overview can be found in Table 1. As can be seen in Table 2, changes in graph measures were not significantly correlated with changes in clinical outcome.    Table 4 shows an overview of the significant (p < .05) findings. A full overview can be found in Appendix C.

| Propagation of effect via functional connections
The   Table 5.
No significant correlations were found. However, only based on findings within nine nodes, a large negative slope was found between the functional connectivity and the effect size of stimulation on the degree. This suggests that the effect of aiTBS on degree depends on the functional connectivity with the stimulation site: higher functional connectivities are linked to higher effect sizes (lower p values). Figure 6 shows an overview of the baseline whole-brain graph measures versus the percentage of change in HDRS score, after versus before verum stimulation. Table 6 shows the statistical values.  effects and the belonging statistics can be found in Table 7.

| DISCUSSION
This study aimed to use graph theoretical analysis to investigate the effects of the relatively new accelerated stimulation protocol to treat MDD patients, namely aiTBS, on the brain's network organization.

| The effect of aiTBS on graph measures
4.1.1 | Whole-brain results On the whole-brain level, no significant differences between the effects of verum stimulation versus sham stimulation were found, and changes in graph measures did not correlate with changes in depression severity symptoms. Previous studies (Ajilore et al., 2014;Lim et al., 2013) found no differences between graph measures clustering

| Nodal results
On the nodal level, some nodes showed significantly different responses to verum and sham stimulation. Because these nodes are   after an accelerated high frequency rTMS paradigm was associated with significant increases of GABA (γ-aminobutyric acid) concentrations in the stimulated area (the same left DLPFC spot that was also targeted here in this study) (Baeken, Lefaucheur, & Van Schuerbeek, 2017b). These GABA increases must be primarily considered as an "excitation" of GABAergic inhibitory neurons and pathways (Lefaucheur, Drouot, Ménard-Lefaucheur, Keravel, & Nguyen, 2006).
Both Kang et al. (2016) and Liston et al. (2014) have reported reductions in connectivity after 10 Hz rTMS, which is also assumed to have excitatory effects. However, one needs to keep in mind that according to Huang et al. (2005), the standard iTBS protocol is thought to result in excitatory effects. The aiTBS protocol is a modified form of the original iTBS protocol, not only in the number of pulses but also in the number of sessions. As it is known that modifications of stimulation protocols are able to reverse the polarity of the after-effects (Gamboa et al., 2011;Gamboa, Antal, Moliadze, & Paulus, 2010;Gentner, Wankerl, Reinsberger, Zeller, & Classen, 2008;Murakami, Müller-Dahlhaus, Lu, & Ziemann, 2012), it remains to be determined whether the net effects in the stimulated and connected areas are excitatory or inhibitory.

| Specific nodal effect
The most significant result (also the only finding that survived Bonferroni correction) was observed in the right supplementary motor area.
Whereas the betweenness centrality increased after verum stimulation, it decreased after sham stimulation. This means that shortest paths between brain regions pass the right supplementary motor area.  Potential of whole-brain graph measures clustering coefficient, global efficiency, modularity, and small-worldness to predict the percentage of clinical change of verum aiTBS. Statistical details can be found in Table 6 [Color figure can be viewed at wileyonlinelibrary.com]  As TMS has been linked to changes in psychomotor performance before in the healthy as well as depressed state, this is of interest to explain to some extent the working mechanisms of this kind of stimulation. For instance, Baeken et al. (2010a) found improved psychomotor performance after high-frequency rTMS treatment in medication resistant depressed patients. Also Hoeppner et al. (2010) showed a trend toward reduction of psychomotor agitation in MDD after high frequency rTMS. Our current findings indicate that left DLPFC aiTBS indeed may affect cortical areas involved in (psycho)motor actions.
In addition, more exploratory analyses revealed that the aiTBS treatment protocol shows effects on several (sub)cortical areas that can be linked to the pathophysiology of depression. For example, the effects of sham and verum aiTBS on degree differ in the left cinguloopercular nodes, which are part of the cingulo-opercular network comprising the bilateral dorsal anterior cingulate cortices (dACC), the anterior insula, anterior prefrontal cortex, and the anterior thalamus (Sylvester et al., 2012). This network integrates visceral, autonomic, and sensory data to assess the homeostatic relevance or "salience" of internal and external stimuli, and the maintenance of tonic alertness or sustained attention (Sadaghiani & D'Esposito, 2015). The network also clears noisy information, suppresses distraction, and keeps cognitive faculties available for current processing demands (Sadaghiani & D'Esposito, 2015). Abnormalities in this network have been reported for obsessive compulsive disorder (OCD) (de Vries et al., 2017), psychosis (Sheffield et al., 2017), and mood and anxiety disorders (de Witte & Mueller, 2016). Of interest, Wu et al. (Wu et al., 2016) showed that depression symptom severity was significantly correlated with the connectivity values of this network. Indeed, increased activity in the dACC or insula during response conflict has been reported during negative mood states (Disner, Beevers, Haigh, & Beck, 2011).
Furthermore, several nodes that showed significantly different effects of verum and sham aiTBS, such as for example the parahippocampal nodes, nodes within the prefrontal cortex, and the posterior cingulum node, belong to the default mode network (DMN). The DMN is found to be activated during resting-state functional imaging and de-activated when performing cognitive tasks (Fox et al., 2005;Smith et al., 2009). When the brain is not engaged in externally driven cognitive processing, self-referential processes are believed to predominate (Gusnard, Akbudak, Shulman, & Raichle, 2001). When clinically depressed, more activity in the DMN is observed (Disner et al., 2011). Changes in DMN activation have earlier been linked to antidepressant responses.

| Spatial distribution of aiTBS effects
Previous studies have already shown distributed "network-effects" of TMS (Fox et al., 2013;Fox et al., 2014). In this study, using nodes showing significantly different effects between verum and sham stimulation, the correlation between effect sizes and functional connectivity strengths did not reach significance. This indicates that the propagation of aiTBS-effect from the stimulated area is not directly linked to the strength of the functional connections. Considering the network-hypothesis, we hypothesize that the indirect effects of TMS occur at different levels. After the activation of brain areas connected to the stimulation site are activated, in the following steps, the brain areas connected to those areas are activated and so on. This could, at least partly, explain the occurrence of increases and decreases in graph measures in distinct areas of the brain.

| Graph measures as biomarker
Clinical improvement was associated with higher baseline clustering coefficient or global efficiency on the whole-brain level. This indicates that all nodes within the whole brain are better integrated. The effect of verum stimulation therefore seems to propagate more easily through the whole-brain via functional connections, also to deeper structures involved in the deregulated neurocircuitry of depression.
On the nodal level, we found that graph measures in multiple nodes showed potential to predict the clinical effect. For example, a positive correlation between the baseline betweenness centrality and clinical effect was found in the left caudate nucleus. So lower betweenness centrality might be advantageous for clinical outcome.
Given that the caudate has neural innervation from amongst others the prefrontal cortex, our left caudate nucleus findings could be linked to the application of left-sided stimulation (Kang et al., 2016). Indeed, stronger connectivity between the dorsal prefrontal cortex and the (dorsal) caudate has been associated with depression severity (Furman, Paul Hamilton, & Gotlib, 2011;Kerestes et al., 2015). Furthermore, observations of increased connectivity with the DLPFC and the more ventral parts of the ACC in MDD was associated with heightened cognitive regulation of affect, usually problematic when clinically depressed; whereas reduced connectivity with the caudate results in worsening symptoms such as anhedonia, reduced motivation, and psychomotor dysfunction (Davey, Harrison, Yücel, & Allen, 2012). Of note, although the sgACC was not implicated in our findings, the structural and functional connections between the striatum (caudate) and the (sg)ACC are well known (Gabbay et al., 2013). In treatment-resistant depression, the sgACC has been proposed as biomarker for response for a variety of interventions, including rTMS treatment (Fox et al., 2012;Fox et al., 2013;Weigand et al., 2017).
However, for the latter application, the functional connectivity findings are not that straightforward (Baeken et al. 2017a;Baeken et al., 2014) and the aiTBS treatment delivered to the left DLPFC may have different neurobiological effects on the reward system (including the caudate), based on the level of anhedonia in the depressive state (Duprat, Wu, De Raedt, & Baeken, 2017). Indeed, it remains to be determined whether the left DLPFC is the best target to stimulate.
Other prefrontal areas, such as the dorsomedial prefrontal cortex have been successfully stimulated in depressed patients (Downar et al., 2014), and alternatively when facing nonresponse, the right orbitofrontal cortex (OFC) was found to be an excellent alternative . The right OFC is considered as a 'non-reward' nexus (Cheng, Rolls, Qiu, Liu, & Tang, 2016) showing reduced functional connectivity in MDD patients. Together with our own findings on clinical improvement combined with baseline striatal (caudate) betweenness centrality, these observations suggest that left DLPFC aiTBS could be successful for a selected cohort of patients.
Furthermore, the degree in the right amygdala was significantly correlated with the clinical effects of verum aiTBS, suggesting that less connections to the right amygdala could be predictive for better clinical responses. Given that the amygdalae are involved in (in)effective emotion regulation in stress-related disorders (Gold & Chrousos, 2002;Perlman et al., 2012) and in particular the right amygdala is implicated when processing negative information stressful events (Baeken et al., 2010b;Mothersill & Donohoe, 2016), it is of interest to note that increased baseline and sustained amygdala activity to antidepressant treatment is associated with clinical nonresponse in major depression (Fonseka, Macqueen, & Kennedy, 2018).

| General limitations
This study has some general limitations that need to be considered.
Notwithstanding that rs-fMRI is a unique and powerful tool to investigate human brain organization, it is based on an inherently ambiguous measure reflecting dynamic couplings that are not yet fully understood. Interscan rs-fMRI data have shown great variability. For example, Ning et al. (2017) aimed to derive the optimal TMS stimulation position based on functional connectivity between the DLPFC and the sgACC and showed different results using resting-state data from same subjects at different time-points. Longer rs-fMRI scans were suggested to reduce this variability. Moreover, various patient-specific factors may also influence the outcome of a stimulation protocol (Silvanto & Pascual-Leone, 2008

MOTION PARAMETERS
Overview of the number of volumes per patient after removing the ones with framewise displacement >.3. Datasets with <200 volumes were discarded from further analysis. Note from Figure A1 that 5 datasets in T1 show excessive motion versus 6 in T2. With four overlapping datasets, this led to removal of 7 datasets for further analyses.

NODE SELECTION
The first 264 nodes are resulting from the Power parcellation, as described in Power et al. (2011)). Thirteen additional nodes were appended, in accordance to Drysdale et al. (2016). An overview can be found in Table A1.   Figure A2 displays an overview of the mean tSNR (averaged over patients) per node. The nodes that were discarded for further analysis are listed in Table A2. Figures B1 and B2 show the distributions of whole-brain graph measures for the subgroup of patients receiving sham and verum aiTBS, respectively.

FIGURE A2
Overview of mean tSNR and standard deviation averaged over all subjects, per node. Nineteen nodes showed tSNR < 40 in more than 10% of the datasets and were excluded for further analysis (Table A2) [Color figure can be viewed at wileyonlinelibrary.com]   Table C1.