Movement-responsive deep brain stimulation for Parkinson’s Disease using a remotely optimized neural decoder

Embedded aDBS functionality and algorithm tuning overview

One subject with PD who had previously received bilateral DBS implants with embedded adaptive capabilities (Medtronic RC+S, Figure 1C) was recruited for the study. The hardware in each hemisphere consisted of a depth lead for stimulating and sensing in the subthalamic nucleus (STN), a sensorimotor electrocorticography (ECoG) strip for sensing over the precentral (PreC, motor) and postcentral (PostC, somatosensory) cortices, and a neurostimulator for driving adaptive stimulation using LFP’s from the three sensed brain areas. The onboard algorithm guiding stimulation used linear discriminant analysis with additional layered operations for controlling the state transition dynamics (Sellers et al., 2021). LFP’s were first transformed into power bands with programmable frequency limits and smoothing properties, then up to four power bands were linearly combined and passed through a series of threshold-based logic operations to define the algorithm state. Our aim was to optimize the algorithm parameters in each stimulator for decoding movement of the contralateral arms and delivering targeted stimulation increases during their respective movements (Figure 1A-B).

A remote data collection platform (Strandquist et al., 2023) was deployed in the participant’s home office to collect model training data and evaluate algorithm performance. In addition to neural data, we recorded accelerometry from bilateral Apple Watches (thresholded to label as moving or stationary), multi-view video recordings, and keylogging from the participant’s home computer (Figure 1C). Five days of data were collected for model training and validation, while a sixth day was dedicated for model testing. Each recording day lasted approximately 25 minutes, during which the participant performed a self-guided series of clinically-standardized movement tasks using both hands (Unified Parkinson’s Disease Rating Scale Part III: rest tremor, finger tapping, hand movements - open-close, pronation-supination) as well as a typing task (Figure 2A,B).

• Download figure
• Open in new tab

Figure 2. Optimizing the RC&S movement classifier.

(a) Six days of data were collected for model training, testing, and validation. During each session the participant performed a series of self-guided clinical tasks: rest recording (assessing tremor/dyskinesia), rhythmic finger tapping, hand open-close, wrist pronation-supination, and typing (pictured in order from left to right). The participant alternated between hands such that one hand was always resting while the other performed the tasks, except for when typing bimanually. (b) Wrist-worn accelerometry was recorded during behavior, and LFP’s were simultaneously recorded from three brain regions: the STN and the Postcentral and Precentral gyri. A three-minute example window is displayed showing contralateral wrist acceleration alongside beta (13-30Hz) and gamma (60-90Hz) power bands recorded from the PreCentral gyrus (PreC). A threshold was applied to the wrist accelerometry (red line) to label binary moving and non-moving states. Movement periods are indicated by orange shading, with movement tasks indicated by symbols corresponding to the cartoons in (a). (c) An automated pipeline for feature selection and parameter optimization was designed for programming the RC+S to optimally classify movement state. Prior to parameter optimization, two sets of candidate power bands were identified that would later be used as stimulation control signals: “Personalized” bands, computed using a modified form of principal components decomposition, and “Canonical” bands, defined by traditional frequency ranges used in existing literature. A progressive narrowing of each candidate pool began with creating separate models for every possible combination of power bands and feeding each through a brief (20 iteration) Bayesian optimization routine (“Broad search”). The top three models underwent an additional 200 iterations in a “Final optimization” stage, after which the single top performer was selected for each hemisphere.

Given the unique architecture of the neurostimulators’ algorithms, a custom machine learning pipeline was developed for optimization (Figure 2C). In brief, we first fit coefficient weights for predicting continuous-valued watch accelerometry from neural inputs. Using Bayesian optimization, we then tuned the parameters that define the algorithm dynamics and translate the accelerometry predictions into binary movement state predictions. Finally, we applied a feature selection routine to search for the optimal set of neural inputs to the algorithm.

Personalizing neural inputs improves algorithm performance

