Continuous deep brain stimulation of the subthalamic nucleus may not modulate beta bursts in patients with Parkinson’s disease

Journal Pre-proof Continuous Deep Brain Stimulation of the Subthalamic Nucleus may not Modulate Beta Bursts in Patients with Parkinson’s Disease Stephen L. Schmidt, Jennifer J. Peters, Dennis A. Turner, Warren M. Grill PII:

S1935-861X(19)30479-6

DOI:

https://doi.org/10.1016/j.brs.2019.12.008

Reference:

BRS 1619

To appear in:

Brain Stimulation

Received Date: 20 May 2019 Revised Date:

19 November 2019

Accepted Date: 10 December 2019

Please cite this article as: Schmidt SL, Peters JJ, Turner DA, Grill WM, Continuous Deep Brain Stimulation of the Subthalamic Nucleus may not Modulate Beta Bursts in Patients with Parkinson’s Disease, Brain Stimulation, https://doi.org/10.1016/j.brs.2019.12.008. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 The Author(s). Published by Elsevier Inc.

Author Contributions:

Stephen L Schmidt: Conceptualization, investigation, software, formal analysis, writing – original draft, visualization. Jennifer J Peters: Investigation, software, writing – review and editing, visualization. Dennis A Turner: Conceptualization, writing – review and editing, funding acquisition. Warren M Grill: Conceptualization, writing – review and editing, supervision, funding acquisition.

BRS-D-19-00453R2 Continuous Deep Brain Stimulation of the Subthalamic Nucleus may not Modulate Beta Bursts in Patients with Parkinson’s Disease Abbreviated Title: Burst Dynamics in PD Stephen L. Schmidt1, Jennifer J. Peters1, Dennis A. Turner1,2, Warren M. Grill1,2,3 * 1

Biomedical Engineering Department, Duke University, Durham, NC, USA Neurobiology and Neurosurgery Departments, Duke University Medical Center, Durham, NC, USA 3 Electrical and Computer Engineering Department, Duke University, Durham, NC, USA *Indicates corresponding author. 2

Warren M. Grill Duke University Department of Biomedical Engineering Fitzpatrick CIEMAS, Room 1427 101 Science Drive Box 90281 Durham NC 27708-0281 [email protected] (919) 660-5276 Phone (919) 684-4488 Fax

1

Abstract Background: Neural oscillations represent synchronous neuronal activation and are ubiquitous throughout the brain. Oscillatory activity often includes brief high-amplitude bursts in addition to background oscillations, and burst activity may predict performance on working memory, motor, and comprehension tasks. Objective: We evaluated beta burst activity as a possible biomarker for motor symptoms in Parkinson’s disease (PD). The relationship between beta amplitude dynamics and motor symptoms is critical for adaptive DBS for treatment of PD. Methods: We applied threshold-based and support vector machine (SVM) analyses of burst parameters to a defined on/off oscillator and to intraoperative recordings of local field potentials from the subthalamic nucleus of 16 awake patients with PD. Results: Filtering and time-frequency analysis techniques critically influenced the accuracy of identifying burst activity. Threshold-based analysis lead to biased results in the presence of changes in long-term beta amplitude and accurate quantification of bursts with thresholds required unknowable a priori knowledge of the time in bursts. We therefore implemented an SVM analysis, and we did not observe changes in burst fraction, rate, or duration with the application of cDBS in the participant data, even though SVM analysis was able to correctly identify bursts of the defined on/off oscillator. Conclusion: Our results suggest that cDBS of the STN may not change beta burst activity. Additionally, threshold-based analysis can bias the fraction of time spent in bursts. Improved analysis strategies for continuous and adaptive DBS may achieve improved symptom control and reduce side-effects.

Keywords Deep Brain Stimulation, Beta Burst, Support Vector Machine, Subthalamic Nucleus, Parkinson’s Disease, Local Field Potential Analysis 2

Abbreviations aDBS

Adaptive deep brain stimulation

cDBS

Continuous deep brain stimulation

FFT

Fast Fourier transform

GP

Globus pallidus

IPG

Implantable pulse generator

LFP

Local field potential

MT

Mismatched thresholds (DBS OFF thresholds applied to DBS ON LFP)

PD

Parkinson’s disease

PSD

Power spectral density

SNR

Signal-to-noise ratio

SVM

Support vector machine

STN

Subthalamic nucleus

Conflict of Interest WMG is Co-Founder, Director and CSO of Deep Brain Innovations LLC. WMG is director and CSO of NDI Healthcare Fund and receives compensation for this position. These relationships are reported to the Conflict of Interest Committee at Duke University.

Acknowledgments The authors thank Brandon Swan, Isaac Cassar, David Brocker, and Chintan Oza for assistance collecting intraoperative data. We thank Gilda Mills and Danielle Degoski for laboratory support and Edgar Peña, Brandon Thio, and Nikki Pelot for their valuable input on the analysis. Funding Sources This work was funded by NIH UH3 NS103468 and NIH R37 NS040894.

3

Introduction Local field potentials (LFPs) reflect the underlying activity of populations of neurons. The presence of short, high amplitude “bursts” (i.e., short periods of high activity distinguishable from baseline) may contribute to dynamic signaling and to the average amplitude of the LFP [1, 2]. For example, bursts in the gamma band may predict better performance on working memory tasks, while bursts in the beta band may predict poor performance on motor and sensory tasks [3-6]. Dysregulated burst activity may also underlie disease pathology. For example, beta oscillatory activity is a candidate biomarker for the motor symptoms of Parkinson’s disease (PD) [7], and long-duration beta bursts appear to contribute to the correlation between the amplitude of beta oscillatory activity and symptoms [8, 9]. Beta amplitude in subcortical regions is correlated with motor symptoms [10, 11], and is reduced during treatment by either levodopa [12, 13] or continuous deep brain stimulation (cDBS) of the subthalamic nucleus (STN) [14, 15]. Since beta amplitude reliably changes with treatment it is also a candidate control signal for adaptive DBS (aDBS) [16-19]. Beta bursts are observed in motor cortex, striatum, STN, and globus pallidus (GP) in PD patients and healthy non-human primates [9, 20, 21]. There is consensus that longer bursts are pathological in PD while shorter bursts are related to normal motor function [5, 22, 23]. For example, beta bursts after a cue predicted longer reaction times [5] and increased time spent in bursts correlated with reduced velocity of movement [24]. Beta bursts are modulated by treatments that affect rigidity, tremor, and freezing of gait [8, 25, 26]. Thus, quantifying the amplitude, duration and rate of beta bursts could facilitate automatic adjustment of DBS parameters as an adjunct to averaged beta power over time [27]. However, quantifying such bursts is more challenging than it appears, and widely used thresholding techniques lead to biased results. A wide range of signal processing methods are currently employed to quantify and define bursts, but, depending on the method, resulting burst durations can vary from ~ 100 ms to 2-3 s for 4

