Proof of concept study to develop a novel connectivity-based electric-field modelling approach for individualized targeting of transcranial magnetic stimulation treatment

Resting state functional connectivity (rsFC) offers promise for individualizing stimulation targets for transcranial magnetic stimulation (TMS) treatments. However, current targeting approaches do not account for non-focal TMS effects or large-scale connectivity patterns. To overcome these limitations, we propose a novel targeting optimization approach that combines whole-brain rsFC and electric-field (e-field) modelling to identify single-subject, symptom-specific TMS targets. In this proof of concept study, we recruited 91 anxious misery (AM) patients and 25 controls. We measured depression symptoms (MADRS/HAMD) and recorded rsFC. We used a PCA regression to predict symptoms from rsFC and estimate the parameter vector, for input into our e-field augmented model. We modeled 17 left dlPFC and 7 M1 sites using 24 equally spaced coil orientations. We computed single-subject predicted ΔMADRS/HAMD scores for each site/orientation using the e-field augmented model, which comprises a linear combination of the following elementwise products (1) the estimated connectivity/symptom coefficients, (2) a vectorized e-field model for site/orientation, (3) rsFC matrix, scaled by a proportionality constant. In AM patients, our connectivity-based model predicted a significant decrease depression for sites near BA9, but not M1 for coil orientations perpendicular to the cortical gyrus. In control subjects, no site/orientation combination showed a significant predicted change. These results corroborate previous work suggesting the efficacy of left dlPFC stimulation for depression treatment, and predict better outcomes with individualized targeting. They also suggest that our novel connectivity-based e-field modelling approach may effectively identify potential TMS treatment responders and individualize TMS targeting to maximize the therapeutic impact.