Since power band inputs may be defined over any frequency range, a method of constraining the search space was required. We used an approach that first identified candidate pools of potentially information-rich power bands, then progressively narrowed the combinatorial field to a single set of optimal inputs. We compared two methods of identifying candidates. The first method selected “canonical” power bands, representing commonly studied components of cortical and subcortical LFP’s related to movement (Buzsáki & Draguhn, 2004): theta (4-7Hz), alpha (8-12Hz), Beta (13-30Hz), and Gamma (60-90Hz). The second method identified “personalized” power bands, determined using a signal decomposition approach aimed at capturing the primary neural signals at each recording location.

The personalization approach first transformed time-domain LFP into a spectral representation using short-time fourier transform (STFT), then performed principal components analysis (PCA) on the z-scored STFT outputs to identify bands of neural activity capturing substantial spectral variance. A peak-finding algorithm was finally applied to the PCA component weights to find the frequency range best approximating each PC (Figure 3A). This process was run on an independent, multi-hour dataset. The resulting bands often overlapped canonically described power bands (Figure 3B). The top two PC’s in all cortical regions reflected personalized versions of beta and gamma, capturing 11-17% and 7-9% of the overall spectral variance, respectively. The top two PC’s in the STN were found in low-beta and gamma ranges, with gamma (11-31%) accounting for a larger proportion of the spectral variance than beta (5-7%). The third and fourth PC’s were more variable across brain areas, typically representing a combined theta-alpha band or low-gamma.

• Download figure
• Open in new tab

Figure 3. Personalizing power bands improves algorithm performance.

(a) To obtain Personalized power bands, principal component (PC) models were fit to spectrogram data from an independent dataset where the participant was freely behaving in their home for 2 hours. The PC weights (displayed for the top two PC’s of each recording channel) were used with a peak-finding algorithm to identify contiguous frequency ranges that approximated each PC in a way that was compatible with the device-embedded controller. Note that frequencies below 4Hz and above 100Hz were excluded due to known noise issues. (b) The top four personalized power bands were selected for each of the three recording locations, resulting in 12 unique biomarkers per hemisphere. The frequency ranges for each Personalized power band are displayed for comparison with Canonical power bands, with the color indicating the PC number (blue: PC0, orange: PC1, green: PC2, red: PC3; colors arbitrary for the canonical bands). (c) The marginal performance contribution for the top four power bands (ranked by mean, indicated by horizontal lines) are displayed for each power band pool. (d) The top 200 models from each of the biomarker pools (Canonical and Personalized) were directly compared to determine if the Personalized power bands provided significant performance benefits as compared to the Canonical power bands. This performance difference relates only to the choice of input power bands. Parameter optimization was performed identically for all models from both pools, therefore isolating the additive impact of power band personalization on top of the performance gains from the rest of the machine learning pipeline. Vertical dotted lines indicate the mean of each distribution. The top Personalized power band models outperformed their matched-rank Canonical models by an average of 6% (left hemisphere) and 2% (right). Asterisks indicate statistical significance at p<0.001 (permutation test - difference in means).

All combinations of power bands from the candidate pools were considered during feature selection. For each of the two candidate pools, combining all three brain regions, four power bands in each region, and at most four total inputs yields 682 models. These were narrowed down to three candidates per pool by selecting the models with the highest F1 scores after 20 iterations of Bayesian optimization (“Broad search”, Figure 2C). The models that advanced from the broad search were further optimized with an additional 200 iterations (“Final optimization”, Figure 2C), and the single top performing model was then selected for final analysis and aDBS experiments.

During optimization, we investigated whether inclusion of cortical recordings and power band personalization added value to the algorithm. We first computed the marginal performance contribution of each individual power band as a measure of feature importance using classifier accuracies from the broad search. The top power bands were uniformly located in cortical locations, indicating that ECoG signals were most valuable for decoding movement (Figure 3C). We then directly compared the two power band identification methods to assess the personalization procedure. Accuracy for the top personalized model was 6% greater than the top canonical model for the left hemisphere, and 3% greater for the right. For both hemispheres, this effect was consistent when comparing the top 200 models from each pool, displaying significantly greater mean accuracy for the personalized bands (p_left,p_right<0.001, permutation test; Figure 3D). Personalized power bands therefore offered significant performance improvements and were selected for subsequent analysis and adaptive therapy.