patients with PD [22, 25, 28]. Herein, we demonstrate that many signal processing methods for quantifying bursts from LFP recordings, particularly the use of thresholds across varying amplitude signals, strongly bias the outcome of burst analysis leading to inaccurate estimation of bursts. We then implement support vector machine (SVM) classifiers, which generate substantially greater fidelity in quantifying bursts synthesized in an oscillator model. Subsequently, we apply the SVMs to LFPs from the STN of persons with PD and find no change in time in bursts, burst rate, or burst duration between cDBS OFF and cDBS ON. The results highlight the perils of threshold-based methods alone to analyze burst activity. The improved understanding of the analysis of dynamic changes in beta activity enabled by other methods is critical to understanding the role of beta in PD symptoms and may lead to better treatments, including improved aDBS controllers. Methods All data were analyzed by custom Matlab (The Mathworks, Natick, MA) scripts. The scripts for synthetic data generation and analysis are included in the research data attachment to this publication. Below are abbreviated methods and more detailed methods are provided in the Supplementary Material. Burst Oscillator Model Determining the effects of signal processing on burst activity required a signal with known characteristics, and a synthetic LFP was generated by combining variable numbers of cycles (i.e., bursts) of a 20 Hz sinusoid with noise and sampled at 1 kHz. Pink (1/f) noise was generated using the arbssnoise function [29] multiplied by 10 to increase the amplitude of the noise. Two parameters controlled the burst rate dynamics: the probability to enter a burst (rising edge) and the probability to leave a burst (falling edge). A vector of uniformly distributed random numbers between 0 and 1 was determined for each 20 Hz cycle of each 300 s duration signal, and we iterated through the random numbers to create bursts of different lengths. The signal-to-noise ratio (SNR) was calculated as the mean of the Hilbert

5

amplitude of the oscillator divided by the Hilbert amplitude of the noise signal. Exact parameters varied between analyses; see Supplementary Material for full details. Burst Identification and Analysis A number of different filtering approaches were implemented to approximate those used in publications that analyzed beta bursts from LFPs and were applied to the synthetic data. In addition, we added a short timescale (150 ms) FIR filter to capture better short burst dynamics. We compared short timescale FIR (150th order) and IIR (4th order) filters combined with the Hilbert transform to IIR, IIR filter combined with rectification and smoothing (200 ms), wavelet transform (7-cycles), and short-term FFT (500 ms). When using thresholds, bursts were identified as periods that exceeded a percentile (from 50th to 95th) of the distribution of signal amplitudes. Bursts were evaluated with minimum durations of 30, 100, and 200 ms. To compare our results to those with longer time scales, we also applied smoothing to the beta amplitude signal ranging from 0 to 1500 ms. As an alternative to threshold-based quantification of burst dynamics, four classification support vector machines (SVMs) were trained using the output of the FIR filter. The first pair of SVMs used the Hilbert transform to estimate the oscillator amplitude for each sample (1 ms). The second pair of SVMs estimated the amplitude for each zero-crossing (i.e., 25 ms for 20 Hz oscillations) of the filtered signal as the maximum or minimum value depending on the sign of the data. The second SVM of each pair was also provided as input the amplitude of the wideband signal calculated using the same method as the burst amplitude. The SVMs therefore find the best point on the beta amplitude axis or the best line on the beta amplitude versus wideband amplitude plane to divide the data (i.e., have one or two features), respectively. The output of each SVM was a label for each time point or half-cycle – one for in burst or zero for not in burst. When using SVMs, bursts were still restricted to periods longer than the minimum burst duration.

6

Burst duration was calculated as the mean time between the rising edge and falling edge of the threshold crossing (or in-burst labeling for SVM analyses) for all bursts in a trial. Similarly, the burst amplitude was the mean of all samples between the rising and falling edges of each burst in a trial. Burst rate was calculated by dividing the number of bursts by the total trial duration. The total amplitude in bursts was the product of the burst rate, mean burst duration, and mean burst amplitude. Human Participant Data Collection The Duke University Health System IRB approved the intraoperative clinical studies, and all participants provided written informed consent. Participants were patients undergoing either DBS lead implantation or implantable pulse generator replacement at Duke University Medical Center. The participants all were awake and under only local anesthesia at the time of data collection. Fourteen of sixteen participants (Table 1) refrained from taking their PD medications for at least 12 hours before surgery. A subset of the data analyzed here was analyzed using different methods and for a different purpose in a previous publication [30]. Trials were either 60 or 300 s in length and included either continuous 130 Hz or 185 Hz stimulation (cDBS ON, at an amplitude previously determined to be therapeutic) preceded by a trial without stimulation (cDBS OFF) of the same length, with a period of several seconds between trials. Biphasic charge-balanced pulses were delivered using a BP Isolator (FHC, Bowdoin, ME) in a monopolar configuration between either contact 1 or 2 of the DBS lead (as determined by the attending neurologist or the patient’s clinical settings) and a return electrode located on the chest. LFPs were recorded differentially between the contacts surrounding the stimulation contact, and amplified by 400 – 1000x, bandpass filtered from 0.1 Hz - 10 kHz (SR560 amplifiers, Stanford Research Systems, Sunnyvale, CA), and digitized at 50-100 kHz. Stimulation artifacts were rejected by diode clipping and input blanking during DBS pulses [31]. STN LFP Analysis

7

The frequency spectrum of the last 20 s of each trial was calculated using the pwelch() function. Beta (12-32 Hz) and broadband (1-55 Hz) power were calculated using the trapz() function. Broadband power was limited at 55 Hz to avoid 60 Hz line noise. The center of the beta band was determined as the maximum difference between the spectra during cDBS OFF and cDBS ON trials between 12 to 32 Hz. Since sampling rates of the LFP varied from 50 to 100 kHz, we first lowpass filtered the data at 400 Hz. We then down sampled the data to 1 kHz and applied a short-duration FIR 6-Hz bandpass filter centered at the most modulated beta band frequency with filtfilt(). The spectrum of the last 20 s of each trial was then recalculated, and the beta amplitude and wideband amplitude were calculated using the methods described for the burst model. The spectrum of each entire trial was then inspected for artifacts in the beta band used for burst analysis. Of the recordings from twenty-one participants, five participants were excluded due to noise in the beta band. Mismatched Threshold (MT) burst data were generated by applying thresholds at absolute amplitudes determined from cDBS OFF to the cDBS ON data. Due to the paired nature of the data and the lack of normality of the residuals, significance tests were performed using the Wilcoxon signed rank test and reported p values are not corrected for multiple comparisons. Synthetic OFF and ON paired data To test further our SVMs, we created randomized data with burst rates and durations matched to those reported in [23]. Our synthetic model was controlled by two parameters the probability of entering a burst and the probability of leaving a burst once entered. We created 16 OFF/ON pairs of mock LFPs to match the number of human participants. To produce variability between the mock participants, each pair of recordings had its own probability to enter and leave bursts scaled between 85% and 115% of the mean value. Results Signal Processing Methods Yield Varying Burst Characteristics from the Same Signals

8