INTRODUCTION
Major depressive disorder (MDD) is characterized by depressed mood or loss of interest or pleasure [1]. Symptoms include feelings of worthlessness; indecisiveness or difficulty concentrating, change in weight, appetite, and sleep; psychomotor agitation; and suicidal ideation [1]. According to the World Health Organization (WHO) World Mental Health Surveys Initiative [2], lifetime prevalence rates within the United States average 16.9%, with women showing higher prevalence than men (20.2% vs. 13.2%) [2]. The most common treatments for MDD include antidepressant medication, psychotherapy/cognitive behavioral therapy, or a combination of the two (https://www.nimh.nih.gov/ health/statistics/major-depression.shtml). Although antidepressant use is prevalent, there are reports of mixed efficacy [3]. Repetitive transcranial magnetic stimulation (rTMS) to the left dorsolateral prefrontal cortex (dlPFC) has been cleared by the US Food and Drug Administration, and has been shown to be effective as an antidepressant treatment in treatment resistant depression, a patient population not responsive to or tolerant of previous antidepressant medication trials [4][5][6]. Although many studies find a significant reduction in depressive symptoms following rTMS treatment [7][8][9][10][11][12][13], there have also been studies finding limited efficacy for major depression [14,15], leaving room for improvement.
One possible explanation for the heterogeneity in TMS treatment response for major depression is that the standard scalp-based targeting approach (i.e., the 5 cm rule) does not account for individual variations in cortical anatomy, and often selects regions outside the dlPFC, which may explain some of the modest therapeutic effects [16]. Image-guided TMS incorporating MRI-neuronavigation, taking into account individual differences in brain anatomy, can target specific functional brain networks with greater accuracy [17,18]. Using structural MRI to locate a specific site at the junction of BA 9 and 46 has also shown greater clinical efficacy in reducing depression than the traditional 5 cm rule [15]. More recent approaches have used functional MRI (fMRI) to target regions that are downstream of the dlPFC, like the subgenual anterior cingulate cortex (sgACC) [19][20][21] that are thought to be dysfunctional in MDD patients [22][23][24][25]. Key to this approach is to use functional connectivity of the downstream site to identify the dlPFC site with the strongest connectivity. Finally, using e-field modelling can improve accuracy even further. Indeed, a recent report estimates the difference between conventional neuronavigated targeting and e-field informed target to be up to 14 mm [26].
Connectivity-based approaches represent the current cuttingedge in the field and are an improvement on the standard scalpbased targeting approach used clinically. However as currently implemented, they typically suffer from two fundamental limitations, which we address in the current work using a novel computational model based on whole brain functional connectivity and electric-field modelling. The first limitation is that by focusing on a single downstream region to define the TMS target they do not account for off-target effects in other downstream regions. To address this limitation, we use electric-field modelling to estimate the changes in connection strength across the entire brain, thus accounting for changes in connectivity that would otherwise be considered "off-target". The second limitation is that most connectivity-based targeting approaches do not typically quantify the link between functional connectivity and symptoms, making it difficult to make a quantitative prediction about symptom change following neuromodulation. To address the second limitation, we use a PCA regression approach to quantify the relationship between functional connectivity and depression symptoms, and use the resulting coefficients to estimate the net effect of TMS on connectivity and thus symptoms. The overall concept is to generate a predicted symptom change score following neuromodulatory TMS at a particular site and stimulation. We then iterate the model across sites and orientations to find the site/orientation combination that leads to the maximal reduction in MDD symptoms. Importantly, because this approach is data-driven it has an added benefit of being both site and symptom independent.
The current work is a proof of concept demonstration of this novel targeting approach using structural and fMRI collected from the Dimensional Connectomics of Anxious Misery dataset (See Seok et al., 2020 [27] for full characterization of cohort). Although we used a dimensional approach to recruit a diagnostically heterogeneous patient sample, it is important to note the corresponding DSM-5 diagnoses represented in the sample all have data supporting therapeutic use of TMS [4][5][6][7][8][9][10][11][12][13][28][29][30][31][32][33][34][35]. The goal is to characterize the relationship between symptoms and connectivity, then model the effect of TMS on connectivity in order to predict the change in symptoms that would be likely to occur given stimulation at a particular site and orientation. In the current work, we apply this modeling approach to two targets, a positive control left dlPFC and a negative control motor cortex (M1), based on targeting guidelines from the literature (See Supplementary Table 1 and 2). We use bootstrapping to estimate the reliability of the targets identified using this approach at the single-subject level, and permutation testing to determine whether this approach makes group-level predictions that are consistent with the literature [36][37][38].

MATERIALS AND METHODS Participants
One hundred eighteen participant (Age (years): M = 28.42, SD = 7.92) entries were pulled from the Dimensional Connectomics of Anxious Misery dataset [27], and included anxious misery (AM) participants (93 total, 64 were females) and healthy controls (25 total, 16 were females; See Supplementary Table 3). Mean (SD) age of the AM cohort was M = 28.29 (8.2) and of healthy controls M = 28.92 (6.91). Group assignment was based on the Neuroticism score on the NEO Five-Factor Inventory (NEO FFI) rather than on diagnosis in keeping with the RDoC concept of finding symptom domains across diagnoses. The population distribution of Neuroticism scores was estimated from a previously collected sample of 635 adults [39]. Different distributions were used for males (M = 19.1 ± 7.1) and females (M = 22.2 ± 7.9). Participants were assigned to the AM group if their Neuroticism scores were at least 1 SD above the mean (males ≥ 26.2 and females ≥ 30.1). The mean Neuroticism score for males and females in the AM group were as follows (Males = 33.04 ± 3.39; Females = 32.82 ± 4.65), while the mean Neuroticism score for males and females in the HC were as follows (Males = 23.33 ± 37.07; Females = 23 ± 3.58). Among the AM group, 24 were medicated, while 69 were unmedicated. In post-hoc analyses, the cross-diagnostic AM cohort included a range of different primary diagnoses, namely MDD (N = 41), persistent depressive disorder (N = 9), social anxiety disorder (N = 5), generalized anxiety (N = 21), and PTSD (N = 17).
Basic Inclusion/Exclusion criteria included MRI contraindications, histories of certain neurological or cognitive disorders and events, histories of exclusionary psychiatric conditions or any other factors that may have affected patient safety (See Seok et al [27]. for complete list). In addition, participants in the AM cohort needed to score a minimum of one (1) standard deviation above a point estimate of the general adult population mean (distributions for males and females were different) in the NEO FFI [40] in order to be included in the study [27]. Conversely, healthy controls needed to score one (1) standard deviation or more below the mean point estimates to be eligible to participate. Subjects were recruited by means of IRB-approved advertisement and phone calls. All participants signed an informed consent form, and the protocol was approved by the Institutional Review Board for human subject research at the University of Pennsylvania. The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.

Procedure
On the initial visit, written consent was obtained from all eligible participants prior to enrolling in the study, followed by a baseline visit with study staff to ascertain their medical history and demographic information. In addition, a research version of the Structured Clinical Interview for DSM-5 (SCID-5-RV) [1] was conducted by a trained member of the staff to document their psychiatric diagnosis and history, and the remaining assessments (See Below) were conducted. In a follow-up visit, a 2 h fMRI scan was conducted, where structural, diffusion, resting state, and task-based fMRI data were collected as part of the larger project.

Assessments
Many clinical measures were used to identify the behavioral and cognitive features of AM. In this project, eligible participants completed 17 selfreport and 3 clinician-administered measures. The clinician-administered measures included the Montgomery-Åsberg Depression Rating Scale (MADRS) [41] and the Hamilton Depression Rating Scale (HAMD) [42] in order to assess depression severity.

Model steps
Image acquisition and preprocessing. T1, T2, and resting state fMRI were acquired on a Siemens Prisma 3T using a 64-channel head coil. SimNIBS (Version 2.1) was used to generate 3D head and coil geometries and compute the e-field models [43]. fMRI data were preprocessed using the Human Connectome Project Pipelines v4.0.1 and the eXtensible Connectivity Pipeline (XCP Engine) [44]. For a full description of the image acquisition and preprocessing steps, see the Supplementary Methods. PCA regression. For the first step in our analysis (See Fig. 1 and Appendix for detailed equations), we used principal components analysis (PCA) regression to model the relationship between functional connectivity and symptoms [45,46,Eq. A2]. We began by conducting a PCA on the functional connectivity data for each group (i.e., AM and control) to reduce noise and minimize collinearity [47]. We restricted attention to the subset of components that explained largest proportions of variability, with the assumption that these contain most of the relevant signal. For this we used a geometric approach that identified the eigenvalue at the elbow of the scree plot (i.e., the eigenvalue furthest from the hypotenuse connecting the first and last eigenvalues). The principal component scores for each of these components were then used to model symptom scores (MADRS and HAMD) in a multiple linear regression. The resulting coefficient estimates were then projected back into feature space (Eq. A3). The result was a vector where each element represents the coefficient for the corresponding functional connection in a linear model predicting symptom score. N.L. Balderston et al.

Selection of sites and orientations.
To characterize the dlPFC, we defined a series of 17 sites along the anterior-to-posterior gradient of the L dlPFC, which was prioritized based on prior clinical efficacy findings (See Supplementary Table 1). This site vector was anchored by three therapeutic sites described in [38] [−44, 40, 29] BA46), and extends 12 mm posterior to the 5 mm site and 12 mm anterior to the BA46 site. We also plotted two additional sites along the vector. First is the site corresponding to the Beam/F3 stimulation site typically identified using the EEG 10/20 system [48]. Second is the site along the vector with the maximal anti-correlation with the sgACC. Note that this second site was calculated from the rsFC data for each subject and plotted on the individual subject heatmaps.
To characterize the hand knob, which was used as a negative control region, we defined a series of seven sites along the medial to lateral axis of the primary motor cortex in the left hemisphere (See Supplementary  Table 2). This anchor was centered on the hand knob, which was identified from Neurosynth using an association test with the search term "hand". The center of mass for the primary motor cortex cluster was used as the Fig. 1 The PCA-Regression data reduction approach used to summarize the relationship between symptoms and connectivity in the current model. A Resting state functional connectivity (rsFC) is calculated using Pearson's correlation across all subjects for all regions in the Gordon atlas [49]. B A principal component analysis (PCA) is used to identify orthogonal components in the rsFC data, and a geometric approach is used to identify a minimal number of components that explain a maximal proportion of the variability. C Component scores for the selected components are extracted and entered into a multiple linear regression to predict symptoms (D). The PCA loadings from the selected components (E) are combined with the coefficient vector from the regression (F) using matrix multiplication to create the output vector (G), which is used to represent multiple regression coefficients projected into the rsFC feature space. E-field connectivity matrices. For each participant, the e-field map was downsampled (i.e., averaged across voxels in each parcel) to the Gordon atlas [49] and concatenated into a 1 × 333 vector. This vector was then sorted by e-field magnitude to resemble a scree plot, and thresholded at the value furthest from the hypotenuse connecting the first and last value in the vector (See Fig. 2). To estimate how stimulation at a specific site and orientation might affect connectivity, the downsampled e-field vector was used to form a 333 × 333 matrix where the values in the matrix represent the average current induced in the ROIs for each connection.

