Connectivity-based selection of optimal deep brain stimulation contacts: A feasibility study.

BACKGROUND: The selection of optimal deep brain stimulation (DBS) parameters is time-consuming, experience-dependent, and best suited when acute effects of stimulation can be observed (e.g., tremor reduction).
OBJECTIVES: To test the hypothesis that optimal stimulation location can be estimated based on the cortical connections of DBS contacts.
METHODS: We analyzed a cohort of 38 patients with Parkinson's disease (24 training, and 14 test cohort). Using whole-brain probabilistic tractography, we first mapped the cortical regions associated with stimulation-induced efficacy (rigidity, bradykinesia, and tremor improvement) and side effects (paresthesia, motor contractions, and visual disturbances). We then trained a support vector machine classifier to categorize DBS contacts into efficacious, defined by a therapeutic window ≥2 V (threshold for side effect minus threshold for efficacy), based on their connections with cortical regions associated with efficacy versus side effects. The connectivity-based classifications were then compared with actual stimulation contacts using receiver-operating characteristics (ROC) curves.
RESULTS: Unique cortical clusters were associated with stimulation-induced efficacy and side effects. In the training dataset, 42 of the 47 stimulation contacts were accurately classified as efficacious, with a therapeutic window of ≥3 V in 31 (66%) and between 2 and 2.9 V in 11 (24%) electrodes. This connectivity-based estimation was successfully replicated in the test cohort with similar accuracy (area under ROC = 0.83).
CONCLUSIONS: Cortical connections can predict the efficacy of DBS contacts and potentially facilitate DBS programming. The clinical utility of this paradigm in optimizing DBS outcomes should be prospectively tested, especially for directional electrodes.

Introduction

Among its other applications, deep brain stimulation (DBS) is often used as an adjunct treatment for advanced Parkinson's disease (PD).1, 2 Over the last three decades, the technique used to select stimulation parameters, or stimulation titration, has remained unchanged. Neurologists often perform a monopolar review under medication ON and OFF conditions3, which involves a step‐wise increase in the stimulation amplitude to elicit acute clinical effects (ACEs) at each DBS contact; these effects may either be therapeutic or unwanted side effects. A therapeutic window (voltage threshold for side effect minus threshold for efficacy) is then calculated to guide the final selection of stimulation contact to maximize efficacy and minimize side effects.4 Stimulation titration can take between 30 and 60 min for each hemisphere and is dependent on the experience of the clinician and active participation by the patient. Participation is required, often while the patient is experiencing unpleasant side effects, such as speech difficulty or mood changes.5 Additionally, this method is best suited for clinical effects that are easily tested and have short latency for improvement upon stimulation adjustment, but it is not practical for complex clinical phenomena with longer latency, such as improvements in gait or dystonia. Overall, this process can be highly experience‐dependent; for example, Okun et al. reported that one‐third of patients labeled as “DBS failures” improved after diligent reprogramming.6 Finally, stimulation titration will likely become more cumbersome with future use of segmented DBS leads, which have more stimulating contacts and enhanced capability for shaping the electrical field. Therefore, novel and potentially automated methods for stimulation titration based on brain connectivity are required to make the programming process easier and more efficient.

The subthalamic nucleus (STN) is the most common target for DBS in patients with PD. The STN has dense connections with the cortex,7 and similar to the functional segmentation of the basal ganglia, has motor, associative, and limbic regions.8 The dorsal STN is rich in kinesthetic cells, with connections to the primary motor cortex, supplementary motor cortex, and premotor cortex, while the medial and ventral segments of the STN are connected to the limbic and associative networks, respectively.9 In recent years, researchers have explored the mechanisms underlying the efficacy of STN DBS. In a landmark study, de Hemptinne and colleagues observed a decrease in phase‐amplitude coupling in the motor cortex coinciding with improvement in motor cardinal symptoms (rigidity, bradykinesia, tremor) after acute STN stimulation.10 This effect may be mediated through the STN's cortical connections to the motor cortex, either directly or indirectly via the thalamus,11 and likely involving the hyper‐direct pathway.12, 13 However, the translation of these findings into clinical practice faces two hurdles. (1) The regions associated with efficacy are often adjacent to those associated with side effects and therefore these must be differentiated; for example, adjacent subregions within the motor cortex may be involved in improvement of rigidity and motor contractions. (2) Distinct cortical areas may be involved in therapeutic efficacy for different clinical symptoms, and these should be identified and incorporated into a connectivity‐based titration algorithm; for example, PD can have different clinical phenotypes and patients may often experience more symptom improvement in one clinical domain with certain stimulation settings (tremor improvement more than bradykinesia, etc.). For these reasons, a connectivity‐based approach for stimulation titration will require the precise identification and modeling of distinct cortical regions associated with the plethora of stimulation‐induced clinical effects, both improvement and side effects.