A variety of amplitude estimation methods are reported to quantify beta amplitude and beta bursts, resulting in a wide range of calculated burst characteristics in patients with PD. These methods typically involve two steps: a step to restrict the signal in time (windowing) or frequency (filtering) and a step to estimate amplitude from the signal (such as the FFT). We first sought to determine if different analysis methods contributed to the wide range of reported burst characteristics. We compared filtering combined with the Hilbert transform, filtering combined with rectification and smoothing, wavelet transform, and windowed FFT to quantify the duration and rate of beta oscillatory bursts (Fig 1A). In addition, we added a short-timescale FIR filter, intended to capture better the short-term burst dynamics. We used the synthesized LFP, where the true characteristics were known, generated with a high SNR by combining an on/off burst oscillator with noise and set the average burst duration to be approximately 3 cycles of a 20 Hz oscillator (Fig. 1B) [32]. The different analysis methods resulted in substantially different beta amplitude time series and amplitude histograms (Fig 1B). For example, the rectified and smoothed IIR filter (purple) and FFT (green) did not separate the third and fourth bursts in the time series due to the short inter-burst interval. Further, the analysis methods yielded different estimates of burst rate, mean burst duration, mean burst amplitude and total amplitude in bursts (Fig. 1 C-F). For all values measured, the FIR filter followed by the Hilbert transform generated results closest to the true values (n = 100 trials, p < 10-3, Wilcoxon signed rank). The IIR followed by the Hilbert transform and wavelet transform also performed relatively well (difference between IIR and wavelet: p = 0.0144 for burst rate, p = 0.490 for burst duration, p < 10-3 for burst amplitude, and p < 10-3 for total amplitude in bursts). Conversely, IIR followed by rectification and smoothing, and windowed FFT did not provide reasonable estimates for the short bursts because their slow time constants (200 and 500 ms, respectively) precluded isolation of bursts with an average duration of less than 200 ms (p < 10-3 for all comparisons). These results clearly

9

demonstrate that different signal processing methods can produce very different burst characteristics from the same underlying data. Unknowable A Priori Knowledge of Time in Burst Is Necessary for Accurate Burst Analysis with Thresholds We applied a range of thresholds and smoothing durations to 100 iterations of the synthetic LFP from the burst oscillator subjected to FIR filtering across a range of burst rates. The specific threshold level and duration of additional smoothing (smoothing beyond the 300 ms Hamming window of our FIR filter) greatly influenced the number of bursts detected from the same underlying data (Fig. 2A). The true average burst duration was held constant so there was a strong correlation between the burst rate and the time spent in bursts (R = 0.996, Pearson’s correlation, Fig. 2B). For any given threshold (i.e., percentile of the distribution of signal amplitudes), the calculated burst rate did not match the true burst rate (Fig. 2C). If the estimated threshold percentile was closely matched to the actual time in bursts (e.g., 75th percentile for 25% of time in bursts) then the calculated burst rate approximated the true burst rate (0 ms smoothing: fitted line m = 0.859, goodness of fit R2 = 0.964, Fig. 2D). However, the proportion of time in bursts cannot be known a priori; therefore, we sought to determine whether accurate burst duration detection could be achieved by scaling the threshold. Reducing the threshold by only 10 percentiles (e.g., 65th percentile for 25% of time in bursts) resulted in poor estimates of burst duration (Fig. 2E). We also varied burst durations while reducing burst rate and again observed that threshold-based quantifications of burst characteristics varied widely for the same underlying data, dependent on the threshold that was selected (Supplementary Figure S1). These results indicate that accurate determination of burst parameters with threshold methods can only be achieved with an appropriate a priori estimate of the time in bursts, which, unfortunately, is unknowable. Threshold-Based Analysis of Beta Dynamics Recorded from Human STN

10

We recorded STN LFPs both with and without cDBS (n = 16 participants, Table 1). Beta amplitudes were estimated by bandpass filtering with a 6 Hz window centered around the beta frequency that was most modulated by cDBS (Fig. 3A-C). We were unable to determine the optimal amplitude thresholds from the amplitude histograms (i.e., there were no distinct modes in the histogram of beta amplitudes) (Fig. 3D, compare to Fig. 1B), and therefore we calculated burst metrics based on multiple threshold values and smoothing durations. Several studies claim that bursts longer than a certain duration are pathological while bursts that are shorter are not. We therefore determined the percent of all bursts that would be presumed pathological (> 500 ms, [8]) from the participant LFP and from white noise (Fig. 3E). The amount of smoothing applied greatly influenced the result, and, critically, threshold analysis of white noise bore striking similarities to the same analysis of participant LFP data. We observed that the product of the burst rate and mean duration were determined not by the underlying data, but rather by the threshold percentile (Supplementary Fig. S2). The biased results produced by threshold-based analyses raise the question of the impact of such analyses on true detection and interpretation of differences in beta oscillatory activity between conditions (e.g., cDBS OFF versus cDBS ON). Many analyses of beta bursts and aDBS methods assign a single threshold dependent upon the signal amplitude across the intervention. Thus, systematic changes in signal amplitude (e.g., reductions resulting from DBS) coupled with a fixed amplitude threshold will lead to different threshold percentiles. As demonstrated in Fig. 2, different threshold percentiles result in different burst characteristics even for the same input signals. Thus, changes in long-term beta amplitude between cDBS OFF and cDBS ON might produce artifactual differences in characteristics dependent on percentile thresholds.

11