Evaluation of offline movement classification

We evaluated performance of the final neural classifiers using a single held-out session to determine whether they would be feasible for guiding adaptive stimulation (Figure 4A). Prediction accuracy for the left and right hemispheres was 83% and 76%, respectively (binomial test H₀: accuracy=0.5, p_left,p_right<0.001). Both True Positive Rate (TPR, 88% and 93%) and True Negative Rate (TNR, 81% and 68%) were significantly above chance (binomial test H₀: TPR, TNR=0.5, p_{left_TPR},p_{right_TPR,}p_{left_TNR},p_{right_TNR}<0.001; Figure 4B), yet there remained modest prediction bias towards the moving state (8% and 18%). Subthreshold acceleration of the inactive arm occasionally occurred during unimanual arm movements and was predicted by the classifier as contralateral movement (Figure 4A, red arrow), contributing to this bias. Accelerometry distributions segregated according to predicted movement state (Figure 4C), illustrating the classifier’s ability to discern true movements. Finally, moment-by-moment stimulation amplitudes were simulated (ramp rates: 1.6mA/sec; state targets: 1.6-2.2mA) to predict neurostimulator activity in later adaptive experiments (Figure 4D). As intended, stimulation was consistently high in each hemisphere during contralateral arm movements (left: 86% of the time, right: 90%), but not when stationary (left: 17%, right: 31%).

• Download figure
• Open in new tab

Figure 4. Final offline model performance.

(a) An example time window of decoder predictions is displayed with the true acceleration (top row), the predicted acceleration (middle row) and an overlay of the true (blue) and predicted (orange) movement states (bottom row). A notable false positive is indicated with red arrows on the right hemisphere plots, highlighting a common observation: even when the neural classifier incorrectly predicted movement, there was often observable sub-threshold acceleration. (b) Confusion matrices are displayed for decoder predictions from a held-out dataset. (c) Distributions of the accelerometry during predicted movement are compared to predicted stationary periods to further assess the discrimination performance of the neural decoder. Acceleration of the left hand (right hemisphere) notably lacked strong bimodality, making binary movement state predictions challenging. (d) The moment-by-moment stimulation current was simulated according to the predictions of the neural classifier and the stimulation parameters that were used in the later online adaptive experiments. The proportion of time spent at each stimulation amplitude is plotted for each true movement state of the offline test data.

Evaluation of movement classification during aDBS

We compared two aDBS algorithms with constant DBS for evaluating the therapeutic impact of movement-responsive stimulation (Figure 5A). The “Movement Responsive” condition consisted of increasing stimulation during predicted movement, while the “Inverted” condition had the inverse association, decreasing stimulation during predicted movement. The Inverted condition, predicted to have the opposite impact of the Movement Responsive condition, acted as a control for the possibility that stimulation intermittency was the driving factor of benefit. Both aDBS conditions varied stimulation between low (1.6mA) and high (2.2mA) levels, titrated prior to experimentation as tolerable for the participant. The third “Constant” condition provided conventional therapy at an intermediate level (1.9mA) regardless of classifier predictions. Data were collected for each of the three conditions in every session using a randomized and patient-blinded ordering. The same movement tasks that were completed during algorithm training were again performed in each block of the aDBS experiments, except for hand open-close, which was replaced with nose-tapping to represent ballistic arm movements involving proximal musculature. 12 sessions of data were collected over a 154 day period, with the last day performed 435 days after the original training data were collected.

• Download figure
• Open in new tab

Figure 5. Performance of neural classifiers is maintained during online adaptive stimulation.