E-field/connectivity symptom estimates
The predicted neuromodulatory effect of TMS on symptoms was modelled using each participant's e-field model matrix [50], the single-subject connectivity matrix, and the estimated group-level coefficients derived from the symptom/connectivity regression, according to Eqs A2-3. These calculations were iterated across all sites and orientations and spatially smoothed using a Gaussian kernel (SD = size(heatmap)/10; See MATLAB function smoothdata) to create the heatmap shown in Fig. 3, where each value in the heatmap represents the predicted relative change in symptoms (sans application of an empirically-determined constant scaling factor C) for that subject given hypothetical TMS administered at that particular site and orientation.

Statistical analysis
Individual subject predictions. One goal of this project was to determine whether the modelling approach above could be used to make targeting decisions for the application of TMS to single subjects. Accordingly, to be effective our model would need to show a systematic effect of site and orientation on predicted symptom change in order to be used for singlesubject targeting. To examine this possibility, we created single-subject heatmaps corresponding to the predicted relative change in symptoms as a function of site and orientation, and identified the point in these heatmaps that yielded the maximum predicted relative decrease in symptom scores, corresponding to the optimal site/orientation of stimulation for that subject.
We then used bootstrapping to estimate the reliability of this targeting approach. For each of the 10,000 iterations of the bootstrapping analysis, we sampled subjects with replacement from the original distribution, recalculated the PCA regression and the corresponding vector (See Appendix: EQ5). Next we recreated the individual subject heatmaps, and identified the optimal stimulation site. Across bootstrap iterations we calculated the Euclidean distance between the single-subject original optimal site/orientations and the bootstrap-calculated optimal site/ orientations. Then for each subject, in order to form a confidence region, we identified the value corresponding to the 95th percentile from the distance metric bootstrap distribution, which was used to display the reliability of the single-subject target in the heatmaps.
Group-level predictions. To determine whether our novel modelling approach predicted significant changes at the group level, for each site/ orientation we conducted one-sample t-test of the predicted change in symptoms (i.e., predicted post treatment symptom score-baseline symptom score) against zero. A significant negative t-score would indicate a predicted decrease in symptoms at the group level. Accordingly, the null hypothesis for each site/orientation combination is that the model predicts no significant group-level change in symptoms following a hypothetical treatment at that site/orientation. We repeated this process for all site/ orientation combinations, and generated heatmaps from the resulting tscores. We used cluster-based permutation testing to correct for multiple comparisons across site and orientation [51]. First we conducted 10,000 permutations where the sign of the predicted symptom change was randomly flipped across subjects and permutations. For each permutation, we calculated t-scores at each site/orientation combination within the heatmap. We thresholded these t-scores using a two-tailed alpha of 0.025. Next we identified clusters of t-scores surviving the threshold, summed the values across elements within the cluster, identified the cluster with the maximum summed t-score, and built a null distribution with these max summed t-score values. We then identified the critical values corresponding to the 1.25th and 98.75th percentiles (with alpha = 0.025 to correct for the number of regions of interest tested) of the max summed t-score distribution generated from the permutation procedure. Finally, a two-step thresholding process was implemented to threshold the group-level statistical maps. First, the individual t-scores were thresholded using an alpha of 0.025. Then, summed t-scores from the surviving clusters from these heatmaps were thresholded using the critical values from the max summed t-score null distribution.