We therefore examined the difference in beta band power between cDBS OFF and cDBS ON (Fig. 4A). Beta power significantly decreased with the application of cDBS (-4.14 [-6.58 – -1.67] % of total PSD from 1 – 55 Hz, median [95% CI], p = 0.0200 Wilcoxon signed rank test, n = 16 participants). There was no difference between the reductions in beta power with cDBS in the tremor- and bradykinesiadominant cohorts (p = 1.00 Wilcoxon ranked sum test, n = 9 tremor and 7 bradykinesia participants). To determine how beta power (over 10s of seconds) can affect thresholds, we examined the correlation between beta power and the square of the threshold amplitude, for amplitude thresholds set to either the 50th or 75th percentile of the signal amplitude distribution. The very strong correlations between the change in the squared threshold amplitude and the beta power indicated that the change in beta power before and after cDBS largely determined the altered threshold amplitudes (Fig. 4B, R = 0.965 and 0.962 for 50th and 75th percentile thresholds respectively, p < 10-3 Pearson correlation, n = 16 participants with one outlier excluded). Combined with the result that burst metrics were highly determined by the threshold value, this analysis indicates that measures of burst characteristics are strongly dependent on the overall (long-term) beta amplitude. We then examined the ratio of the threshold amplitudes between cDBS ON and cDBS OFF conditions. The mean ratio was relatively consistent across variations in the threshold and smoothing parameters, ranging from a ratio of 0.60 to 0.64 (Fig. 4C). This indicates that comparisons of burst characteristics between groups will be affected by the underlying decrease in the long-term beta amplitude on the time scale of seconds regardless of the amplitude threshold and smoothing duration that are chosen. One observation from prior beta burst analysis is an apparent increase in the proportion of short duration bursts (i.e., physiological) compared to long duration bursts (i.e., pathological) with either therapeutic DBS (both aDBS and cDBS) or medication. Our analyses suggest that this observation may be influenced by the thresholds for burst detection during DBS OFF and DBS ON being set at the same absolute amplitude, rather than the same percentile of the absolute signal amplitude (which differs 12

considerably between DBS OFF and DBS ON). Therefore, we quantified the effect of simply reducing the synthetic beta amplitude by a scalar factor intended to mimic the reduction in long-term beta during cDBS ON (Supplementary Fig. S3, Supplementary Table 2). We observed changes in burst duration and burst rate, in addition to anticipated changes in burst amplitude, as the beta amplitude was scaled, even though the underlying burst characteristics were unchanged. To determine whether this could bias results from human data, we examined in further detail the duration of bursts estimated from STN LFPs. We observed many fewer bursts when the threshold was the 95th percentile of cDBS OFF amplitude distribution as compared to the 50th, and the duration of bursts was shifted to favor short bursts (Fig. 5A). We then calculated the average burst duration for each threshold for cDBS OFF, cDBS ON and threshold amplitudes calculated by cDBS OFF applied to cDBS ON (MT, Fig. 5B). For lower thresholds we observed a decrease in burst duration with cDBS but observed a still greater decrease with MT. However, as the total number of bursts and MT trials with any detected bursts decreased, this difference was lost (for all thresholds see Supplementary Table 1). Similarly, when considering the fraction of total bursts in the shortest bin (30 – 100 ms), the fraction increased with cDBS but the proportion of short bursts was exaggerated by applying cDBS OFF thresholds to cDBS ON data (Fig. 5C). These data demonstrate that changes in beta burst characteristics between conditions may be exaggerated by application of single thresholds. While we observed a significant difference in the average burst duration, we also note that the average burst was shorter than the smoothing caused by the FIR filter. We also examined the same metrics with “DC correction” and averaged thresholds between conditions [8, 25] (Supplementary Figure S4). With DC correction, the beta amplitude is effectively normalized by mean subtraction. Because the range of the beta amplitudes is not also normalized, there remains a possibility of differences in measured burst activity between conditions. In our participant 13

cohort we observed significant differences in burst duration and fraction of bursts between 30 – 100 ms for thresholds from the 55th to 90th percentiles (Supplementary Table 2). Collectively, these results highlight that there may be significant differences in beta burst characteristics between conditions at a given threshold, however these differences can be exaggerated when the threshold is implemented by raw amplitude rather than by percentiles of the respective amplitude distributions. SVM Analysis of Burst Activity Threshold based analysis proved highly susceptible to systematic errors and may exaggerate differences in beta burst characteristics between conditions. We therefore implemented several SVMs to classify LFP as burst or non-burst. Generally, SVMs classify data by finding the hyperplane which best divides the labeled data. We implemented SVMs that used either continuous or cycle-by-cycle estimates of the beta band LFP amplitude. Half of the SVMs also included amplitude data from a wider band around the beta band, resulting in 4 different SVMs. Therefore, for each type of amplitude estimation, we implemented both a 1D (beta amplitude only) and 2D (beta amplitude and wideband amplitude) SVM. We used our burst oscillator model to generate labeled data to train the SVMs. The SVMs were trained on 50 s of labeled data with a burst rate of 1.08 Hz, an average duration of 346 ms, and a SNR of 4.18. We then compared the output of the SVMs to threshold analysis across 161 trials of synthetic LFP with randomly assigned burst rates, durations, amplitudes and noise levels (Fig. 6). The optimal threshold was applied at 62.7th percentile to match exactly the 50 s of training data. We observed that the accuracy of labeling was very high for all SVMs (> 95%, compared to 89% for threshold) and improved with increasing SNR. Unlike threshold analysis, the fraction in burst was free to vary and closely followed the underlying true data (fitted lines: m > 0.83, R2 > 0.89 for all SVMs). In contrast to threshold analysis, the mean burst duration appeared to be relatively independent of the time constant of the filter (Fig. 6C).

14

We then used each of the SVMs to label bursts in the intraoperative participant data (Fig. 7 shows representative data from the Hilbert amplitude only SVM, and results from all 4 SVMs are provided in Supplementary Table 3). We observed that the amplitude of bursts decreased with cDBS. Unlike threshold analysis, when using the SVMs the time in bursts could vary for the data from each participant (resulting in an effective threshold for each). We therefore examined the ratio of the effective threshold amplitude squared and the ratio of the power during cDBS ON to cDBS OFF. These ratios were significantly correlated, but to a substantially lesser extent than those for threshold analysis (r = 0.462, p = 0.0152). The duration of bursts was relatively short with no difference between cDBS OFF and cDBS ON (102 [90.9 - 111] ms and 94.2 [86.6 - 96.4] ms, cDBS OFF and cDBS ON respectively, p = 0.234 Wilcoxon signed rank). Burst rates were highly variable between participants but not different with cDBS (1.79 [1.51 - 1.84] Hz and 1.71 [1.38 - 1.87] Hz, cDBS OFF and cDBS ON respectively, p = 0.756 Wilcoxon signed rank). No differences in fraction in burst, burst duration nor burst rate were detected with any SVM (Supplementary Table 3). When considering all participants in the study we did not find significant improvement in intraoperative assessments of motor symptoms with cDBS at therapeutic frequencies (p = 0.313 for tremor power, p = 0.938 for bradykinesia scores, Wilcoxon signed rank, n = 9 for tremor and 7 for bradykinesia). It is likely that the microlesion effect occurring just after DBS electrode placement and lack of repeated trials contributed to this. We therefore examined bursts only in participants whose motor symptoms did improve with intraoperative cDBS (n = 9 participants, Supplementary Fig. 5). As with the full cohort, there was a change in burst amplitude driven by the decrease in long term beta amplitude, but no differences were detected in burst duration, burst rate, or fraction in burst (Supplementary Table 4). We also examined participants considered tremor dominant separately from those considered bradykinesia dominant (Supplementary Fig. 6). We observed a significant decrease in

15

burst amplitude for tremor cohort but not for bradykinesia dominant cohort, but no changes in burst duration (Supplementary Table 5). Differences between conditions are discoverable by SVMs We did not observe differences in burst activity that were independent of long-term beta amplitude with STN stimulation using SVM analysis. We considered two possibilities, 1) the SVM could not accurately determine differences between groups or 2) there are no differences in burst activity with cDBS in the STN. Using our burst oscillator model we generated 16 ON/OFF pairs of trials with mean burst durations and burst rates matched to those reported in [23] – changes in burst activity in the globus pallidus following administration of levodopa – and applied thresholds and SVM analyses to determine fraction in burst, burst durations and burst rate (Fig. 8). The fraction in bursts was not different between the ON and OFF conditions and observed by all analyses (Supplementary Table 6). All SVM analyses detected the difference in burst durations and burst rates, except for the cycle-by-cycle amplitude plus wideband SVM. This SVM found only a trend in the difference of burst rates. Conversely threshold-based analysis was unable to identify a difference in either mean burst duration or burst rate. Thus, SVMs can indeed isolate differences between groups, and the lack of differences detected in our participant data indicate that it is likely that STN cDBS did not modulate burst activity in this cohort. Discussion Changes in beta oscillatory activity correlate with motor symptoms and treatment of PD, and therefore provide a candidate biomarker for aDBS [7]. Many recent studies also considered burst dynamics as a candidate complementary biomarker. To our knowledge, this is the first study to address systematically the quantification of beta bursts. Typically, as beta bursts have been defined as intermittent periods of high amplitude, and some method is required to separate bursts from baseline activity. Threshold analysis is a common technique with separation based on some (often arbitrary) percentile or raw amplitude. However, we demonstrated that threshold analysis is unlikely to yield