(a) For evaluating movement-responsive DBS, three test conditions were used: “Movement responsive” increased stimulation during sensed movement, “Inverted” decreased stimulation during sensed movement, and “Constant” kept stimulation at an intermediate level. (b) Example data from the Movement responsive block of a single session displays the wearable accelerometry (top row) and neural predictions of the accelerometry (second row) together with the corresponding discrete movement states (third row) and the responsive stimulation adjustments (bottom row). (c) Confusion matrices depict the state-specific prediction accuracy from the neural classifier averaged across conditions. (d) The proportion of time at the low, high, and intermediate stimulation amplitudes was computed separately for each movement state and test condition. The goal of the algorithm was to maximize the time spent at high stim when in the moving state, while minimizing the time spent at high stim when in the stationary state (and vice-versa for the Inverted condition).

The movement classifier continued to perform well during aDBS (Figure 5B, Supplementary Video 1). Performance was stable for many months at a time; however, changes in neural signal magnitude were observed on two occasions more than a year after training data collection (once in each hemisphere, days 384 and 429). These were addressed through manual adjustment of algorithm thresholds without modifying other parameters, which restored performance to previous levels (Supplementary Figure 1). Prediction accuracy for the left and right hemispheres across all testing sessions was 82% and 76%, respectively (binomial test H₀: accuracy=0.5, p_left,p_right<0.001), and did not differ substantially across stimulation conditions (Supplementary Figure 2). The classifiers displayed mild prediction bias. An 8% bias towards predicting movement was observed for the left hemisphere, while a 14% bias towards no movement was observed for the right (Figure 5C). Overall, intended stimulation contingencies were achieved for each condition (Figure 5D).

Self-perceived therapeutic effects

We first obtained subjective therapeutic assessments by asking the participant to blindly score each of the blocks on perceived quality of movement (Figure 6A). One-way ANCOVA with performance order as a covariate revealed a significant effect of condition on self-scores (p=3.1e-4). Pairwise T-tests showed significantly reduced self-scores for the Inverted condition as compared to both the Movement Responsive (p=5.7e-4) and Constant (p=1.9e-4) conditions. Adjusting stimulation with respect to movement state therefore influenced subjective therapeutic quality, though primarily driven by reduction in quality reported during the Inverted condition.

• Download figure
• Open in new tab

Figure 6. Movement responsive stimulation impacts self-perceived therapeutic quality, movement speed, and dyskinesia in an accuracy-dependent manner.

(a) Blinded self-scores were reported by the subject at the conclusion of all 12 sessions to quantify the participant’s self-perceived therapeutic quality of each stimulation condition. The order in which the conditions were performed in each session was randomized in balanced fashion and was included as a covariate in all across-condition statistical testing. (b) The impact of classifier performance on therapeutic outcomes was assessed by regressing self-scores on the F1 score on a per-session basis. A single self-score was provided for both hands; the reported F1 score represents the average over the two hemispheric classifiers. Only data from the Movement responsive condition are shown. (c) As a control analysis to determine if therapeutic outcomes were simply better with higher average stimulation (e.g., from classifier biases), self-scores were regressed on mean stimulation amplitude across the two hemispheric stimulators on a per-session basis. (d-e) Quantifiable features of movement from representative motor tasks are shown independently for each hand, with the left bar and datapoints in each pair representing the left hand. (d) Movement speed was assessed using the repetition rate during two repetitive motor tasks: wrist rotations (left) and nose-tapping (right). Note that the right hemisphere (left hand) classifier displayed a bias towards predicting stationarity during movement. (e) Dyskinesia was assessed during the bimanual rest period where movement was assumed unintentional. Low-frequency power in the wrist-worn accelerometry data was used to quantify dyskinesia. Note that both hemispheric classifiers were highly accurate at predicting the stationary state. P-values were determined by one-way (a) and two-way (d-e) ANCOVA, and Pearson correlation (b-c). Asterisks indicate statistical significance at p<0.05. All box plots display mean, interquartile range, minimum, and maximum.