PCA regression
Separate PCA regression analyses were conducted for the AM group and the control subjects. For the AM group, 11 signal components were selected, and these components explained 31% of the variability in the functional connectivity signal, 16.8% of the   Table 1 for PCR component β-values for MADRS and HAMD scores and Supplementary Fig. 1 for component item  loadings). For the control group, six components were selected, and these components explained 40% of the variability in the functional connectivity signal, 33.7% of the variability in the MADRS scores, and 44.3% of the variability in the HAMD scores (See Supplementary Fig. 2).
The stability of the PCA was evaluated using results from the bootstrapping analysis described below. For each of the 10,000 iterations of the bootstrapping analysis, we sampled subjects with replacement from the original distribution, recalculated the PCA and calculated both the number of selected components as well as the total variability explained by the selected components. For the AM group, there was an average of 15.02 ± 2.89 components selected across the bootstrapping iterations and these components explained an average of 53.71 ± 6.63% of the variability in FC across the sample. For the HC group, there was an average of 8.99 ± 4.11 components selected across the bootstrapping iterations and these components explained an average of 75.14 ± 18.22% of the variability in FC across the sample.

Individual subject predictions
Single subject heatmaps from the AM group show clear systematic patterns of predicted symptom change as a function of both site and orientation (See Fig. 4), suggesting that this modelling approach can be used to simultaneously optimize these two parameters at the single subject level. These patterns are characterized by punctate "cool spots" in the heatmap that identify a clear site/orientation combination that yields maximal predicted symptom reduction. To characterize the variability in our targeting approach, we conducted 10,000 bootstrapping iterations and identified the optimal site/orientation combination on each. We then calculated the Euclidean distance between the bootstrap-derived values and the original values and plotted a circle with an area that corresponds to the 95th percentile in the bootstrapped Euclidean distance distribution on the single subject graphs shown in Fig. 4 and Supplementary Fig. 3. Importantly, this pattern seems to be specific to those in the AM group as control subjects do not show the characteristic "cool spots" in their heatmaps (See Supplementary Fig. 3).