16

accurate quantification of burst characteristics and leads to biased results in the face of changes in longterm beta amplitude. We therefore used an SVM-based approach which mathematically separates high amplitude bursts from background activity. Our results suggest that changes in beta burst parameters may be a redundant measure correlated to long-term beta amplitude and thus may not be an independent biomarker. Common signal processing methods for the extraction of bursts from LFP recordings bias the measurement of burst activity. Specifically, thresholding exaggerates the difference in burst activity between conditions that modulate the long-term signal amplitude, such as the application of cDBS, administration of levodopa, or with sustained movement. We applied a range of signal processing methods to quantify dynamic beta signal amplitudes and burst characteristics using both data synthesized from an oscillator model and STN LFPs recorded from participants with PD. Filtering, smoothing, amplitude estimation method, and selection of amplitude thresholds all significantly influenced the characteristics of detected bursts. Inappropriate selection of amplitude threshold and filter characteristics can lead to inaccurate results, and thresholds based on raw amplitudes, rather than normalized percentiles before and after PD treatment, exaggerated differences between conditions. The reduction of the average (long-term) beta amplitude by DBS primarily determines the differences in burst amplitude, duration, and rate between the ON and OFF conditions. Subsequently, we trained SVMs on model data and observed that the SVM analyses did not have the same limitations as threshold analysis. When applied to human data, we observed no differences in burst dynamics with the application of cDBS, even though our SVM analyses were indeed able to differentiate between groups using parameters reported from a cohort of PD patients ON and OFF levodopa [23]. The Effects of Filters and Thresholds

17

Filtering may profoundly change the result of data analysis [33], and at a minimum, filtering imposes some degree of smoothing [34]. In conjunction with the principle of matched filtering [35], bursts with durations approximately equal to the duration of the filter are preferentially identified, and thus our results are related to the duration of the FIR filter. Depending on the smoothing applied in the estimation of beta amplitude, the resulting burst durations may arbitrarily vary (Fig. 1, Fig. 4E). However, SVM analyses did not preferentially select any particular duration of bursts (Fig. 6C). It remains unclear what the minimum duration of a burst should be. Our short latency FIR filter was tuned to detect short-term bursts on the scale of ~100 ms because they closely align with spiking data [9, 32, 36, 37]. The minimum burst duration was set to slightly shorter than a single cycle of the fastest frequency (30 ms) to account for the smoothing of the filter. We observed differences in raw values of burst rate and duration with minimum durations of 100 and 200 ms. However, the effects of threshold percentile and additional smoothing remain at these minimum durations (Supplementary Figure S2). It is possible that repeated short bursts are mis-identified as a single sustained burst when smoothed over a longer duration and are relevant in PD. Indeed, prototype aDBS in PD, where contingent or modified stimulation is applied based on beta amplitude, was effective in a number of studies [17-19, 38-41]. However, if clinically relevant bursts are dependent on filtering, it is unclear what timescale and analysis of beta dynamics may lead to the best results for patients with PD, or if the timescale is a parameter that should vary between patients. Although threshold amplitude has a critical role in the analysis of burst activity, this method remains a main determinant of aDBS contingent stimulation and determines the amount of time that aDBS is ON. Beta amplitude dynamics and burst activity during the transition between the aDBS OFF and aDBS ON are of particular interest for aDBS [24]. Indeed, assessment of threshold after application of DBS (with a resulting reduction in long-term beta amplitude) was used to preset the amount of time aDBS is active [16]. Similarly, the amplitude of the stimulation may be scaled with the beta amplitude rather than switching the stimulation on or off, 18

and studies of aDBS implementing proportional control of stimulation show promising results [17, 19, 41, 42]. Beta Bursts and Long-Term Beta Power Because bursts are both transient and relatively high amplitude it is tempting to think that longterm beta power is determined by the relative presence of bursts. The presence of high-amplitude activity will increase the average power, but it does not necessarily follow that bursts are the primary driver of total power. Conceptually, beta power may be divided into burst and background power. When using a threshold, it is not possible to determine whether treatment by DBS or levodopa decreased the amplitude of bursts or only the amount of time spent in bursts. SVM analyses did not reveal any difference in time spent in bursts between cDBS OFF and cDBS ON (Fig. 7B). The decrease in burst amplitude (Fig. 7A) shows that burst amplitude is not agnostic to treatment. However, the correlation in Fig. 7C, while significant, does not explain most of the variance. Together these data demonstrate that both background and burst amplitudes were decreased by cDBS within the cohort tested here. It should be noted that while burst amplitude changes with stimulation it is not clear that this is a complementary biomarker for disease state that is independent of long-term beta power. Analysis of Short-Term Amplitude Dynamics SVMs are not the only possible analysis method that should be less susceptible to error than amplitude threshold. Indeed, Markov models and other advanced burst analysis methods are beginning to appear in neuroscience [3, 43, 44]. However statistical and machine-learning methods require training data which, in the case of DBS, may require long-term recording capability. Here we generated model data to train the SVM. Action potentials are more likely to be synchronized during beta and hippocampal theta bursts [32, 36, 45]. Therefore, training may be performed with high confidence on

19

LFPs when action potentials are simultaneously recorded. The algorithms may then be adapted to macroelectrode LFP recordings. Spectral analysis in participants with PD is usually referenced to the power or amplitude in a broad band (here 1-55 Hz). However, burst analysis typically depends on isolation of a single frequency band. In the case of burst analysis in subjects with PD, the beta band is isolated from the other bands of the spectrum by filtering. Inter-patient variability is then normalized by application of thresholds based on percentiles of the amplitude distribution. However, this approach rescales the amplitudes without consideration of the underlying system. For example, beta amplitudes may be high because the participant has potentiated beta oscillations that are correlated to pronounced motor symptoms, or because the sensing electrodes were placed in an optimal environment for sensing. Here we tried to capture this with the two SVMs additionally trained with wideband data. However, the gains in accuracy with this additional input were limited. Limitations Intraoperative and immediate postoperative data collection offer unique, short-term opportunities for research but also several disadvantages [46]. Intraoperative data collection is limited to short periods with brief (usually less than 10 s) pauses between OFF and ON and the electrode conditions may not necessarily be stable (although more so with patients undergoing IPG replacements years after DBS lead implantation). In addition, while participants were awake and had withdrawn from PD-related medications before surgery, some changes to LFP may outlast the sedation. Due to the limited time, we were only able to perform one method of behavioral assessment and opted for tremor when it was present to be able to conduct a greater number of trials. Therefore, the classification of ‘tremor dominant’ does not preclude the presence of bradykinesia. We did not observe a correlation between long term beta power and improvement in accelerometer or tapping scores for these cohorts