Furthermore, we sought to determine if self-scores depended upon algorithm performance or other confounding factors. A positive correlation between self-scores and classification performance is predicted under the foundational hypothesis that increasing stimulation during movement provides therapeutic benefit. We indeed observed a significant correlation between classifier performance and self-scores across sessions (Figure 6B; Pearson r=0.52, one-sided p=0.043). Additionally, we investigated whether greater average stimulation was simply driving perceived differences in therapeutic quality for the two aDBS conditions. Counter to this explanation, a significant negative correlation was observed between mean stimulation amplitude and self-score (Figure 6C; Pearson r=-0.44, two-sided p=0.032). Collectively, these results suggest that self-perceived motor improvements resulted from accurately targeting movement states and not indiscriminately increasing stimulation (e.g. from classification biases).

Movement speed and stimulation side-effects

We next tested the prediction that movement-locked increases in stimulation would enable faster movements, and stimulation decreases would slow movement. We quantified movement speed as the repetition rate for wrist pronation-supination and nose-tapping (Figure 6D). Two-way ANCOVA’s were performed with performance order as a covariate, revealing significant main effects of condition (p=1.3e-8), and hand (p=9.2e-3) for wrist rotations, and condition for nose-tapping (p=1.3e-3). A significant interaction between hand and condition was observed for both wrist rotations (p=1.1e-9) and nose-tapping (p=0.019). Pairwise T-tests showed that, for both movement types, speeds were slower in the Inverted condition compared to other stimulation conditions for the right hand only (p<0.01). This suggests that movement speed was modulated by stimulation condition for the right dominant hand alone. Importantly, classification of the moving state specifically for the left non-dominant hand was lower than any other state predictions (Figure 5C, Right Hemisphere). This resulted in stimulation being lower on average during left hand movements as compared to the right. However, stimulation was still consistently low during rest for both.

We further predicted that dyskinesia, a common side-effect of dopaminergic medication and DBS (Maciel et al., 2021), would be reduced during rest periods of the Movement Responsive condition when stimulation amplitudes were low. We quantified dyskinesia using low-frequency power (1-4 Hz) of the watch accelerometry (Rodríguez-Molinero et al., 2019; Figure 6E). Two-way ANCOVA with performance order as a covariate revealed significant main effects of condition (p=1.0e-6), and hand (p=0.018), but no significant interaction (p=0.61). Data for the two hands were therefore grouped for pairwise T-tests that showed significantly reduced dyskinesia in the Movement Responsive condition compared to other stimulation conditions (p_{vs_constant}=8.2e-6, p_{vs_inverted}=1.1e-6). We also evaluated tremor, quantified using accelerometry power in the 4-7 Hz range (Dai et al., 2015), which was previously identified as containing tremor markers for this participant at low stimulation (data not shown). We found no evidence that Movement Responsive aDBS worsened tremor, as power in the 4-7 Hz band was also reduced (Supplementary Figure 3). In summary, Movement Responsive stimulation provided a bilateral reduction in resting state dyskinesia, and a unilateral increase in movement speed impacting the dominant hand.

Naturalistic typing performance

We finally investigated movement quality in a bimanual setting by analyzing typing performance, predicting that movement-locked stimulation increases would enable faster action without compromising dexterity. Three metrics were analyzed (Figure 7): mean keypress duration (quickness of individual finger movements), overall typing speed, and backspace rate (error rate). For keypress duration, one-way ANCOVA with performance order as a covariate revealed a significant effect of condition (p=4.6e-7). Pairwise T-tests showed significantly shorter keypress durations for the Movement Responsive condition as compared to both other conditions (p_{vs_cconstant}=8.2e-6, p_{vs_inverted}=1.1e-6). For typing speed, a one-way ANCOVA again revealed a significant effect of condition (p=0.031). A significant pairwise difference was observed between the Movement Responsive and Inverted conditions (p=8.5e-3), with Movement Responsive typing being faster by 0.33 keypresses/s on average. No significant differences were observed for backspace rate (ANCOVA, p=0.092). Movement-responsive stimulation therefore impacted both quickness of individual finger movements and overall typing speed as predicted without significantly compromising dexterity.

• Download figure
• Open in new tab

Figure 7. Movement adaptive stimulation increases keypress responsiveness during naturalistic typing.

