Probabilistic Mapping Reveals Optimal Stimulation Site in Essential Tremor

The objective of this study was to obtain individual clinical and neuroimaging data of patients undergoing deep brain stimulation (DBS) for essential tremor (ET) from 5 different European centers to identify predictors of outcome and to identify an optimal stimulation site.

substrate of stimulation-induced tremor alleviation. [5][6][7] These findings could implicate that instead of a confined anatomic (sub)area, stimulation at any point along the CTT network will elicit optimal tremor suppression. Depending on the intended target site of stimulation, reported follow-up, and outcome measure, thalamic and subthalamic DBS have been reported to achieve an on average 48% to 73.8% tremor reduction. 8,9 Apart from the observed variability of outcomes between studies, there is also a considerable variability across studies about suboptimal and nonresponders that have ranked up to 22%. Furthermore, inconsistent and partly contradictory results were published regarding predictors of outcomes, such as the preoperative tremor severity, age, and stereotactic coordinates of the stimulating electrode. [8][9][10] Hence, after almost 30 years of DBS for ET, it remains unclear if there are solid predictors of outcome and if there is an optimal stimulation site. Previous studies were mostly based on smaller and monocentric series of patients with presumably lower variability in electrode locations. However, a certain degree of variability in electrode placement and measured outcomes is a prerequisite for studies that aim to investigate relationships between stimulation sites and clinical outcomes. To this end, analysis of aggregated data from multiple centers with differing electrode targeting approaches would be appropriate. In this study, we gathered clinical and neuroimaging data of a large cohort of patients that underwent DBS for ET in 5 different centers to identify robust predictors of outcome and to identify an optimal stimulation site by applying probabilistic mapping.