Group-level predictions
To determine whether the model made significant predictions at the group level, we calculated the predicted change in depression symptoms for each site/orientation combination by subtracting the observed pre-TMS MADRS/HAMD from the predicted post-TMS MADRS/HAMD scores. These difference scores represent the predicted change in depression symptoms, given a specific site/ orientation of stimulation. We then compared these values to 0 (i.e., no effect of TMS) using a single sample t-test. We corrected for multiple comparisons across sites/orientations the permutation testing procedure described above. As shown in Fig. 5 our model predicts a significant decrease in MADRS and HAMD scores in AM subjects with both left and right dlPFC stimulation. In the left hemisphere, the cool spots are adjacent to BA9 at 30-45°o rientation. In the right hemisphere, the cool spots are positioned slightly anterior to the cool spots in the left hemisphere, and strongest at 120 and 135°orientation. Consistent with the results at the single-subject level, these results suggest that there is a systematic relationship between both site and orientation and predicted symptom change. Importantly, these results are limited to the AM subjects, as the model predicts no significant changes in depression symptoms in the healthy control subjects (See Supplementary Fig. 4). In contrast to the data from the dlPFC, the model does not predict any significant changes in depression symptoms for site/orientation combinations centered on the hand region of M1 for either measures (MADRS/HAMD) in AM patients ( Supplementary Fig. 5) or healthy controls ( Supplementary Fig. 6).
It is important to assess the degree to which variability in image quality affects the predictions made by our model. Accordingly, we calculated tSNR and average motion (RMS across TRs) and correlated these values with the predicted symptom change scores for MADRS and HAMD for both the AM and HC groups (See Supplementary Table 4). For the AM group, image quality metrics were uncorrelated with symptom predictions generated by our model, suggesting that differences in image quality cannot account for our effects in this group. In contrast, for the HC group, average motion across TRs was marginally correlated with predicted changes in both MADRS and HAMD, suggesting that image quality may have played a part in the predictions for this group. Given that the predicted symptom scores reflect the variability in both the symptom and the imaging data, this group difference could be due to the substantially different levels of variability in MADRS/HAMD scores in the AM group compared to the HC group.