Three metrics of typing performance were evaluated across stimulation conditions. (a) Movement responsive stimulation resulted in a significant decrease in mean keypress duration relative to both of the other conditions. (b) Faster individual keypresses translated to an overall increase in typing speed when comparing only the Movement Responsive and Inverted conditions. (c) No significant differences in backspace rate (errors) were observed. KP: keypress; BS: backspace. Asterisks indicate statistical significance at p<0.05 (one-way ANCOVA). All box plots display mean, interquartile range, minimum, and maximum.

DBS implant and patient information

This study was reviewed by the University of California, San Francisco Institutional Review Board and registered on clinicaltrials.gov (IDE G180097; NCT03582891). One participant with PD was recruited from a parent trial testing slow aDBS and enrolled in the study after providing written informed consent. The participant had been previously implanted with one pair of electrode leads in each hemisphere, consisting of (A) quadripolar depth electrodes (Medtronic model 3389) in the STN, and (B) a quadripolar subdural ECoG array (Medtronic model 0913025) spanning the arm/hand regions of the precentral and postcentral gyri (motor and somatosensory areas). Each pair of electrode leads was connected to an investigational sensing-enabled INS (Summit RC+S, Medtronic) in the ipsilateral subclavicular space. The INS provided therapeutic stimulation to the STN through the depth electrodes, recorded LFP signals from the depth and ECoG electrodes, and recorded accelerometry data through an accelerometer in the INS case. All stimulation changes were self-initiated by the participant, and participant safety was ensured by allowing them to revert to their standard, pre-programmed DBS settings at any time using a handheld patient programmer device.

Neural data sensing configuration and signal processing

All neural recordings were performed during monopolar stimulation. Bipolar recordings in the STN were collected using a “sandwich” configuration (electrodes on either side of the stimulating contact) for common-mode rejection of the stimulation artifact. Two channels of neural data were recorded along the ECoG array, using the first and second electrodes for one channel (postcentral gyrus) and the third and fourth electrodes for the other (precentral gyrus). Neural time-domain data were recorded at a sampling rate of 500Hz. Device-embedded signal processing consisted of short-time Fourier transforms (STFT) with an FFT size of 256 pt and an interval of 50 ms.

Behavioral data streaming, signal processing, and synchronization

Detailed descriptions of the behavioral data collection platform can be found in Strandquist et al., 2023. Briefly, arm movement data were collected using Apple Watches (Apple, Inc) worn on each wrist and uploaded to a remotely accessible third-party data repository (Rune Labs, Inc.). These data were sampled at 50 Hz, and power was calculated using STFT with an FFT size of 64 pt and interval of 100 ms. Power was calculated in the 0-5 Hz range for movement decoding, 1-4 Hz for quantifying dyskinesia (Rodríguez-Molinero et al., 2019), and 4-7 Hz for quantifying tremor (Dai et al., 2015). Video data for participant observation and kinematic pose estimation were recorded on three cameras arranged at varied angles in the participant’s home office. A custom keylogging application was built to provide statistics regarding typing performance without compromising privacy by removing the letter identity of each keypress. This application was deployed on the participant’s home office computer with patient knowledge and consent. All data was securely transferred to remotely accessible databases over a virtual private network (VPN). To overcome drift and time discrepancies across system clocks and align the data streams, a synchronization event was performed at the start of all recordings. The participant tapped each side of their chest over the implanted INS while standing in view of the cameras. This movement created identifiable spikes in the INS and Apple Watch accelerometry signals and could be discerned on the video, allowing realignment of all data to a common time.

Data collection protocols

Three different datasets were collected remotely for this study. For identifying personalized power band candidates, a single two hour session of free behavior was collected during monopolar stimulation at 1.6 mA. The participant was instructed to go about their regular daily routine while occasionally taking breaks to perform hand open-close movements and rest their hands on their lap. Since this dataset was only used to identify the most prominent components of the neural signals through PCA, there were not firm constraints on the participant’s actions during data collection other than ensuring that there were periods of both movement and non-movement. Separate datasets were used for identifying the personalized power bands and evaluating their use in decoding movement to ensure that we avoided any risk of overfitting and represented generalizable performance.