We hypothesized that optimal stimulation location can be estimated based on the cortical connections of DBS contacts. To test this hypothesis, we first sought to determine which cortical regions were associated with stimulation‐induced therapeutic effects (rigidity, bradykinesia, and tremor improvements), and side effects (paresthesia, motor contractions, and visual disturbances). We analyzed acute stimulation data from 24 patients who underwent STN DBS and identified cortical connections using probabilistic tractography. We then created a connectivity‐based algorithm to estimate the therapeutic window (voltage threshold for side effect minus therapeutic effect) to classify DBS contacts into efficacious (therapeutic window ≥2 V) or nonefficacious (therapeutic window <2 V). We subsequently compared this classification to actual stimulation settings. Finally, we tested the accuracy of the estimations using an independent test cohort.

Materials and Methods

Patient selection

All patients provided informed consent and the study was approved by The Ohio State University Wexner Medical Center Ethics Board. The work described in this article has been carried out in accordance with the Code of Ethics of the World Medical Association (Declaration of Helsinki). We retrospectively analyzed clinical and imaging data from 38 patients with a clinical diagnosis of akinetic rigid PD who underwent STN DBS with Medtronic 3389 lead (Medtronic Inc., Minneapolis, MN). Prior to surgery, all participants completed a multidisciplinary evaluation with neuropsychology, radiology, movement disorders neurology, and functional neurosurgery. Intraoperative microelectrode recordings and macrostimulation were used for localization of therapeutic target. Using a wide bipolar configuration (0–3+), the therapeutic window of the entire range of DBS contacts was tested before final implantation. All patients also completed a 1‐year follow‐up.

Imaging acquisition

Each patient underwent preoperative T1 (1 mm isovoxel, MPRAGE) and diffusion‐weighted imaging (61 directions, 2 mm isovoxel, 71 axial slices encompassing the whole brain; diffusion‐weighting with b = 1000 s/mm2 was applied along 64 directions, uniformly distributed on the sphere, and one b = 0 sec volume was also acquired), using a Philips Achieva dStream 3 Tesla magnet (Philips Healthcare, Best, The Netherlands) and a padded 32‐channel birdcage coil to minimize discomfort and head motion. Postoperatively, all patients underwent a CT scan (64 channels, slice thickness 0.7 mm, no interslice gap) of the head, 4–6 weeks after implantation of the DBS electrodes. Details of the surgical procedure and postoperative imaging have been previously described.14

Recording of ACEs

Four to 6 weeks following the surgery, each patient returned for a motor assessment and programming of the DBS electrodes. The data from these visits were collected in a prospective database and retrospectively reviewed for this analysis. A movement disorder specialist increased the stimulation amplitude in a step‐wise manner while maintaining frequency at 130 Hz and pulse width at 90 μsec. The appearance of ACEs was noted separately for each implanted electrode contact. We classified ACEs into two categories (“improvement” or “side effects”). The improvement category comprised the following domains: (1) rigidity, (2) bradykinesia, and (3) tremor improvements; side effects comprised: (1) sensory, (2) motor, and (3) eye deviation or vision changes. We then combined each of the voxel maps for improvement in motor cardinal symptoms (rigidity, bradykinesia, tremor) into “improvement” and the side effects voxel maps into “side effects” to generate relevant voxel maps for training a classifier to distinguish DBS contacts into efficacious versus nonefficacious.

Imaging analysis