20

using methods described in [30]. Considering the limitations of our intraoperative data we did not investigate a correlation between burst activity and change in motor symptoms. It is possible that a larger number of participants, multiple trials per participant and behavioral assessments on roughly the same timescale as the bursts will prove useful to find optimal analysis settings. Given the inaccuracy in burst rate and duration calculations by suboptimal thresholds (Fig. 2E), we suspect that any optimal threshold analysis will be determined on an individual level, perhaps changing in short intervals, as ongoing beta activity is continuously re-analyzed. The SVMs that take only the beta amplitude estimation as input select a static threshold and apply it to the entire recording. We assumed that beta dynamics were different between DBS ON and DBS OFF but relatively similar within during a single recording (60 or 300 s of DBS ON or DBS OFF). We showed that changes in the long-term beta amplitude alter the resulting burst dynamics when a single amplitude threshold is applied. However, in certain cases it may not be possible to separate periods of high and low amplitude beta oscillations. For example, the amount of beta in the STN is influenced by contralateral stepping [47]. In such cases, it will be necessary to train SVMs or other machine learning algorithms with appropriate controlling variables (i.e. heel position). Burst identification with SVMs may be improved by adding previous estimates of the amplitude as independent features rather than the single estimate used here. We avoided such an approach because it would apply smoothing, in addition to that of the filter, to the output (as a weighted average over the time of the features). It should also be noted that for aDBS, the SVM using multiple previous amplitude estimates as features would introduce delay in burst identification. The number of past amplitude features should therefore be chosen such that the controller can respond in time to reduce bursts effectively. Our computational model was limited by the simplicity of the generated LFP. We attempted to mitigate the effects on the SVM by adding noise, randomizing trial parameters, and randomizing individual burst parameters. However, no oscillator model can capture the complexity of in vivo LFP. It is 21

unclear in what way this simplicity biases the SVM when applied to participant data. Future studies could utilize training based on LFP labeled by simultaneous spiking or training by unsupervised machine learning algorithms to mitigate this limitation. Figure Captions Figure 1. Comparison of signal processing methods to analyze bursts present in local field potentials. A) Diagram of several commonly used methods to estimate instantaneous amplitude from the LFP. The 2nd row (FIR Filter, etc.) corresponds to methods to isolate data in time or frequency while the 3rd row (Hilbert, etc.) corresponds to amplitude estimation methods. Wavelet convolution (yellow) can be considered in both rows, as it both isolates in time and frequency as well as estimates the amplitude. B) Example amplitude estimation. Top: a five second trace of the synthetic signal generated from the on/off oscillator (black) with noise added (red) on a normalized amplitude scale. Middle: normalized histograms of signal amplitude for each analysis method with the 85th percentile indicated (red lines). Bottom: The resulting amplitude time series for each of the five methods. C-F). The burst rate, mean burst duration, mean burst amplitude and total amplitude in bursts estimates for each method and the true values. Error bars represent the 5-95% confidence interval of the median across 100 replications of 300 s epochs of the synthetic LFP. Figure 2. Threshold analysis. A) Top: a 5 second trace of the output of FIR filter (blue) with the thresholds for 50th to 95th percentiles of the signal amplitude indicated by the horizontal lines. The number and duration of bursts are highly dependent on the percentile chosen. Center: the underlying oscillator data for the same period. Bottom: amplitudes and 69.3rd percentile thresholds for 4 different levels of smoothing. The number and duration of bursts are highly dependent on the chosen smoothing duration. B) Scatter plot of the burst rate and total fraction of time in bursts across trials, and the best linear fit. Because the average true burst duration was relatively constant, the burst rate dictated the

22

change in the total time spent in bursts (R = 0.996, Pearson correlation). C) Calculated vs true burst rates for 50th, 75th and 95th threshold percentiles for all trials. Accurate estimations of the burst rate lie on the unity line (black). D) Calculated vs true burst rate when the selected threshold was determined by the fraction of time spent in bursts. Accuracy was much improved compared to (C). E) Calculated vs true burst rate when the selected threshold was 10 percentiles less than that determined by the fraction of time in bursts. Accuracy is low compared to the matched data in (D). Add. Smooth.: additional smoothing beyond that applied by the FIR filter. Perc: Percentile. Fig. 3 Analysis of beta oscillatory activity in the STN of participants with Parkinson’s disease. A) Power spectra of the final 20 s epoch of cDBS OFF (blue) and cDBS ON (orange) for one participant. B) The difference in spectral density caused by the application of cDBS. Beta was greatly reduced at approximately 13 Hz. C). Spectra of the final 20 s of the cDBS OFF trial before (blue) and after (green) the application of the FIR filter centered at the most attenuated beta frequency. D). Histogram of the beta oscillation amplitudes with 50th, 75th and 95th percentiles indicated. Compare to the histograms of the oscillator model in Fig 1B. E) Percentage of bursts deemed pathological (> 500 ms duration) for cDBS OFF (left) and white noise (right) with a minimum burst duration of 100 ms. Perc: percentile. PSD: power spectral density. Fig. 4. Effects of cDBS on burst analysis. A) The percent of the broadband power spectral density contained in the beta band (12-32 Hz) for cDBS OFF and cDBS ON for tremor dominant (blue circles) and bradykinesia dominant (orange diamonds) participants. Application of cDBS generally reduced the percent of the total power within the beta band. There was not a difference in the decreases in beta power between tremor- and bradykinesia-dominant cohorts. B). Scatter plots of the ratio between cDBS ON to cDBS OFF of beta power and ratio of the square of the threshold amplitude set to either 50th or 75th percentile of the signal amplitude distribution. One participant not plotted as an outlier due to a change in beta power of ~8-fold. Red lines are best linear fits of the data. The correlation between the 23