A second dataset was collected for optimizing the movement classifier and performing feature selection from the identified power band candidates. This dataset included six days of data where the participant performed a prescribed set of motor tasks at three different levels of monopolar stimulation (1.2 mA, 2.4 mA, 2.6 mA). Each day consisted of approximately 25 minutes of total recording time. The participant was given an instruction sheet with the list of motor tasks, which included standardized movement tasks from the Unified Parkinson’s Disease Rating Scale Part III (UPDRS; sections 3.17-3.18 rest tremor, 3.4 finger tapping, 3.5 hand movements - open-close, 3.6 hand movements - pronation-supination) as well as blocks of text for typing that were each approximately 100 words long (Figure 2A,B). The participant was a proficient typer prior to beginning the study.

For evaluating the movement responsive aDBS algorithm, 12 days of data were collected while the participant performed a similar set of tasks in three different stimulation conditions (Movement-responsive, Inverted, Constant). Each day again consisted of approximately 25 minutes of total recording time. The open-close hand movements (UPDRS section 3.5) from the initial six day dataset were replaced with finger-to-nose tapping, and no movements of the lower body were performed. All movements were instructed to last 20 seconds, and the rest phase 10 seconds.

While optimal monopolar stimulation amplitudes were titrated to allow noise-free STN recordings with the sandwich sensing configuration, the participant’s standard therapeutic settings were conventional stimulation (cDBS) in a bipolar arrangement. In the monopolar setting, stimulation amplitudes were titrated separately for each of the three collected datasets before the first session. The high (2.2 mA) and low (1.6 mA) stimulation amplitudes for adaptive recordings were selected to be balanced around the constant stimulation amplitude (1.9 mA). These targets were also conservatively defined to ensure that the Inverted condition remained tolerable for the participant, since the inverted stimulation contingency was perceived as worse therapeutically. The main Movement Responsive test condition was therefore theoretically limited in how positive of an effect it could have, since its reciprocal needed to remain both counterbalanced and tolerable.

Personalized power band identification and assessment

Personalized power bands were identified by performing Principal Component Analysis (PCA) on the power spectral data. A custom peak-finding algorithm was then used on each PC weight vector to approximate it with a single contiguous frequency range that was compatible with RC+S signal processing capabilities (Figure 3A). The resulting power bands therefore (A) aligned with the dimensions of the neural data that captured the most variance, and (B) produced features that were approximately linearly uncorrelated to facilitate downstream linear combination for classification models. The steps are briefly outlined below:

1. Compute a [n_samples x L/2] spectrogram of the time-domain neural data using the rcssim signal processing with the following parameters:

   1. Fs_td = 500 Hz (time-domain sampling rate)

   2. L = 256 (FFT size)

   3. interval = 200 ms

2. Perform PCA on the z-scored spectrogram data and retain the top 4 PCs. Take the absolute values of each PC weight vector, which will be size [L/2 x 1] and represent the weight that a given PC places on each frequency bin in the previously performed FFT.

3. Run the peak-finding algorithm on the PC weight vectors to identify a contiguous range of frequency bins that best approximates each PC (restrict searchable range to 4-100Hz to avoid known noise issues):

   1. Excluding frequency bins that have been included in power bands that were identified prior, find the frequency bin with the largest weight for the current PC (i.e., the peak of the weight vector)

   2. Extend the current power band away from the peak in both directions until one of the following conditions are met:

      1. The frequency bin has been included in a previously identified power band (prevents power band overlap)

      2. The weight on the frequency bin drops below 50% of the peak weight (excludes frequencies that are not well-represented in the PC)

      3. A trough has been encountered (prevents multi-peak power bands)

      4. The frequency bin is more than 25Hz from the peak (limits power bands to 50Hz full range)

   3. Repeat for a total of 4 PC’s, independently for each anatomical recording area and brain hemisphere