The images were preprocessed using customized MATLAB scripts interfacing with different software packages, including FSL,15 SPM,16 AfNI17 and lead‐DBS (Figure 1).18 Each diffusion scan was corrected for distortions due to eddy currents and head motion using affine transformations (using EDDY in FSL). Next, we computed a nonlinear transform from diffusion to MNI ICBM 2009b asymmetric template (0.5 mm) space, initialized by an affine transform to each patient's T1 image, using FLIRT and FNIRT in FSL. BEDPOSTX software (part of the FSL suite) was used for Bayesian estimation of the diffusion parameters at each voxel (3 fiber model). We determined the spatial location of each DBS contact on the postoperative CT, co‐registered the CT with the preoperative T1‐weighted MR images using Advanced Normalization Tools (ANTs) routines, and then normalized both images to the MNI ICBM 2009b asymmetric template (0.5 mm) using ANTs.

Figure 1

Schematic representation of the processing steps implemented to identify the cortical voxels associated with each clinical domain, followed by the steps required to train and test the support vector machine (SVM) classifier for the prediction of optimal combination of contact location and amplitude of stimulation.

We generated individual models of the volume of tissue activation (VTA) for each patient and stimulation pairs (voltage without clinical effect and voltage with clinical effect) associated with each ACE (FEM‐based VTA model in Lead‐DBS).19 Using the VTA as seeds, we then performed whole‐brain probabilistic tractography in each patient's native diffusion space, using Probtrackx2 in FSL. The following settings were used: 5000 samples per voxel, curvature threshold of 0.2, maximum number of steps per sample of 2000, step length of 0.5 mm, subsidiary fiber volume threshold of 0.01, loopback check, and no exclusion mask (“masking” in tractography refers to the inclusion or exclusion of specific regions in the brain, based on prior understanding of white matter pathways). To identify significant voxels associated with each clinical domain, we applied a nonparametric permutation inference approach (using paired t‐test between no‐effect/effect tractography maps with the Randomize function in FSL and threshold‐free cluster enhancement with 5000 permutations) and selected voxels with P < 0.0001 for group‐level analysis.

Identification of cortical voxels specific for each domain

To determine the differences in cortical voxels associated with different clinical domains, we performed a group‐level multivariate analysis. Voxel maps for each domain were thresholded at the maximum voxel intensity (P < 0.0001) and masked to include voxels at the cortex–white matter interface (using a 3 mm cortical ribbon dilation mask in the ICBM template space). Using 3dMVM in AfNI, “hemisphere” was selected as a between‐subject variable and “domain” as a within‐subject variable. We then performed individual linear t‐tests to compare the specific magnitude of effects associated with each clinical domain. The resulting voxel maps were then thresholded by the highest false discovery rate (FDR) value (q < 0.0001) and binarized for subsequent analysis. A high‐resolution multi‐modal cortical atlas20 was used to spatially localize significant voxels for therapeutic effects and side effects.

SVM classifier training and validation

To test the study hypothesis, we developed a classifier to distinguish the cortical clusters associated with therapeutic and side effect categories using support vector machine (SVM) analysis. SVMs are supervised learning models based on feature classification and regression analysis that are able to assign observations into different labeled categories (usually two conditions). From each individual electrode, we created VTAs at 0.5 V increments between 1 and 4 V. To help define the training space, we performed probabilistic tractography using 100 samples21 from each VTA for each subject, and used the cortical voxels associated with motor improvement and side effects as termination masks, respectively. After normalization by total number of streamlines (waytotal in FSL) and thresholding of the connectivity maps to an index of 0.01, we clustered the surviving voxels by intensity using a finite normal mixture modeling approach, which led to the extraction of nine features for each entry. An SVM algorithm based on 10‐fold repeated cross‐validation with 100 permutations was trained on the data using the Caret library in R, reaching an accuracy of 51.8% in the error matrix (95% CI = 0.5001–0.5761).22 The testing dataset used similar parameters, but used the whole‐brain tractography from each VTA without restricting it to a subset of cortical masks. Each data point consisted of a contact number and voltage specification; DBS contacts were classified as efficacious based on the span of the therapeutic window (calculated by subtracting the voltage threshold for side effect and the voltage threshold for efficacy). The efficacious contacts were further divided into two categories based on the therapeutic window: excellent (≥3 V) and good (2–2.9 V). All others were classified as “suboptimal” or “non‐efficacious”. Estimated efficacious contacts were then compared with the actual stimulation contacts used for each patient at 1‐year follow‐up. We used empirical receiver operating characteristic (ROC) curves to examine the sensitivity and specificity of this classifier at multiple cutoffs.23