Patients
This study retrospectively enrolled datasets of 119 patients with ET syndrome treated chronically with unilateral or bilateral DBS and operated at 5 different European DBS centers. Diagnosis of ET was assessed by specialized movement disorder neurologists according to recommended guidelines and indication for surgery was discussed at multidisciplinary boards. 11 Large parts of the data were recently published (13 patients from a randomized controlled crossover trial 12 and 81 from 4 retrospective studies [13][14][15][16]. Unpublished data from 25 additional patients were included where full datasets were available. Individual datasets were included in the present study if they contained: (1) preoperative tremor scores of validated clinical rating scales, (2) available baseline covariates (sex, age, and disease duration), (3) 12 months of postoperative tremor scores with stimulation on, (4) pre-and postoperative neuroimaging allowing reconstruction of lead location or individual anterior commissure-posterior commissure coordinates of the reconstructed stimulating electrode, and (4) stimulation parameters at 12 months following surgery. Exclusion criteria were a diagnosis of other tremor syndromes (Holmes tremor, dystonic tremor, and Parkinson tremor) or a magnetic resonance imaging (MRI)-verified lesion in the cerebellum, thalamus, or brain-stem. The study was conducted in accordance with the Declaration of Helsinki and approved by the institutional review board of the University Bern (KEK 2018-00841).

Surgical Procedure and Clinical Evaluation
In all patients, electrodes (model 3389 or 3387, Medtronic; Boston Cartesia Vercise directional leads, Boston Scientific, or St. Jude 6148, St. Jude Medical) were implanted stereotactically into the Vim or PSA. The neurostimulation parameters were programmed according to best clinical practice by the local DBS neurologist. Pre-operative "medication-off" and postoperative "stimulation-on, medication-off" tremor severities were assessed "nonblinded" based on the Fahn-Tolosa-Marin tremor rating scale (TRS) in 4 centers and based on the Bain's TRS in 1 center by a specialized movement disorder neurologists at 12 months postoperatively. 17,18

Study Outcomes
The primary outcome measure was the percentage reduction of the part A subscore for upper and lower extremities of the TRS or Bain's TRS per hemibody at 12 months following DBS lead implantation.

Lead Reconstruction and Estimation of Stimulation Volumes
DBS electrodes were reconstructed using the Lead-DBS toolbox (version 2.3.2) in MATLAB 2019b (The MathWorks, Natick, MA). 19 Postoperative computed tomography (CT) or MRI scans were linearly coregistered to the pre-operative MRI using advanced normalization tools (ANTs) with brain shift correction as implemented in lead-DBS. Multispectral normalization to ICBM 2009b MNI space was carried out by applying the ANTs SyN Diffeomorphic Mapping. 20 This method was shown to segment the subthalamic nucleus (STN) region with high precision comparable to manual expert segmentations in a recent comparative study. 21 DBS electrodes were automatically prereconstructed using the phantom-validated and fully automated physical, affective, cognitive, environmental, and relationship (PACER) method. 22 All reconstructed electrodes were individually checked visually for plausibility by examining the postoperative electrode position on the original postoperative image in relation to reliably detectable landmarks by a neurosurgical fellow (author S.B.) and a senior functional neurosurgeon (author A.N.). In a minor number of cases (11 patients and 19 electrodes), lead-DBS failed to detect electrodes automatically or provided an aberrant reconstruction far outside the thalamosubthalamic region. In these specific cases, we compared the reconstruction result with the electrode position from the postoperative MRI or CT scan and manually corrected the electrode position according to the artifact based on the postoperative images in lead-DBS. For segmented electrodes, the orientation was determined with the directional orientation detection algorithm that was validated in phantom models and clinical studies before. 23,24 All volumes of tissue activation (VTAs) were estimated with the SimBio/Fieldtrip pipeline directly in template space with the DISTAL atlas. 25 We applied the same VTA model and settings for all leads to facilitate computation. Default values were used for white matter and grey matter conductivities, 0.14 and 0.33 S/m, respectively, and a threshold of 0.2 V/mm.

Stimulation Map Construction
We computed a probabilistic stimulation map by aggregating individual stimulation volumes with their associated tremor improvement of the contralateral hemibody at 12 months. All VTAs were estimated in MNI space and left hemispheric VTAs were nonlinearly warped to the right hemisphere.
First, we calculated an N-image of the cohort by aggregating VTAs across the entire cohort. The N-image represents a heatmap and can be thresholded according to the number of times (n) each voxel in the map gets activated. To display the entire stimulation area of the cohort, we calculated the N-image to illustrate all voxels that were at least activated 4 times (representing 10% of small cohorts).
Second, we calculated a mean improvement map. To this end, each voxel was assigned the mean tremor improvement (part A) of all associated stimulation volumes activating this voxel. Each voxel had to be activated by at least 4 different stimulation volumes to be included in the analysis.
Third, we computed a significant improvement map adapted from Reich et al. 26 In a first approach and closely akin to the method by Reich et al, the tremor improvement values associated with a given voxel were tested against the remaining improvement values not associated with this voxel by applying a 2-sided t test. With a significance level of 0.05, we yielded a better improvement cluster as well as a worse improvement cluster. However, both clusters were relatively small and did not pass multiple testing correction with a false discovery rate of 0.05 (data not shown). The reasons for this probably comprise the non-normally distributed outcome data with a disproportionate number of excellent responders as well as the inhomogeneous stimulation map with a predominant number of voxels that get activated only a few times by different VTAs resembling a lognormal distribution. In a second approach, we shifted from a binary better-worse division to a division into suboptimal, good, and excellent improvement clusters according to the 33%, 66%, and 100% percentile. We intended to separate the suboptimal improvement cluster (corresponding to voxels with a <50% tremor reduction) from a combined good and excellent improvement cluster (representing voxels with >50% tremor reduction). Therefore, we applied a right-tailed Wilcoxon signed rank test. The null hypothesis was that the tremor improvement values associated with a given voxel were from a distribution with a median m. We chose m to be the 33% percentile of all tremor improvement values. The alternate hypothesis was therefore that the associated tremor improvement values of a given voxel were from a distribution with a median >33% percentile. With a significance level of 0.05, we identified 2 different clusters and multiple testing correction was applied with a false discovery rate of 0.05. This resulted in a significant improvement map highlighting a statistically significant improvement cluster or statistically significant stimulation "sweet spot." The suboptimal improvement cluster did not pass multiple testing correction and was removed from further analysis.
The N-image, mean improvement map together with their suboptimal and excellent responder clusters and significant improvement cluster were displayed in the 3-dimensional DIS-TAL human brain atlas as well as superimposed onto the 2-dimensional Schaltenbrand and Wahren stereotactic human brain atlas. 27

Structural Connectivity Analysis
The mean improvement map was displayed in relation to an average normative CTT based on high-quality diffusion imaging data of 32 healthy subjects from the Human Connectome Project. The method was recently described in detail in Dembek et al. 6 In summary, 62 CTTs were calculated using BEDPOSTX for estimating the local fiber orientation and PROBTRACKX2 for probabilistic tractography in FSL. 28 The contralateral dentate nucleus was chosen as the seed region, whereas the contralateral superior cerebellar peduncle, the ipsilateral red nucleus, and the ipsilateral precentral gyrus served as waypoints. Each individual track frequency map was transformed into a track probability map into the MNI space and probability values per voxel were averaged across all subjects.

Outcome Prediction Model
Finally, we validated the significant improvement map. To this end, we applied a stratified 10-fold cross-validation (ie, the data set of 213 stimulation volumes was randomly split into 90% training data [192 stimulation volumes] and 10% test data [21 stimulation volumes]). The significant improvement map or sweet spot was computed with the training data. With the test data, we calculated the Euclidean distance from the centroid of a test stimulation volume to the centroid of the sweet spot. We also calculated the overlap volume between the test stimulation volume and the CTT. This was repeated 10 times (ie, each stimulation volume was a test volume once). Of note, the sweet spot was spatially stable across the stratified 10-fold cross-validation. A stratified 5-fold cross-validation did not pass multiple testing correction to calculate the significant improvement map. Here, 80% of the data set, or 170 stimulation volumes, were not sufficient.
First, we intended to predict continuous tremor improvement. The factors "distance from VTA to the sweet spot" and the "probabilistic overlap volume between VTA and CTT" were fed as predictor variables into a support vector machine regression model. Tremor improvement was the response variable. Our assumption was that both proximity to the sweet spot and activation of the CTT would be meaningful predictors. Of note, a VTA may be close to the sweet spot but not activate the CTT and such a VTA would be assumed to result in less tremor improvement than a VTA with equal distance to the sweet spot but that activates the CTT. With respect to the model, we used a nonlinear regression model with a Gaussian Kernel function. Conceptually, this transforms the 2 predictors "distance" and "overlap volume" to a higher dimensional space. In that space, the algorithm then computes a linear function that is close to the response variable, while being as flat as possible (MATLAB command fitrsvm).
Second, we pursued another approach and intended to classify tremor improvement. Again, the combination of factors distance and probabilistic overlap volume were used and fed into a classification ensemble (MATLAB command fitcensemble with automatic hyperparameter optimization). The intention was to predict whether stimulation volumes would result in greater or less than 50% tremor improvement.

Statistical Analysis
Data were analyzed by applying descriptive nonparametric statistics after normality testing and visual inspection of the QQ-plots using SigmaPlot (Systat Software, San Jose, CA) and MATLAB (The MathWorks, Natick, MA) for baseline characteristics, tremor scores, and outcome data. Nested linear mixed-effect models were used to test for the influence of different baseline covariates (sex, age, pre-operative tremor score, center, total electrical energy delivered as fixed factor, and electrodes as random factor) on outcome (percentage change of part A tremor subscores per hemisphere/ contralateral hemibody). For statistical testing of the mean improvement map, we applied a one-sided Wilcoxon signed rank test (see above) and multiple testing correction with false discovery rate of 0.05. Data are presented as mean AE 95% confidence interval (CI) if not indicated otherwise. Statistical tests were 2-tailed and a p value <0.05 was considered statistically significant.

Clinical Data
A total number of 119 patients from 5 European centers were included in the study. Table 1 shows the baseline demographic data for each center's patient cohort. In all patients, the part A subscore was available after 12 months follow-up, which corresponded to 237 hemispheres. The mean tremor reduction (part A) per hemisphere across all centers was 64% (61-68%, 95% CI). Scatter plots of outcomes per center are shown in Figure 1. The results of the nested linear mixed effect model with age, gender, pre-operative tremor intensity, total electrical energy delivered by DBS and center as fixed factors, and patient and lead nested in the patient as a random factor yielded only the 2 factors "pre-operative tremor intensity" and "center" as significantly associated with outcome ( Table 2). Individual tremor improvement varied substantially within the entire cohort (range = À25 to 100% tremor improvement). The outcome distribution of the entire cohort was leftskewed. The 33% percentile corresponded to a tremor improvement of 50% and the 67% percentile corresponding to an 80% tremor reduction.

Stimulation Map
Electrodes could be reconstructed in 107 patients (213 hemispheres). Twelve patients (16 hemispheres) had  to be excluded from the image analysis as insufficient data quality of the postoperative CT or MRI scan did not allow a sufficiently precise detection of the electrode. Figure 2 displays the spatial distribution of the implanted leads in the Vim and PSA. The stimulation area covers large parts of the motor thalamus (beyond the Vim) and large parts of the PSA lateral to the red nucleus (RN), covering the caudal zona incerta (cZI) and extending into the entire posterior STN. As the results of the linear mixed effect model suggest the factor "center" to be significantly associated with outcome (with patients from the Bern cohort having a significantly worse outcome compared to  the other centers), we compared the electrode implantation sites among centers. Based on both visual inspection and statistical testing, the electrodes from the Bern cohort are on average placed more medial compared to the other centers (2.2 AE 1.45 mm, p < 0.001 analysis of variance [ANOVA]). For this reason, we decided not to correct for the factor "center" when constructing the stimulation map and outcome prediction models as we assume the difference between centers to result from different electrode placements. The mean improvement map represents voxels that are associated with a certain clinical outcome. Displaying the whole map is visually not very informative, as areas associated with excellent tremor reduction are covered by areas of suboptimal tremor reduction, we thresholded the map according to our statistical approach by applying a partition into 3 clusters according to the 33%, 67%, and 100% percentiles. We display the suboptimal responder cluster corresponding to a <33% tremor reduction and an excellent responder cluster corresponding to a >67% tremor reduction along with the surface volume of the entire unthresholded map in Figure 3. Whereas the area of excellent tremor reduction extends from the posterior aspect of the PSA to the central part of the Vim, the suboptimal responder voxels are located in the anterior and medial aspect of the PSA and extend to the ventralis oralis posterior nucleus and the internal capsule laterally. The excellent responder cluster aligns visually well with the tractography-based CTT (see Fig 4). Furthermore, we display the location of the identified significant good responder cluster that projects onto the margin of the posterior aspect of the PSA and Vim. The mean position of implanted DBS electrodes was X = 12.66, Y = À15.02, and Z = À3. 40. In comparison, the centroid coordinates of the statistically significant good responder cluster were X = 13.63, Y = À15.13, and Z = À2.58 mm in the MNI space. A projection of the stimulation map onto the stereotactic atlas of Schaltenbrand and Wahren is shown in Figure 5.

Structural Connectivity
We displayed the suboptimal and excellent responder clusters as well as the significant good responder cluster in relation to the averaged normative CTT template (see Fig 4). The CTT passed through the significant good responder cluster. This visual impression was confirmed by statistical testing. Voxels within the significant good responder cluster were more likely to contain the CTT as expressed by higher tract probability values compared to the rest of the map (0.50 AE 0.13 vs 0.21 AE 0.08; p < 0.0001, t test).

Outcome Prediction Model
Finally, we evaluated the predictive value of this map by applying a stratified 10-fold cross-validation. The clinical outcome predicted from the significant improvement map by applying the support vector machine was associated with observed clinical outcome (mean R 2 = 0.14, p = 0.02). Thus, the model estimations explained 14% of the variance of tremor score improvement. Using the distance to the sweet spot alone resulted in an R 2 of 0.06, while using the probabilistic overlap volume alone yielded an R 2 of À0.17 (ie, distance to the sweet spot and probabilistic overlap volume alone would not be a good predictor). In comparison, an ordinary linear regression model with both distance and probabilistic overlap volume would yield an R 2 of 0.037, underlining the better performance of the support vector machine. When applying a classification ensemble on the distance and overlap volume, we could correctly classify good and suboptimal responders (corresponding to >50% and <50% tremor reduction, respectively) with a sensitivity of 89% a specificity of 57%. The overall accuracy was 72%.

Discussion
To our knowledge, this is the largest multicenter study to integrate and compare clinical and neuroimaging data of patients with ET undergoing DBS. There was a considerable anatomic dispersion of implanted leads across included patients, which was reflected by an overall widespread area of stimulation covering large parts of the motor thalamus and the PSA. By applying voxel-wise statistical testing we could point out a "sweet spot" area of stimulation in the posterior part of the PSA and the inferior aspect of the Vim. The area of this "good responder cluster" coincided with the area of highest likelihood to contain the tractography-based normative CTT. Taken together, these findings suggest that stimulation of the posterior PSA and the Vim along the course of the CTT is associated with the highest likelihood of tremor suppression. The anatomic location of individual stimulation volumes was the overall best predictor of stimulation-induced tremor suppression with a sensitivity of 89% and a specificity of 57% to correctly identify responders with a >50% tremor reduction based on a leave-one-out cross validation.
The 5 included centers placed their electrodes by slightly different targeting areas and approaches. Of note, the identified statistically significant improvement cluster (stimulation sweet spot) is closely related to the mean position of the implanted DBS electrodes (stimulation hotspot). The results can be interpreted in different ways. First, the exact location and extent of the identified significant improvement cluster is too restricted as it depends on the distribution of data and our chosen statistical approach by applying voxel-wise statistical testing and multiple testing correction. Areas of the map that contain more datapoints (which is the case at the stimulation hotspot) will more likely yield significant results and survive multiple testing correction. The fact that the stimulation area, which is associated with a good outcome, is more widespread than the identified significant stimulation sweet spot is in favor of this explanation. On the other hand, the identified significant improvement cluster could indeed reflect the "true" sweet spot and different the different centers just target around this optimal stimulation point based on years of clinical experience and targeting optimization. In line with this view, our findings are based on chronic stimulation settings that reflect empirical stimulation parameter adjustment to optimize the stimulation outcome.

Anatomic and Functional Considerations
The CTT projects from the deep cerebellar nuclei mainly to the contralateral thalamus. [29][30][31] The observation that several previous studies of DBS for ET reported satisfactory tremor reduction upon stimulation of slightly different thalamic and subthalamic targets including the PSA, cZI, or prelemniscal radiations led to the formulation that the CTT might be the underlying neuroanatomic substrate stimulation of which would suppress tremor. 3,5,7,32,33 Since then, multiple studies from different groups studied the relationship among electrode locations, VTAs, and the tractography-based CTT or functional connectivity patterns of good and poor responders but came to different and sometimes contradictory conclusions. For instance, Akram et al and Al-Fatly et al identified a thalamic sweet spot of optimal connectivity to the contralateral cerebellar dentate nucleus and the ipsilateral primary motor cortex (M1), whereas Middlebrooks identified a thalamic region that was connected to the supplementary motor area (SMA) and premotor area to be correlated with better tremor reduction. 13,34,35 The results of the present work indicate that the area of stimulation which corresponds to efficient tremor reduction coincides with a high probability to contain the CTT, which is structurally connected to M1 and the contralateral dentate nucleus of the cerebellum (ie, support the view of Akram's and Al-Fatly's work as well as the one by Dembek et al). 6 The present findings of the location of the probabilistic stimulation sweet spot at the transition between the ventral Vim and the PSA are in line with previous findings from an intra-operative electrophysiological study by Milosevic et al. 36 In their intra-operative study, the authors investigated the effects of high-frequency microstimulation on both neuronal firing and tremor suppression simultaneously by 2 closely placed microelectrodes. They found that 200 Hz high-frequency stimulation was significantly more effective to reduce tremor reduction and spontaneous cell firing in thalamic neurons compared to 100 Hz stimulation. Notably, they observed that the most ventroposterior stimulation sites had the best effect on tremors. The authors concluded that thalamic neuronal inhibition seems necessary for tremor reduction and may function in effect as a thalamic filter to uncouple thalamocortical from corticospinal reflex loops.
Although there is no question about the existence of cerebellothalamic connections from the deep cerebellar nuclei to the motor thalamus and their relevance in modulating movement, there is some controversy about which deep cerebellar nucleus is predominantly involved in tremor genesis. Based on anatomic tract tracer studies both the dentate nucleus as well as the interposed nuclei were shown to be connected to the motor thalamus and frontal regions including M1 and the SMA depending on the applied methodology and examined species. 30,31,[37][38][39] The present study results suggest that high-frequency stimulation of cerebellothalamic fibers from the dentate nucleus to the Vim suppress tremors based on probabilistic tractography findings. However, we do not claim that fibers from the interposed nuclei also contribute to this effect. A definitive answer to this controversy can only be provided by well-designed post-mortem anatomic studies in humans.
The presented results point toward a possible causal relationship of CTT stimulation for tremor suppression. However, an involvement of the cZI for tremor genesis and stimulation-induced suppression cannot be ruled out. The cZI is bounded by the medial lemniscus posteriorly, the ascending cerebellothalamic fibers in the prelemniscal radiation anteriorly and the superior part of the STN laterally. As already pointed out by Plaha, the cZI provides a unique GABAergic link between the basal ganglia output nuclei and the cerebellothalamocortical loop which places it in a key position to transmit synchronized oscillations generating tremor into these loops. 33 The close anatomic proximity of the CTT and the cZI makes it difficult, if not impossible to segregate DBS induced stimulation effects on these two structures with current neuroimaging and VTA-modeling approaches.

Limitations
Limitations of this study are the retrospective nature of the magnitude of used primary outcome data to calculate the probabilistic stimulation map and the two different outcome scales used across centers. However, both clinical TRSs are validated and by using the percentage change of tremor, outcome data become normalized and reflect tremor improvement independent from the scale used. Second, the algorithm for image analysis, including image normalization, lead reconstruction, and electrode orientation is inherently prone to error. Although, most of the applied methods have been validated in phantom models before, without histological confirmation, DBS electrode reconstruction always remains presumptive. These inherent errors are very likely to limit the overall predictive value of the applied model. Third, we used a VTA model that assumes an isotropic environment, although the leads had been implanted in a region with moderate anisotropy. Diffusion weighted imaging and an expansion of the VTA model would be necessary to provide better patientspecific fiber tracking and to accommodate anisotropy. Moreover, each algorithm dedicated to co-registration and normalization involves inherent errors and impact the overall calculation precision. Fourth, we used normative connectome data to estimate the average CTT. Although these normative connectome atlases do not represent patient-specific connectivity, they in turn have the benefit of high signal-to-noise ratios. Furthermore, the results of a recent study by Wang et al put into perspective the limitations of normative versus patient-specific connectivity analysis. 40 In addition, we did not integrate side effects into our probabilistic modeling approach. One major reason is the retrospective design of the study and missing data. However, it was regular clinical practice across all centers to titrate and optimize the stimulation parameters by balancing the best possible tremor suppression while avoiding limiting stimulation-induced side effects. It is therefore reasonable to assume that the probabilistic stimulation map in a certain way contains indirect information about side effects. Nevertheless, future research needs to focus on objectively and quantitatively assessed side effects.
Last, we tested different outcome-prediction models based on different assumptions and variables to get to the final results that were not controlled or corrected for these multiple approaches. Still, the optimized algorithm underlying the outcome-prediction model can only explain 14% of the variance of the observed outcome if we seek to predict individual percentage tremor reduction. We suspect that the overall widespread stimulation area and the noisy nature of both clinical outcome and imaging data are the main reason that our applied models fail to predict individual continuous tremor improvement more accurately. Because of the low R-squared values, we further tested a classification with a random forest to predict a binarized outcome instead. The results of the cross-validation and their predictive value need to be interpreted with caution and the presented results are mainly descriptive until confirmed in independent datasets, ideally through prospective testing.