Feature importance was assessed for each power band using its marginal performance contribution, which was computed as the performance increase (or decrease, as a result of overfitting) when a given power band was included as an additional feature in a model. Since models using every possible combination of four or less power bands were independently optimized in Step 1 of the feature selection process, the marginal performance contribution for a given power band was represented as a distribution of 196 values, one for each combination that included that power band.

Movement classifier and optimization pipeline

The Summit RC+S embedded aDBS has four primary components to its onboard operations. (1) Neural signals are first processed into power band outputs given a set of sensing parameters that include the frequency edges dictating each band. (2) Up to four power bands are streamed, averaged across a number of subsequent samples dictated by the update rate, and linearly combined via feature weights into a continuous valued linear discriminant output (LD). (3) LD outputs are then passed through a series of logic operations, which include comparison to a threshold value and requisite hold times above (onset duration) or below (termination duration) that value to trigger a state change. (4) Upon state change, stimulation is ramped towards the target amplitude associated with the new state. Note that this is a non-exhaustive description of the device capabilities that is focused on the most important features and parameters for the present study.

Because of the unique operations of the RC+S, a custom optimization pipeline was required for training a movement classifier that leverages its full functionality. In brief, each loop of GPBO simulated the embedded algorithm for a given set of parameters and neural inputs, evaluated the algorithm performance against labeled movement data, and modeled the relationship between parameters and performance to support an intelligent and efficient search-based optimization. LD weights were the only parameters not modeled and optimized through the GPBO search method, but instead were fit by linearly regressing the continuous-valued wearable accelerometry on the neural signals within each GPBO loop. By fitting these weights in closed form rather than as additional parameters in the GPBO gaussian process model, the GPBO routine was greatly simplified. Notably, the thresholds were not fit in closed form (e.g., through linear discriminant analysis), since it was assumed that interactions between thresholds and parameters such as onset duration would require that they be jointly optimized. Using the regression technique, LD outputs therefore represented an estimate of the continuous-valued accelerometry, which was then processed through the remaining algorithmic steps to produce binary state outputs given the set of parameters defined in the current GPBO loop. These state outputs were then compared against the true movement state determined from wearable accelerometry to calculate an F1-score. The F1-score was used as the objective function of the GPBO search, which was selected as the performance metric to promote unbiased and accurate predictions. Prediction bias reported elsewhere in the results was calculated as the difference between the frequency of labels predicted by the classifier and the frequency of labels in the actual data. Optimization was performed using five-fold cross-validation, using a separate day of data for each fold. All optimization code was written in Python, using the rcssim package for simulating RC+S operations and sklearn and skopt for model fitting and optimization.

As a minor stage in this process, a state change blanking variable dictating a refractory period to allow stabilization of artifacts associated with changes in stimulation was also optimized. This simple step consisted of modeling the artifact, which was different at cortical and subthalamic recording locations, and automatically assigning the appropriate duration of blanking given the neural inputs present in the model (i.e., whether inputs were from the cortex or subthalamic nucleus).

In summary, the following sensing parameters were optimized for each hemisphere:

1. Four power bands (though only using three was found to be optimal in the right hemisphere), each consisting of:

   1. Frequency lower limit

   2. Frequency upper limit

And the following algorithmic parameters were optimized for each hemisphere:

1. LD weights (four total, one for each power band)

2. Update rate

3. Threshold

4. Onset duration

5. Termination duration

6. Lag (not a device-programmable parameter, but allowed for time delay between neural signals and movement)

7. State change blanking

In total, 16 sensing parameters and 20 algorithmic parameters were therefore optimized.

Statistical testing

All statistical testing and data visualization was performed using the Python programming language and common statistical packages. Basic statistics were performed using scipy, and ANCOVA models were performed using statsmodels. Following ANCOVA, T-tests were performed as planned contrasts only in the event of statistically significant ANCOVA effects, thus no adjustments were made for multiple comparisons. The T-tests were covariate-adjusted for all statistically significant covariates identified through ANCOVA.