SVM classifier testing

The classifier was tested on a separate cohort of 14 STN DBS patients. Eight patients were imaged using the same MRI scanner and protocol of the training cohort. Six additional patients were imaged using different scanners with comparable structural and diffusion‐weighted protocols. (see Data S1 for the DTI sequence parameters used in this cohort).

Results

Demographics and clinical measures

The mean age of the patients at the time of surgery was 65 years (SD: 8 years). All patients showed significant improvement in motor cardinal symptoms (rigidity, bradykinesia, tremor) after DBS (preoperative OFF medication, UPDRS‐III: 35.1 ± 8.8; postoperative OFF medication and ON DBS, UPDRS‐III: 16.1 ± 8.2), corresponding with a significant and consistent reduction in Levodopa‐equivalent dose (preoperative: 1590.8 ± 642.2 mg vs. postoperative: 730.7 ± 409.4 mg). Twenty‐four patients received DBS implantation in 38 hemispheres (14 bilateral and 10 unilateral).

Cortical clusters associated with efficacy and side effects

Distinct cortical clusters are associated with stimulation‐induced improvements in motor symptoms (rigidity, bradykinesia, tremor) versus side effects

Figure 2 displays the cortical areas associated with stimulation‐induced improvements in motor symptoms: superior frontal gyrus, dorsolateral prefrontal cortex (corresponding atlas notation 8Ad, 9a, 9p), supplementary motor cortex (atlas notation 6 m), and medial prefrontal (atlas notation 8BL). The voxel clusters in the premotor and supplementary motor cortex (Brodmann area 6) were significantly associated with stimulation induced motor improvement than involuntary muscle contractions. Overall, while comparing stimulation‐induced efficacy with side effects (sensory, motor, and vision changes), the voxels in the superior frontal gyrus and frontal eye field (atlas areas 8BL), bilaterally, were associated with the efficacy domains20 (Table S1). In contrast, the middle and inferior frontal gyrus (middle and inferior premotor cortex; atlas notation area 9 m and i6‐8) and frontopolar cortex (atlas notation p47r) were associated with side effects.

Figure 2

Summary of cortical voxels associated with efficacy (A; improvements in rigidity, bradykinesia, and tremor) and side effects (B; sensory, motor, and visual). The voxels are overlaid on the Conte69 standard surface in left and right lateral hemispheres. Below each hemispheric representation are waterfall plots displaying the cortical areas with the highest connections (number of voxels connected), using the cortical parcellation and nomenclature published previously in Glasser et al. (20).

Distinct cortical clusters are associated with stimulation‐induced improvements in motor symptoms (rigidity, bradykinesia, tremor)

In addition to the shared cortical regions, improvements in rigidity, bradykinesia, and tremor were also associated with distinct cortical connections (mixed effect model for each individual domain, with FDR correction, q < 0.0001). The complete results are presented in Table S1. The motor cortex (atlas area 4) was associated with improvement in both bradykinesia and tremor, while the premotor cortex and medial prefrontal cortex (atlas area 9 m) were uniquely associated with tremor and bradykinesia, respectively. Both rigidity and bradykinesia improvement, but not tremor improvement, were associated with connections to the middle frontal gyrus (area 8Ad).

Estimation of efficacious contacts

To train the SVM classifier, 794 out 1064 entries were classified as “efficacious” and 270 as “side effects”, based on nine selected features. The graph visualizing the importance of individual variables estimated by a learning vector quantization model24 is available in Figure S1.

Out of 152 implanted contacts, 137 (90.1%) were estimated to have a therapeutic window between 2 and 3 V. At 1‐year follow‐up, 47 stimulation contacts were being used in 38 implanted hemispheres. The SVM algorithm classified 42 out of these 47 stimulation contacts (89.4%) as efficacious with a therapeutic window of either ≥3 V (n = 31) or 2–2.9 V (n = 11). Five stimulation contacts were incorrectly classified as either nonefficacious or within the therapeutic window of 1–2 V. The ROC curve is displayed in Figure 3; the AUC (area under curve) is 0.8506 (95% CI = 0.7026–0.9987). The sensitivity of this estimation was 89.4% (95% CI = 76.9–96.5%) and the specificity 5% (95% CI = 1.6–11.3%). The contacts estimated for each patient are outlined in Table 1.