change in beta power spectral density and squared threshold amplitude was very high (r > 0.95 for all thresholds) indicating that thresholds are largely a measure of long-term beta power. C) The change in amplitude threshold value for all combinations of threshold percentile and applied smoothing. Fig. 5. Amplitude threshold analysis of STN LFP. A) The total number of bursts within 100 ms duration bins for 50th percentile thresholds (left) and 95th percentile threshold (right). CDBS ON, cDBS OFF, and MT are included for each bin. Application of thresholds calculated during cDBS OFF to cDBS ON data (MT, yellow) resulted in fewer bursts generally. Increased threshold percentiles resulted in fewer bursts. B) The mean burst duration across conditions for 50th and 95th percentile threshold. At the 50th percentile, bursts were significantly shorter for cDBS ON than cDBS OFF, however MT was significantly shorter still. This result indicates that single thresholds may exaggerate the difference between conditions. No significant difference was found at the 95th percentile threshold, likely due to the decreased total number of bursts. C) Similar results were observed for the fraction of all bursts between 30 – 100 ms duration. MT: mismatched threshold. Fig. 6. SVM analysis of model data. A) Scatter plot of in and out of burst label accuracy versus SNR for each trial. SVMs were highly accurate with better performance as the SNR increased. B-E) Calculated versus true fraction in burst, burst duration, burst rate and burst amplitude. SVM analyses were free to vary in each of these metrics resulting in calculations nearly matching the true values. Threshold analysis was limited to a single fraction in burst and therefore burst duration and rate were constrained. Cont. Amp: SVMs trained on Hilbert amplitude of the input signals. Cycle: SVMs trained on the maximum amplitude of each zero-crossing of the input signals. WB: SVMs trained on both the oscillation and a wideband amplitude input. Fig. 7. SVM analysis of participant LFP. A) Burst amplitude during cDBS OFF and cDBS ON. B) Fraction in burst during cDBS OFF and cDBS ON. C) Scatter plot of the ratio of squared effective threshold amplitude

24

(the amplitude corresponding to a threshold matching the fraction in bursts) versus the ratio of beta power between cDBS ON and cDBS OFF. The red line represents the best linear fit of the data. D). The median burst duration was ~ 100 ms for both cDBS OFF and cDBS ON. E) No significant difference in burst rates between cDBS OFF and cDBS ON. N.S: not significant (p > 0.05, Wilcoxon ranked sum). Data shown were calculated with the Hilbert amplitude SVM but are representative of all SVMs tested. Fig. 8. SVM analyses can distinguish between condition groups. A) Fraction of time spent in bursts for each condition (blue or orange dots, OFF or ON respectively) and analysis group for n = 16 trial pairs. There was no significant difference between the underlying (True) data which was detected by all analyses. B) Many analyses were able to determine differences in burst durations between OFF and ON, but not thresholding and Cycle with Wideband. C) All analyses except thresholding were able to distinguish significant differences between burst rates of OFF and ON in the underlying data. *: p < 0.05 (Wilcoxon signed rank, no correction for multiple comparisons). Thresh: threshold analysis at the 62.7th percentile (the optimal threshold for the training data). Amp: SVMs trained on Hilbert amplitude of the input signals. Cycle: SVMs trained on the maximum amplitude of each zero-crossing of the input signals. WB: SVMs trained on both the oscillation and a wideband amplitude input. Table 1. Participant Characteristics (N = 16).

25

References [1] Jones SR. When brain rhythms aren't ‘rhythmic’: implication for their mechanisms and meaning. Current Opinion in Neurobiology 2016;40:72-80. [2] van Ede F, Quinn AJ, Woolrich MW, Nobre AC. Neural Oscillations: Sustained Rhythms or Transient Burst-Events? Trends Neurosci 2018;41(7):415-7. [3] Becker R, Vidaurre D, Quinn AJ, Abeysuriya R, Parker Jones O, Jbabdi S, et al. Transient spectral events in resting state MEG predict individual time-frequency task responses. bioRxiv 2018:419374. [4] Lundqvist M, Rose J, Herman P, Brincat SL, Buschman TJ, Miller EK. Gamma and Beta Bursts Underlie Working Memory. Neuron 2016;90(1):152-64. [5] Torrecillos F, Tinkhauser G, Fischer P, Green AL, Aziz TZ, Foltynie T, et al. Modulation of Beta Bursts in the Subthalamic Nucleus Predicts Motor Performance. The Journal of neuroscience : the official journal of the Society for Neuroscience 2018;38(41):8905-17. [6] Shin H, Law R, Tsutsui S, Moore CI, Jones SR. The rate of transient beta frequency events predicts behavior across tasks and species. Elife 2017;6. [7] Hoang KB, Cassar IR, Grill WM, Turner DA. Biomarkers and Stimulation Algorithms for Adaptive Brain Stimulation. Front Neurosci 2017;11:564. [8] Tinkhauser G, Pogosyan A, Little S, Beudel M, Herz DM, Tan H, et al. The modulatory effect of adaptive deep brain stimulation on beta bursts in Parkinson's disease. Brain : a journal of neurology 2017;140(4):1053-67. [9] Feingold J, Gibson DJ, DePasquale B, Graybiel AM. Bursts of beta oscillation differentiate postperformance activity in the striatum and motor cortex of monkeys performing movement tasks. Proc Natl Acad Sci U S A 2015;112(44):13687-92. [10] Kuhn AA, Kupsch A, Schneider GH, Brown P. Reduction in subthalamic 8-35 Hz oscillatory activity correlates with clinical improvement in Parkinson's disease. Eur J Neurosci 2006;23(7):1956-60. [11] Pogosyan A, Yoshida F, Chen CC, Martinez-Torres I, Foltynie T, Limousin P, et al. Parkinsonian impairment correlates with spatially extensive subthalamic oscillatory synchronization. Neuroscience 2010;171(1):245-57. [12] Marceglia S, Foffani G, Bianchi AM, Baselli G, Tamma F, Egidi M, et al. Dopamine-dependent non-linear correlation between subthalamic rhythms in Parkinson's disease. J Physiol 2006;571(Pt 3):579-91. [13] Priori A, Foffani G, Pesenti A, Tamma F, Bianchi AM, Pellegrini M, et al. Rhythm-specific pharmacological modulation of subthalamic activity in Parkinson's disease. Experimental neurology 2004;189(2):369-79. [14] Eusebio A, Cagnan H, Brown P. Does suppression of oscillatory synchronisation mediate some of the therapeutic effects of DBS in patients with Parkinson's disease? Front Integr Neurosci 2012;6:47. [15] Whitmer D, de Solages C, Hill B, Yu H, Henderson JM, Bronte-Stewart H. High frequency deep brain stimulation attenuates subthalamic and cortical rhythms in Parkinson's disease. Frontiers in human neuroscience 2012;6:155. [16] Little S, Pogosyan A, Kuhn AA, Brown P. beta band stability over time correlates with Parkinsonian rigidity and bradykinesia. Experimental neurology 2012;236(2):383-8. [17] Arlotti M, Marceglia S, Foffani G, Volkmann J, Lozano AM, Moro E, et al. Eight-hours adaptive deep brain stimulation in patients with Parkinson disease. Neurology 2018;90(11):e971-e6. [18] Priori A, Foffani G, Rossi L, Marceglia S. Adaptive deep brain stimulation (aDBS) controlled by local field potential oscillations. Experimental neurology 2013;245:77-86. [19] Rosa M, Arlotti M, Ardolino G, Cogiamanian F, Marceglia S, Di Fonzo A, et al. Adaptive deep brain stimulation in a freely moving Parkinsonian patient. Movement disorders : official journal of the Movement Disorder Society 2015;30(7):1003-5. 26