DISCUSSION
Here we present a proof of concept methodological study where we propose a model that uses functional connectivity, symptom scores, and electric-field modelling to predict symptom change following a hypothetical course of neuromodulatory treatment with TMS. Our model calculates the relationship between symptoms and connectivity, as well as the effect of TMS on connectivity to yield a predicted change score in the symptom, assuming TMS is delivered at the site and orientation specified in the model. We then tested the model systematically across sites and coil orientations in two cortical locations serving as positive (dlPFC) and negative (M1) controls. For the dlPFC, our model suggests that stimulation at this site would lead to a significant reduction in depression symptoms (MADRS/HAMD) in AM patients but not controls, and the magnitude of this potential effect should vary systematically as a function of both site and orientation. Importantly, this observation is specific to the dlPFC, as our model  suggests that stimulation at M1 would likely not result in a significant decrease in MADRS or HAMD, regardless of position or orientation. In addition, these observations seem to be reliable at both the single-subject and the group level, suggesting that this approach has sufficient ability to identify subject-specific TMS targets to maximize symptom reduction. Consistent with previous work implicating the junction of BA9 and BA46 as a potentially optimal site of stimulation for depression symptoms [52][53][54][55][56], our model suggests that stimulation marginally posterior to this site would maximally reduce both MADRS and HAMD scores across subjects compared to the other sites tested [36]. Our model adds to the work suggesting that rsFC networks can inform TMS targeting [37], and our results are consistent with the findings that individualized rsFC maps may be most informative [38]. Previous research has shown that regions of the dlPFC that are more strongly anticorrelated with sgACC tend to show better clinical efficacy when targeted with therapeutic TMS [52,57], and the closer the stimulation site is to this optimally anticorrelated dlPFC location, the better the clinical outcome [55,56,58]. Importantly, active rTMS to the l dlPFC has been shown to reduce anticorrelation between the dlPFC and sgACC [54,59,60], and prospective targeting based on individualized dlPFC-sgACC connectivity leads to a greater reduction in symptoms compared to traditional targeting [54,55]. In addition, there is some research to suggest that individualized targeting of BA46 through a parcel guided approach can improve treatment response in patients who have previously failed standard TMS targeted at the dlPFC using the 5 cm rule [36]. Based on this work, it could be argued that the current state-of-the-art for individualized TMS targeting for depression is to target sites in the left dlPFC with maximal connectivity with the subgenual cingulate cortex.
However, a review of 33 studies with pre-post resting state shows that although TMS induces robust changes in rsFC, the majority of effects are outside of the stimulated functional network [61]. Indeed, the changes in functional connectivity following rTMS treatment for depression can be seen throughout the brain in regions that are important for affective responding (e.g., insula, amygdala, inferior parietal lobule, etc.), suggesting that the net effects on symptoms may be driven by these largescale network connectivity changes [62]. Consistent with this idea, treatment response has also been linked to rsFC changes outside of the sgACC [63,64]. Together these results suggest that a whole brain connectivity model may outperform single connection targeting approaches, like the ones described above.
In keeping with these results, our whole brain connectivity approach has several potential benefits over single connection targeting approaches [52][53][54][55][56]. First, our model uses e-field modelling to account for the variability in electric field spread due to individual differences in head and cortical anatomy, and uses this information to inform targeting. Second, our approach accounts for other connections in the brain that (1) may be Fig. 5 Group-level heatmaps plotting dlPFC predictions for the anxious misery group. A Heatmap representing the predicted MADRS scores following a hypothetical course of TMS treatment to the left dlPFC. B Heatmap representing the predicted MADRS scores following a hypothetical course of TMS treatment to the right dlPFC. C Heatmap representing the predicted HAMD scores following a hypothetical course of TMS treatment to the left dlPFC. D Heatmap representing the predicted HAMD scores following a hypothetical course of TMS treatment to the right dlPFC. Y-axis represents coil orientation. Red circles represent sites where the change in MADRS/HAMD scores was not statistically different from 0. 5 cm = Site in vector corresponding to the 5 cm rule target commonly used in therapeutic applications for depression [38]. BA9 = Site in vector corresponding to Broadmann area 9 [38]. F3 = Site in vector closest to the BEAM/F3 target commonly used in therapeutic applications for depression [48]. BA46 = Site in vector corresponding to Broadmann area 46 [38].
impacted by TMS treatment, and (2) contribute to symptom reduction. Importantly, our model uses a data driven approach to weight these connections by the amount of variability they explain in the depression symptoms. Finally, our approach uses agnostic data-driven estimates between behavioral characteristics and functional connectivity to drive targeting. Therefore, this model can be extended to other symptoms/behavioral characteristics where reliable measures of behavioral markers are available (e.g., addiction, anxiety, etc.).
The main focus of the paper was depressive symptoms, which we captured using scores on the MADRS and HAMD. However, it is important to note that the model described here can be applied to other symptom measures as well. In previous work, we have used various symptom clustering approaches to extract distinct symptom measures from available clinical data [65], resulting in several distinct symptom clusters. Although out of the scope of the current work, future studies should determine the degree to which these other distinct symptom clusters could potentially be used to direct TMS targeting. Indeed, defining biotypes for depression and anxiety starting with neural circuit dysfunction [66], may yield additional symptom-specific targets. Accordingly, the model presented here could be used not only to individualize targeting at the single subject level, but it could potentially be extended to individualize targeting at the single session level, where each TMS session in a course of treatment is targeted at a specific symptom or set of symptoms, thereby maximizing the therapeutic effect.
Our approach differs from other connectivity-based targeting approaches, including that of similar paper by Siddiqi et al (2020) [67], in two fundamental ways. First, these studies typically calculate connectivity with the stimulation site using a voxelwise approach. This approach assumes a direct correspondence between the TMS effect and the functional connectivity. That is, it doesn't take into consideration the focality of the TMS-induced electric field. Our approach weights the connectivity changes proportional to the efield model, which should yield a more accurate estimate of the TMS effects. Second, despite data showing that the orientation of the coil with respect to the underlying cortical anatomy impacts the electric field at the target site [68], most targeting approaches do not control or account for coil orientation. With our approach, this is built in to the model at the e-field modelling step, and multiple iterations of the model are calculated at each target site with varying yaw orientations. Although real-world tests should be conducted to establish how much improvement is to be gained by controlling for coil orientation, our modelling data suggests that orientation accounts for as much variability in the simulated TMS response as does location.
Our approach is similar to other connectivity-based targeting approaches in that it assumes that pre-stimulation resting state functional connectivity (rsFC) will be predictive of treatment outcome. There are data suggesting that this may be the case. For instance, TBS induces changes in rsFC in networks that are functionally connected to stimulation site [69], and changes in connectivity are correlated with task performance related to the stimulation site [70]. Interleaved TMS/fMRI using resting statebased targeting shows that TMS induces network specific changes in activity in regions strongly vs. weakly connected to the stimulation site [21]. Similarly, baseline resting state connectivity predicts TMS-evoked BOLD responses in targets downstream from the stimulation site [71]. In terms of treatment outcome, baseline dlPFC/striatal connectivity predicts improvement in depression symptoms following a standard rTMS treatment protocol [64]. Similarly, sgACC connectivity at baseline predicted treatment outcome in a recent accelerated iTBS trial, and changes in sgACC connectivity were also associated with reduced depressive symptoms [72]. Finally, connectivity between the dlPFC and the anterior insula at baseline predicted treatment outcome following treatment with both 10 Hz and iTBS [73].
From a broader perspective (presented by Fitzqerald [2021] [74]), it is important to consider how data-driven targeting approaches can impact clinical practice. The majority of the studies evaluating the efficacy of rTMS for the treatment of depression used the 5 cm rule to target the dlPFC. Although there are issues with this approach, widespread adoption of individualized targeting in the clinic means that these patients are being given a treatment regimen that potentially differs from that which has been previously shown to be effective. In addition, it remains to be seen whether rsFC based targeting leads to substantially improved clinical outcomes. Importantly, to establish the superiority of such a method would require large samples, as the results would rely on the sub-population of patients whose individualized site is located substantially distal relative to the control site.