Figure 3

The ROC curve comparing the sensitivity and 1‐specificity of the SVM classifier in the 24 patients (training) cohort.

Table 1

The estimated therapeutic window in the ≥3 V and 2–2.9 V categories for the training cohort.

Subject ID	Estimated therapeutic window
	Left hemisphere								Right hemisphere
	≥3 V				2–2.9 V				≥ 3V				2–2.9 V
1	0	1	2	3					0	3			2	1
2					3	0	1		1	3			2	0
3	0	1	2	3					0	1	3		2
4	0	1	3	2
5	0	1	3		2				0	1	2	3
6	0	1	2	3					0	1	2	3
7	0	1	2	3
8	1	2	3
9					1	0	2	3	0	1	2	3
10	0	1	2	3					1	2	0	3
11	0	1	2		3
12	0	3	1		2
13	1								1	2			3
14	1	2	3		0
15
16	0	2											3
17									0	2			3	1
18	0	3	1		2
19	0	1	2	3					0	1	2	3
20	1	2	0		3				0	1	2		3
21	0	1	3		2				0	1	2	3
22	1	3			0	2			1	2	3		0
23									0	3			2	1
24	0	1	3	2					0	1	2	3

The actual stimulation contacts at 1 year are highlighted in green (≥3 V) and orange (2–2.9 V categories).

In the test cohort, the classifier was able to match the long‐term stimulation settings for 19 out of 37 stimulated contacts – Table S3). The ROC curve is displayed in Figure S3; the AUC is 0.8333.

Discussion

Here, we demonstrate the feasibility of identifying efficacious STN DBS contacts based on diffusion MRI (dMRI) connectivity in PD patients. Using probabilistic tractography, we first created cortical connectivity maps associated with improvements and side effects. An automated algorithm was then created to identify efficacious stimulation location based on the preferential connectivity of DBS contacts to the cortical regions of interest.

Cortical connectivity patterns

Cortical regions associated with therapeutic effects of STN DBS have previously been investigated using dMRI.12 In line with these reports, we also observed that improvements in bradykinesia, rigidity, and tremor were associated with preferential white matter connectivity to motor areas in the premotor, supplementary motor, and motor cortex. We were also able to distinguish clusters associated with side effects from the neighboring regions associated with efficacy; for example, adjacent voxels within the motor cortex were associated with stimulation‐induced efficacy and motor side effects. Similarly, paresthesia was associated with voxel clusters in the sensory cortex. In addition, we found unique cortical regions associated with efficacious clinical domains; for example, tremor and bradykinesia were associated with significant connections to the primary motor cortex (M1). Similar findings were recently reported by Horn et al.19 The pathophysiology of resting tremor is thought to be distinct from levodopa‐responsive bradykinesia and rigidity mainly because the associative motor areas, instead of the primary motor (M1), are more involved with execution than rest.25 Additionally, the involvement of nonmotor areas including the lateral frontal, paracentral, inferior‐parietal, and parieto‐occipital regions, has been previously linked to the severity of bradykinesia.26

It is important to emphasize that whole‐brain probabilistic tractography analysis has high sensitivity but poor specificity, implying that some cortical clusters identified in this analysis may be false positives. One strategy to avoid this involves restricting the target masks to predefined regions of interest; for example, based on the proposed existence of the hyper‐direct pathway between the STN and the motor cortex, the analysis can be restricted to the motor and premotor cortices.27 We identified cortical regions using whole‐brain probabilistic tractography in an effort to increase sensitivity and to maximize the likelihood of identifying all potential cortical voxels13 without restricting the analysis based on the known existence of specific white matter pathways.28 Similarly, we did not limit the tractography analysis to one hemisphere because there may be some overlap of white matter pathways at the brainstem level, and their modulation may result in DBS‐induced ACEs.29 To increase the specificity, we used nonparametric statistics to threshold the results, based on the established statistical methods from functional connectivity literature.30 To predict therapeutic contacts, we aimed to eliminate false positive results and therefore restricted target masks to voxels significantly associated with clinical efficacy or side effects during the SVM training.

Implications for stimulation titration