[20] Pina-Fuentes D, van Zijl JC, van Dijk JMC, Little S, Tinkhauser G, Oterdoom DLM, et al. The characteristics of pallidal low-frequency and beta bursts could help implementing adaptive brain stimulation in the parkinsonian and dystonic internal globus pallidus. Neurobiology of disease 2019;121:47-57. [21] Tinkhauser G, Torrecillos F, Duclos Y, Tan H, Pogosyan A, Fischer P, et al. Beta burst coupling across the motor circuit in Parkinson's disease. Neurobiology of disease 2018;117:217-25. [22] Deffains M, Iskhakova L, Katabi S, Israel Z, Bergman H. Longer beta oscillatory episodes reliably identify pathological subthalamic activity in Parkinsonism. Movement disorders : official journal of the Movement Disorder Society 2018;33(10):1609-18. [23] Lofredi R, Neumann WJ, Brucke C, Huebl J, Krauss JK, Schneider GH, et al. Pallidal beta bursts in Parkinson's disease and dystonia. Movement disorders : official journal of the Movement Disorder Society 2019;34(3):420-4. [24] Lofredi R, Tan H, Neumann WJ, Yeh CH, Schneider GH, Kuhn AA, et al. Beta bursts during continuous movements accompany the velocity decrement in Parkinson's disease patients. Neurobiology of disease 2019. [25] Tinkhauser G, Pogosyan A, Tan H, Herz DM, Kuhn AA, Brown P. Beta burst dynamics in Parkinson's disease OFF and ON dopaminergic medication. Brain : a journal of neurology 2017;140(11):2968-81. [26] Anidi C, O'Day JJ, Anderson RW, Afzal MF, Syrkin-Nikolau J, Velisar A, et al. Neuromodulation targets pathological not physiological beta bursts during gait in Parkinson's disease. Neurobiology of disease 2018;120:107-17. [27] Moraud EM, Tinkhauser G, Agrawal M, Brown P, Bogacz R. Predicting beta bursts from local field potentials to improve closed-loop DBS paradigms in Parkinson's patients. Conf Proc IEEE Eng Med Biol Soc 2018;2018:3766-96. [28] Wang DD, de Hemptinne C, Miocinovic S, Ostrem JL, Galifianakis NB, San Luciano M, et al. Pallidal Deep-Brain Stimulation Disrupts Pallidal Beta Oscillations and Coherence with Primary Motor Cortex in Parkinson's Disease. The Journal of neuroscience : the official journal of the Society for Neuroscience 2018;38(19):4556-68. [29] Zhivomirov H. A Method for Colored Noise Generation. Rom J Acoust Vib 2018;15(1):14-9. [30] Brocker DT, Swan BD, So RQ, Turner DA, Gross RE, Grill WM. Optimized temporal pattern of brain stimulation designed by computational evolution. Science translational medicine 2017;9(371). [31] Kent AR, Grill WM. Recording evoked potentials during deep brain stimulation: development and validation of instrumentation to suppress the stimulus artefact. Journal of neural engineering 2012;9(3):036004. [32] Murthy VN, Fetz EE. Coherent 25- to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys. Proc Natl Acad Sci U S A 1992;89(12):5670-4. [33] de Cheveigné A, Nelken I. Filters: When, Why, and How (Not) to Use Them. Neuron 2019;102(2):280-93. [34] Gabor D. Theory of communication. Part 1: The analysis of information. Journal of the Institution of Electrical Engineers - Part III: Radio and Communication Engineering. 93. 1946:429-41. [35] North DO. An Analysis of the factors which determine signal/noise discrimination in pulsedcarrier systems. Proceedings of the IEEE 1963;51(7):1016-27. [36] Murthy VN, Fetz EE. Synchronization of neurons during local field potential oscillations in sensorimotor cortex of awake monkeys. Journal of neurophysiology 1996;76(6):3968-82. [37] Cagnan H, Mallet N, Moll CKE, Gulberti A, Holt AB, Westphal M, et al. Temporal evolution of beta bursts in the parkinsonian cortical and basal ganglia network. Proc Natl Acad Sci U S A 2019;116(32):16095-104.

27

[38] Swann NC, de Hemptinne C, Thompson MC, Miocinovic S, Miller AM, Gilron R, et al. Adaptive deep brain stimulation for Parkinson's disease using motor cortex sensing. Journal of neural engineering 2018;15(4):046006. [39] Little S, Beudel M, Zrinzo L, Foltynie T, Limousin P, Hariz M, et al. Bilateral adaptive deep brain stimulation is effective in Parkinson's disease. Journal of neurology, neurosurgery, and psychiatry 2016;87(7):717-21. [40] Little S, Pogosyan A, Neal S, Zavala B, Zrinzo L, Hariz M, et al. Adaptive deep brain stimulation in advanced Parkinson disease. Ann Neurol 2013;74(3):449-57. [41] Velisar A, Syrkin-Nikolau J, Blumenfeld Z, Trager MH, Afzal MF, Prabhakar V, et al. Dual threshold neural closed loop deep brain stimulation in Parkinson disease patients. Brain stimulation 2019. [42] Rosa M, Arlotti M, Marceglia S, Cogiamanian F, Ardolino G, Fonzo AD, et al. Adaptive deep brain stimulation controls levodopa-induced side effects in Parkinsonian patients. Movement disorders : official journal of the Movement Disorder Society 2017;32(4):628-9. [43] Hirschmann J, Baillet S, Schnitzler A, Woolrich M, Vidaurre D, Florin E. Spontaneous network activity accounts for variability in stimulus-induced gamma responses. bioRxiv 2018:381236. [44] Cole S, Voytek B. Cycle-by-cycle analysis of neural oscillations. Journal of neurophysiology 2019. [45] Cole SR, Voytek B. Hippocampal theta bursting and waveform shape reflect CA1 spiking patterns. bioRxiv 2018:452987. [46] Swan BD, Grill WM, Turner DA. Investigation of deep brain stimulation mechanisms during implantable pulse generator replacement surgery. Neuromodulation : journal of the International Neuromodulation Society 2014;17(5):419-24; discussion 24. [47] Fischer P, Pogosyan A, Green AL, Aziz TZ, Hyam J, Foltynie T, et al. Beta synchrony in the corticobasal ganglia network during regulation of force control on and off dopamine. Neurobiology of disease 2019;127:253-63.

28

Sex male 14 (87.5) female 2 (12.5) Race white 18 (100) Ethnicity Hispanic 1 (6.25) non-Hispanic 15 (93.75) Age mean years +/- S.D. 63.7 +/- 7.4 range 51 - 74 Motor Symptom Evaluated tremor 9 (56.25) bradykinesia 7 (43.75) Surgery Type Implant 10 (62.5) Implantable Pulse Generator change 5 (31.25) Lead Revision 1 (6.25) Medication Exceptions 1 subject: 1mg Azilect, 9.25 hours before research 1 subject: 1.5 tabs Sinemet 25/100, 11.5 hours before research.

Highlights Filtering and amplitude estimation techniques cause bias estimates of burst activity. Threshold analysis requires (unknowable) a priori knowledge of the fraction in bursts. Support vector machines (SVMs) enable unbiased analysis of burst characteristics. SVM analysis found no differences in burst dynamics with continuous DBS.