Strengths and limitations
The strength of the work described here stems from the generality of its applicability. The model described in this work is both site and symptom independent. Therefore, the methods described here can be used to optimize TMS targeting for a variety of disorders shown to be effectively treated with TMS [75]. This approach could be used prospectively to screen for potential TMS effects within new disorders or for novel symptom dimensions. The model is also able to make predictions at the single-subject level regarding the optimal stimulation site for a given symptom, allowing for single-subject optimization of TMS targeting. In all cases, the predictions of the model can be directly tested by applying TMS according to the model predictions. Another strength is that our results are consistent for MADRS and HAMD scores, suggesting that our model is capable of identifying stimulation sites to target depressive symptoms. However, it should be noted MADRS and HAMD scores are highly correlated across subjects (AM, r = 0.857; Control, r = 0.664), which means that the similarity in results across the two measures is likely driven by this shared variability. Accordingly, it is important not to characterize the results from these two scores as independent results, but rather different characterizations of the same underlying relationship between depression and rsFC.
The primary weakness of this study is that it is a proof of concept study, rather than a confirmatory study, meaning we do not actually administer TMS and measure the pre/post symptom changes in the current work. We understand validation of the current model is the necessary next step. Publication of the model in its initial state will allow multiple labs to test the predictions of the model, allowing for a more rapid validation/refinement of the work. To facilitate this process, we will make the code for our model publicly available.
Another potential weakness of the current work is that we used Pearson's correlations, a nondirectional measure, to model functional connectivity. However, it is likely that TMS affects upstream and downstream connections differentially [76]. Accordingly, it may be more appropriate to model the functional connections using effective connectivity, and differentially weight the connections based on whether they are upstream or downstream of the stimulation site. Of course this becomes more complicated in our model, given that the stimulation "site" is modelled as a thresholded map of the current induced across the whole brain, given a TMS coil position and orientation. Accordingly, to factor in differential effects for upstream and downstream connections, it would be necessary to generate an effective connectivity map of the entire brain. While techniques like dynamic causal modelling (DCM) can be used to generate realistic effective connectivity models with a limited number of sites and connections [77][78][79], it remains to be seen whether this approach can be reliably used to generate whole brain effective connectivity models for use at the single-subject level [80].
It should also be noted that the current analysis was limited to a single vector of points along the anterior to posterior gradient of the dlPFC. This limited search space could help explain the small effects when comparing across site. In the future we plan to expand this method to the whole brain to identify potentially novel therapeutic stimulation sites for depression and other psychiatric disorders.
Another limitation is that the model assumes that changes in connectivity tend to scale linearly with the magnitude of the electric field induced, given a particular site/orientation of stimulation. This may be another explanation for the small effects when comparing across site. This mathematical assumption likely holds for little if any in vivo applications of TMS. Indeed, it is well known that the effects of TMS on outcome variables like cortical excitability [81,82], functional connectivity [54,60], and behavior [83,84] are dependent on the stimulation parameters used [61,85]. Indeed, amplitude [86], frequency [87], pattern [88], number of sessions [6], and intervals [88] between sessions are all known to be important determinants of observed TMS effects. We understand that the linear assumption put forth in our model is a massive oversimplification, but we see this as a starting point. Accordingly, we have included a proportionality constant in our model to accommodate parameter-dependent scaling in connectivity changes. Our goal is to state the initial model as parsimoniously as possible, to allow the model to be tested with a variety of experimental and clinical TMS protocols. As more information becomes known about the relationship between rTMS parameters and connectivity changes, researchers will be able to adjust the proportionality constant of the model to suit their protocol.

CONCLUSIONS
Individualized targeting is one key to maximizing the efficacy of clinical rTMS, and mapping TMS-induced changes in rsFC to changes in clinical symptoms is critically important for this goal. The model proposed is one potential approach to accomplish this goal. The primary strength of this model is the generality of our approach, as the model proposed is both site and symptom independent. In this initial work, we show that the model makes predictions at both positive and negative control sites that are consistent with observed symptom changes following therapeutic applications of rTMS for depression [52][53][54][55][56]. The next step for this work is to test the predictions of this model using a variety of therapeutic TMS protocols. Critically, it remains to be seen whether our approach, along with other data driven, yet expensive approaches are superior to less expensive empirical approaches. What is needed is a head-to-head comparison of outcomes using these approaches in a real-world clinical sample. Such comparisons are critically needed to justify large-scale adoption of these data-driven approaches in the clinic.