Stimulation titration may be challenging, even in patients with optimally‐placed electrodes. This aspect will become even more important with the wider use of multi‐contact directional DBS leads,31 which exponentially increase the number of available options for testing and programming. Investigations into the modeling of electrical fields associated with DBS electrodes have improved visualization and identification of potential neural elements associated with therapeutic effects of DBS. Assuming that each patient has a unique connectivity fingerprint,32 we used probabilistic tractography to formulate an alternative approach for selecting optimal stimulation amplitude & contact. We trained an SVM classifier to categorize connectivity patterns associated with either “improvement” or “side effect”. We then tested the classifier to label whole‐brain connectivity maps from each simulated VTA. Future strategies to improve the sensitivity and specificity may include analysis of connectivity to subcortical structures, clinical information (e.g., levodopa response), and potentially using advanced machine learning algorithms. Recent advances in GPU‐based statistical programs improve computational time, which makes it possible to accomplish extensive repetitive tasks in a relatively short amount of time. Using a dedicated server with a 16‐core CPU and two GPU graphics cards, the prediction of optimal stimulation settings in a bilateral STN DBS cases takes 30–45 min from the stages of imaging import and normalization to probabilistic tractography and building the SVM feature table required to classify the contacts.

Limitations

This was a retrospective analysis using data from 38 patients with PD. Given the constraints of acquisition times in a clinical population, we chose a dMRI sequence with 2 × 2 × 2 mm voxel size. However, this scheme limits the number of voxels sampled within STN, for example, we will sample 12 voxels in the STN assuming its size to be 4 × 4 × 6 mm. The feasibility of improving voxel resolution to 1.5 mm iso‐voxel should be explored in future studies to further increase the number of voxels sampled in STN. Despite the rigorous criteria for patient selection and DBS programming, the results may not be immediately generalizable due to the limited sample size. Additionally, we did not analyze the degree of improvement or severity of side effects. The model used for VTA estimation is based on a FEM model that does not accurately capture the different impedances of cellular populations in white matter versus gray matter.18 It is also worth mentioning that electrode impedances change over time, affecting the stimulation volume.33 The proposed solution does not take into consideration the spatial location of the DBS electrode relative to the STN boundaries and the sampling of LFPs from each contact, that can potentially improve the predictive value of the model. We limited our analysis to only cortical targets because probabilistic tractography is limited in its ability to accurately capture local connectivity in subcortical gray matter structures.34 A larger sample size is desirable to rank the contacts based on their likelihood for individual symptom improvement (e.g., tremor more than rigidity or bradykinesia). While this algorithm is able to predict optimal stimulation contacts, the prediction of other stimulation parameters like stimulation frequency and pulse width will require further investigation. The processing time needed for prediction currently, can potentially be reduced to a manageable scale by the use of advanced machine learning algorithms.

Author Contributions

1. Research Project: A. Conception, B. Organization, C. Execution

FS, VK: A, B, C

AR, YN: B

QR: C

2. Statistical Analysis: A. Design, B. Execution, C. Review and Critique

FS, VK, QR, YN: A, B, C

3. Manuscript Preparation: A. Writing of the first draft, B. Review and Critique

FS, VK, YN: A, B

PD, BC, MD, MK, AR: B

Conflict of Interest

The authors have no conflicts of interest to declare while working on this manuscript.

Figure S1. The graph visualizing the importance of individual variables in the SVM model estimated by a learning vector quantization approach.

Click here for additional data file.

Figure S2. The group parcellation of the subthalamic nucleus (A–C) based on their connections with cortical masks associated with rigidity (violet), bradykinesia (yellow), rigidity (green), involuntary motor contractions (red), paresthesias (blue), and visual disturbances (sea green). The parcellation with the actual location of STN (pink) as identified with intraoperative physiology (D–F).

Click here for additional data file.

Figure S3. The ROC curve comparing the sensitivity and 1‐specificity of the SVM classifier in the 14 patients (testing) cohort.

Click here for additional data file.

Table S1. Full list of cortical areas (in both Glasser et al. and Brodmann naming convention) associated with each clinical domain.

Click here for additional data file.

Table S2. The MNI coordinates of the clusters in the subthalamic nucleus.

Click here for additional data file.

Table S3. Test cohort results. For each patient the correctly predicted option for each side is highlighted in bold.

Click here for additional data file.

Data S1. Supplementary methods.

Click here for additional data file.