Concurrent decoding of distinct neurophysiological fingerprints of tremor and bradykinesia in Parkinson’s disease

  1. Peter M Lauro  Is a corresponding author
  2. Shane Lee
  3. Daniel E Amaya
  4. David D Liu
  5. Umer Akbar
  6. Wael F Asaad  Is a corresponding author
  1. Department of Neuroscience, Brown University, United States
  2. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, United States
  3. The Warren Alpert Medical School, Brown University, United States
  4. Norman Prince Neurosciences Institute, Rhode Island Hospital, United States
  5. Department of Neurosurgery, Rhode Island Hospital, United States
  6. Department of Neurosurgery, Brigham and Women’s Hospital, United States
  7. Department of Neurology, Rhode Island Hospital, United States

Abstract

Parkinson’s disease (PD) is characterized by distinct motor phenomena that are expressed asynchronously. Understanding the neurophysiological correlates of these motor states could facilitate monitoring of disease progression and allow improved assessments of therapeutic efficacy, as well as enable optimal closed-loop neuromodulation. We examined neural activity in the basal ganglia and cortex of 31 subjects with PD during a quantitative motor task to decode tremor and bradykinesia – two cardinal motor signs of PD – and relatively asymptomatic periods of behavior. Support vector regression analysis of microelectrode and electrocorticography recordings revealed that tremor and bradykinesia had nearly opposite neural signatures, while effective motor control displayed unique, differentiating features. The neurophysiological signatures of these motor states depended on the signal type and location. Cortical decoding generally outperformed subcortical decoding. Within the subthalamic nucleus (STN), tremor and bradykinesia were better decoded from distinct subregions. These results demonstrate how to leverage neurophysiology to more precisely treat PD.

Editor's evaluation

This important study advances our understanding of Parkinson's by identifying micro and macro scale signatures linked to critical symptoms (e.g., tremor and slowness of movement), and effective motor control. The evidence supporting the conclusions is solid, and leverages a rich dataset obtained during naturalistic movement. The work will be of interest to neuroscientists, neurologists, and biomedical engineers.

https://doi.org/10.7554/eLife.84135.sa0

Introduction

Parkinson’s disease (PD) is a common and complex neurodegenerative disorder characterized by the dynamic expression of particular motor features such as tremor and bradykinesia (Armstrong and Okun, 2020; Parkinson, 2002). These distinct motor signs are expressed variably across patients and may respond differently to dopamine replacement therapy; their differential expression is often used to classify patients into phenotypic subtypes (Koller, 1986; Sethi, 2008). Despite this heterogeneity, both of these motor features (and both tremor-dominant and non-tremor-dominant patient subtypes) respond to high-frequency deep brain stimulation (DBS) applied to the subthalamic nucleus (STN) (Katz et al., 2015; Limousin et al., 1998).

DBS delivered in a closed-loop fashion (i.e., in response to neurophysiological biomarkers) has shown promising therapeutic potential primarily toward alleviating bradykinesia (Little et al., 2016a; Little et al., 2016b), but current efforts focusing on β frequency oscillations (15–30 Hz) have been shown to inadequately treat or worsen tremor in some cases (Piña-Fuentes et al., 2020; Velisar et al., 2019). Thus, tremor may be better signaled by different components within the local field potential (LFP) spectrum, and closed-loop DBS could benefit from a clearer understanding of the neurophysiological biomarkers that differentiate these motor signs from each other, and from more optimal motor performance in the absence of these impairments.

To this point, STN LFP recordings from patients with different PD subtypes have revealed distinct patterns of oscillatory activity (Telkes et al., 2018). In addition to spectral variability, specific stimulation sites within the STN have been associated with the preferential reduction of individual motor signs (Akram et al., 2017). Moreover, these STN sites were associated with specific patterns of anatomical connectivity with cortical structures (Haynes and Haber, 2013). Much like how overlapping subdivisions of basal-ganglia-cortical circuits have been found to encode separate aspects of movement (Mosher et al., 2021; Neumann et al., 2018), separate motor features may be mediated by different sub-circuits involving the STN and sensorimotor cortex (Gibson et al., 2021).

In order to better reveal the functional and anatomical substrates of distinct PD motor states, we enlisted patients with PD undergoing awake DBS electrode implantation to perform a continuous visual-motor task that allowed rigorous, concurrent measurement of different motor metrics while we acquired STN (micro- and macroelectrode) and cortical (electrocorticography [ECoG]) recordings. Prior studies have not attempted to simultaneously decode different aspects of disease expression, contrast these measures with symptom-free performance, and examine disease expression on the short timescales relevant to that varying expression. While our group has previously demonstrated the ability to decode global PD motor dysfunction from STN recordings on short timescales (Ahn et al., 2020; Sanderson et al., 2020), we focus here on individual motor features and their specific neurophysiological manifestations. Specifically, we trained machine learning models to directly decode tremor or slowness from neural recordings to reveal the spectral and anatomical fingerprints of these cardinal motor features of PD.

Results

Motor behavior during the target tracking task

Twenty-seven patients with PD undergoing STN DBS implantation and 17 age-matched controls performed a visual-motor task in which they followed an on-screen target with a cursor controlled by either a joystick or a stylus and tablet (Figure 1A). Twenty-three patients (and 12 control subjects) performed a version of the task with fixed patterns of target movement, while four patients (and five control subjects) performed a version with randomly generated target paths. Each patient performed 1–4 sessions of the task during the procedure for a total of 69 sessions, while control subjects each performed 1 session extra-operatively for a total of 17 sessions. Tremor amplitude and cursor speed – task metrics calculated to reflect the expression of tremor and bradykinesia – were quantified from the cursor traces. These behavioral metric data were then averaged into 7 s non-overlapping epochs. To compare metrics across subject populations while considering epochs, trials, and sessions as repeated measurements within individuals, linear mixed models (LMMs) were used (see Materials and methods). The resulting metric distributions for PD vs. control subjects demonstrated increased tremor for PD patients (n=6498 epochs across 44 subjects, PD vs. control, LMM coefficient/LMM β=0.337, Z=2.169, p=0.030), but only a trend for decreased speed (PD vs. control, LMM β=−0.592, Z=−1.194, p=0.232) (Figure 1B and C).

Tremor and movement speed calculated from fixed- and random-pattern intraoperative visual-motor tasks.

(A) Left: Schematic of task target (green) and cursor (gray) traces from a single trial of the fixed- (top) or random- (bottom) pattern task. Center-top: Bandpass filtered cursor traces from a task trial. Ca refers to the amplitude of the analytic signal (a) of the cursor trace (C). Center-bottom: Lowpass filtered cursor traces from a task trial. Right-top: One-dimensional projection of bandpass filtered traces (black), with tremor amplitude measured from the envelope (orange). Right-bottom: Cursor speed measured from lowpass filtered traces (black). Figure adapted from Figure 1 of Ahn et al., 2020. (B, C) Distributions of 7 s tremor amplitude (top) and cursor speed (bottom) epochs for control subject and Parkinson's disease (PD) patient populations in the fixed-pattern (B) (n=5375 epochs across 35 subjects) and random-pattern (C) (n=1123 epochs across 9 subjects) task. ° – degrees of visual angle. (D, E) Task-based tremor amplitude (D) and slowness (E) corresponded to UPDRS measures of tremor (D) or bradykinesia (E) (n=24 subjects). ρ=Spearman correlation statistic.

Tremor distributions compiled across task versions revealed that while subjects with PD spent a substantial fraction of time without tremor, they also exhibited a large range of tremor expression not present in control subjects (Figure 1B and C, top). On the other hand, the two task versions generated different movement speed distributions (Figure 1B and C, bottom). While the fixed-pattern version of the task elicited a bimodal distribution (reflecting slower turns and faster straight path segments) in control subjects, the random-pattern version elicited a single peak corresponding to the fixed target speed used in that task. Nonetheless, in the random-pattern version, the PD cursor speed distribution was shifted to the left (i.e., slower) relative to control subjects (n=1123 epochs across 9 subjects, PD vs. control, LMM β=−1.004, Z=−2.210, p=0.027). Cursor speed was converted to ‘slowness’ in order to control for target trajectory/speed variability by normalizing each session’s distribution of cursor speed to its minimum and maximum values (0: highest speed, 1: lowest speed). To determine whether these two motor metrics reflected distinct components of PD motor dysfunction, each patient’s metric distribution medians were correlated against their UPDRS III motor subscores. Indeed, tremor amplitude positively correlated with the resting tremor subscore (ρ=0.607, p=0.001, n=24 PD subjects, Spearman correlation) (Figure 1D) and slowness positively correlated with the hand open/close subscore (a subscore used in part to assess bradykinesia) (ρ=0.406, p=0.035, Spearman correlation) (Figure 1E). However, the opposite correlations were not significant (tremor amplitude – UPDRS hand open/close: ρ=0.022, p=0.446; slowness – UPDRS resting tremor: ρ=0.120, p=0.219). In addition, action/postural tremor subscores trended toward a positive correlation with tremor amplitude (ρ=0.354, p=0.057), while resting and action/postural tremor subscores were positively correlated within subjects (ρ=0.602, p=0.002). Although postural tremor can correlate with resting tremor when patients with PD are measured by the UPDRS for the former and the Washington Heights-Inwood Genetic Study of Essential Tremor (WHIGET) Rating Scale for the latter, resting tremor is thought to be more specific to the PD pathophysiology (Louis et al., 2001). Therefore, tremor and movement speed were considered to reflect two key aspects of PD motor dysfunction.

Tremor and slowness were distinct and opposing symptomatic states

Relative to each other, tremor and slowness typically did not co-occur but rather were inversely expressed in time (n=27 subjects, tremor × slowness, LMM β=−0.584, Z=−19.351, p=2.00*10–83) (Figure 2A and B). To understand whether this anti-correlation may have been due in part to motor features manifesting on different timescales, autocorrelograms were computed for each metric with 100 ms epochs. Here, tremor was typically expressed continuously for longer periods (autocorrelogram full-width half-maximum [FWHM], 0.898 s) as opposed to slowness (0.297 s) (Figure 2C). Using this FWHM as the minimum, we calculated periods of time where metrics were sustained above control levels (i.e., symptomatic periods) across subjects. The median duration of symptomatic tremor episodes was 2.000 s, and slowness episodes lasted for 0.500 s (Figure 2D). With these differing timescales and anti-correlated presence, tremor and slowness appeared to represent distinct symptomatic states.

Tremor and slowness represented two non-overlapping motor states with differing timescales.

(A) Examples from three individual subjects of cursor (solid lines) and target (translucent lines) traces (top row) and calculated motor metrics (bottom three rows) within single trials. Periods of increased expression of individual motor metrics are highlighted by their respective color. (B) (Left) Scatter plot of all cursor speed and tremor measurements in 7 s epochs across subjects. (Right) Histogram of subject-wide behavioral Spearman correlation with tremor and slowness metrics (n=27 subjects). (C) Autocorrelograms of symptomatic (tremor, slowness) and non-symptomatic (effective motor control) metrics. Colored vertical dashed lines indicate full-width half-maximum (FWHM) for each metric. Top-left inset depicts a zoomed-in window of the autocorrelogram. (D) Histogram of sustained motor metric period duration (i.e., symptomatic state duration) across subjects. Solid lines indicate gamma distribution fit to each motor metric state histogram, while dashed vertical lines indicate the median symptomatic state length for each metric.

Because PD produces a fluctuating motor deficit such that there can be moments of normal-appearing motor behavior (Mazzoni et al., 2007), we labeled epochs without motor dysfunction as ‘effective’ motor states. Specifically, epochs with lower tremor and/or higher movement speeds were assigned values closer to 1 while more symptomatic epochs (high tremor and/or slower movement speeds) were assigned values closer to 0. Compared to other metrics, effective motor control was expressed on longer timescales (FWHM = 3.784 s, median state length = 7.900 s) (Figure 2A, C and D).

Tremor and slowness had distinct representations within the STN

A total of 203 microelectrode and 176 macroelectrode recordings (microelectrode tips and macroelectrode contacts separated by 3 mm on the same electrodes) were acquired from the STN as patients performed the task. To assess whether tremor or slowness could be decoded from these recordings, spectral estimates of power from 3 to 400 Hz were obtained using a wavelet convolution. Narrowband power estimates were grouped into six broad frequency bands (θ/α,β,γlow,γmid,γhigh,hfo) with 7 sub-bands each, for a total of 42 neural ‘features’ per 7 s epoch (Ahn et al., 2020). Neural decoding models (support vector regression [SVR] with a linear kernel and 100-fold cross-validation) were trained directly on the epoch’s average metric (tremor or slowness values averaged within each epoch), and their performance was assessed with squared Pearson’s r (r2) between observed and decoded metrics. Across each subject’s best-performing microelectrode recordings (MER), tremor decoding performance (r2=0.232±0.200) was superior to slowness decoding (r2=0.125±0.108) (n=203 MER models, tremor v. slowness, LMM β=0.051, Z=5.477, p=4.33*10–8). No such difference was observed across macroelectrode recordings (n=176 macroelectrode recordings models, tremor v. slowness, r2=0.209±0.174 v. r2=0.198±0.147, LMM β=−0.005, Z=−0.496, p=0.620). To determine if tremor and slowness had distinct neurophysiological signatures, SVR model feature weights were aggregated for each metric. To understand which spectral features were used consistently across models, feature weights were compared to null distributions generated from models where motor metric values were shuffled with respect to the corresponding spectral features. Microelectrode tremor decoding models positively weighted low-frequency features (θ,α,β; 4–21 Hz; p<0.001, permutation test, see Materials and methods). Hfo (275–375 Hz) weights were also positively associated with tremor decoding (p<0.014, permutation test). In contrast, macroelectrode tremor decoding models negatively weighed β power (14–41 Hz; p<0.026, permutation test) while positively weighing γ/hfo activity (60–375 Hz; p<0.011, permutation test). In other words, optimal macroelectrode tremor decoding relied on decreased β power and increased γ/hfo power.

For slowness, microelectrode decoding models had negative θ, γlow , and hfo weights (5–12 Hz, 33–56 Hz, 200–375 Hz) (p<0.012, permutation test). Macroelectrode decoding models positively weighted β frequencies (12–30 Hz; p<0.006, permutation test) along with negative γ/hfo weights (33–375 Hz; p<0.001, permutation test). Tremor and slowness model features differed when compared directly, with hfo frequencies (225–375 Hz) being elevated during tremor in both micro/macroelectrode recordings (Figure 3A). Overall, just as tremor and slowness represented two distinct, anti-correlated symptomatic states of PD, tremor and slowness decoding models from the STN revealed distinguishable patterns of underlying neural activity.

Subthalamic tremor decoding models emphasized lower frequencies whereas slowness models emphasized higher frequencies.

(A) Average tremor decoding and slowness model coefficients for all subthalamic nucleus (STN) microelectrode (left) (n=203 microelectrode recordings, 27 subjects) and macroelectrode (right) recordings (n=176 macroelectrode recordings, 27 subjects). Solid lines indicate average weights, with positive/negative values reflecting a positive or negative relationship with the metric. Error bars indicate s.e.m. across subjects. Black lines (top) represented contiguous spectral features that significantly differed between tremor and slowness decoding models. (B) Average model coefficients for effective motor control and tremor for all STN microelectrode (left) and macroelectrode (right) recordings. (C) Average model coefficients for effective motor control and slowness for all STN microelectrode (left) and macroelectrode (right) recordings.

However, in order to rule out the possibility that the alternating patterns of relevant neural decoding features simply reflected the anti-correlated nature of tremor and slowness, we tested whether decoding models trained for tremor could accurately decode slowness. When directly comparing tremor and slowness decoding performance on tremor-trained models, slowness decoding was inferior for both microelectrode (tremor v. slowness decoding, r2=0.232±0.197 v. 0.002±0.001; LMM β=0.101, Z=11.242, p=2.54*10–29) and macroelectrode (tremor v. slowness decoding, r2=0.205±0.172 v. 0.002±0.001; LMM β=0.086, Z=10.294, p=7.53*10–25) recordings. If decoding features for tremor and slowness were simply inverted, applying models for decoding another metric would result in significantly negative r2-values. Thus, the neural features used for individual metric decoding likely reflected a unique spectral state or ‘fingerprint’.

To further validate that our approach was able to decode motor dysfunction in a symptom-specific fashion, we examined the relationship between individual tremor expression and tremor decoding performance (as not all patients with PD exhibit tremor). Here, we found that task-based tremor distribution medians positively correlated with individual’s highest decoding performance (n=24 patients, ρ=0.442, p=0.031). While UPDRS measures of tremor did not directly correlate with decoding performance (p<0.878) – likely due to a variety of factors distinguishing standard UPDRS tremor assessment from intraoperative task performance – the dynamic, continuous, low-velocity-biased expression of tremor on the naturalistic task (which itself did correlate with UPDRS resting tremor scores) provided an immediate, ground-truth behavior for patient-specific neurophysiological decoding.

Effective motor control had characteristic neural signatures

Effective motor control was similarly decoded from both micro- (r2=0.140±0.104) and macroelectrode (r2=0.204±0.097) recordings. Effective motor control decoding was characterized by the absence of β (10–28 Hz) power in both micro- and macroelectrode recordings (p<0.006, permutation test), while macroelectrodes also exhibited positive γ power weights (30–175 Hz; p<0.020, permutation test). Power in γlow frequencies (30–48 Hz) in particular was significantly increased during effective motor control relative to both tremor and slowness decoding models (p<0.006, permutation test) (Figure 3B–C, right).

In total, STN activity contained specific features that distinguished symptomatic from non-symptomatic motor states. Tremor was characterized by lower frequencies (θ/α) in microelectrodes, slowness by β frequencies in macroelectrodes, and effective motor performance was uniquely characterized by γlow frequencies from both recording types.

Full-spectrum neural decoding outperformed beta-band decoding

To directly test whether each behavior model used neural features across the spectrum, we compared the relative ability of full-spectrum and canonical band (β, 12–30 Hz) models. Full-spectrum decoding had significantly greater performance for macroelectrode (full vs. beta-only decoding, LMM β=0.018–0.035, Z=2.388–3.949, p<0.017) and microelectrode (full vs. beta-only decoding, LMM β=0.014–0.017, Z=2.241–3.154, p<0.025) for all three metrics, with the exception of microelectrode-tremor decoding (full vs. beta-only decoding, LMM β=0.014, Z=1.548, p=0.122).

Optimal subthalamic tremor decoding sites were dorsolateral to optimal slowness decoding sites across patients

To investigate whether tremor and slowness were more optimally decoded from distinct areas within the STN, recording sites for each session were reconstructed using subject-specific neuroimaging (peak MER density in MNI space: x=-12,y=-10,z=-6.0) (Figure 4A; peak macroelectrode recording density: x=-12,y=-9,z=-3.0). For each recording site, the corresponding decoding model performance for each metric was plotted (Figure 4B and C).

Optimal subthalamic tremor decoding sites were dorsolateral to optimal slowness decoding sites.

(A) Recording density of stationary microelectrode recordings across patients (n=182 microelectrode recording sites, 25 subjects) and task sessions overlaid on an MNI reference volume (approximate outline of the subthalamic nucleus [STN] in bolded black, zona incerta outlined above, substantia nigra outlined below). L: left. y-value corresponds to coronal slice in MNI space. (B) Tremor decoding model r2-values for stationary microelectrode recordings. (C) Slowness decoding model r2-values for stationary microelectrode recordings. (D) Difference in tremor vs. slowness decoding r2-values for stationary microelectrode recordings. Warmer colors indicate voxels where tremor decoding was superior, whereas cooler colors indicate where slowness decoding was superior. (E) Recording density of moving microelectrode recordings across all patients and task sessions overlaid on an MNI reference volume. (F) Tremor decoding model r2-values for high-density STN survey recordings. (G) Slowness decoding model r2-values for high-density STN survey recordings. (H) Difference in tremor vs. slowness decoding r2-values for high-density STN survey recordings. r2-Values depicted here are site-specific r2-values generated from the whole-STN model applied to individual depth recordings. Warmer colors indicate voxels where tremor decoding was superior, whereas cooler colors indicate where slowness decoding was superior.

We then compared the voxel-wise relative performance between tremor and slowness throughout all recorded STN voxels by using a modified 3D t-test with spatially based permutation shuffling (see Materials and methods). Tremor was better decoded in recordings from dorsolateral STN (n=182 MER sites, x=-14.0,y=-13.0,z=-5.0; Z=2.116,p=0.017), whereas slowness was better decoded from recordings in central/ventromedial STN (x=-12.0,y=-14.0,z=-6.0; Z=1.911,p=0.028) (Figure 4D).

Optimal locations for tremor and slowness decoding were not found to differ significantly by macroelectrode location (n=176 macroelectrode recording sites, p>0.05). Moreover, the locus of optimal effective motor control decoding was not observed to differ from those of tremor or speed using either micro- and macroelectrode recordings (p>0.05). Nevertheless, we found that that differences in metric decoding were not only related to the frequencies present, but also to an electrode’s location within the STN, as assessed over the entire study PD population.

Optimal subthalamic tremor decoding sites were dorsolateral to optimal slowness decoding sites within individual patients

To verify the spatial relationship of optimal tremor and slowness decoding within patients, five additional right-handed patients (70.0±8.9 years of age; 2F, 3M; UPDRS III: 45.2±9.5) underwent a modified version of the random-pattern task. Rather than acquiring recordings from a stationary site, here we surveyed the entire length of the STN by systematically moving the electrodes between task trials in small, discrete steps using automatic, computer-driven microdrive control (see Materials and methods, High-density STN survey).

SVR models for tremor and slowness were then calculated by incorporating recording data across all sites/trials within a single trajectory. Although decoding performance of models derived from multi-site data exhibited a trend of lower performance than models trained on single-site recordings (tremor: r2=0.232±0.200 v. 0.073±0.052; slowness: r2=0.125±0.108 v. 0.061±0.079, effective motor control: r2=0.140±0.104 v. 0.043±0.058), these differences were not significant (n=203 stationary MER sites, n=17 moving MER trajectories, moving v. stationary data, LMM β=−0.075 to –0.075, Z=−0.945 to –1.291, p=0.197–0.344). Despite less data at each recording site, whole-STN models demonstrated above-chance decoding performance for all three metrics (tremor: 8/17 trajectories, slowness: 6/17, effective motor control: 6/17).

Recording sites along each trajectory were reconstructed using imaging (Figure 4E), and site-specific metric decoding r2-values were calculated by applying the whole-STN SVR model to individual site recordings (Figure 4F and G) (see Materials and methods). Decoding performance was then compared across patients (Figure 4H). Across these recordings, tremor was optimally decoded at (x=13.0,y=13.0,z=5.0; Z=1.911,p=0.028), while slowness was optimally decoded at (x=12.0,y=15.0,z=8.0; Z=1.937,p=0.026). Within individual subjects, tremor was again found to be decoded dorsolaterally to slowness.

Cortical recordings also revealed distinct representations of tremor, slowness, and effective motor control

Ten subjects additionally had ECoG recordings from sensorimotor cortex (motor cortex: n=16 contacts, somatosensory cortex: n=15, see Materials and methods). SVR models for metric decoding were similarly trained on ECoG signals. ECoG decoding performance did not differ between tremor or slowness (n=85 ECoG recordings across 27 sessions and 10 subjects, tremor: r2=0.323±0.153, slowness: r2=0.314±0.143, tremor v. slowness, LMM β=0.010, Z=0.578, p=0.563).

To understand which spectral features contributed to cortical motor metric decoding, SVR model weights were aggregated across all patients and recordings and compared to metric-shuffled models. When compared directly, cortical tremor and slowness models had opposing relationships in α/β (8–40 Hz, p<0.027, permutation test), γmid (45–125 Hz, p<0.023, permutation test), and γhigh (150–225 Hz, p<0.015, permutation test) frequency bands (Figure 5A). Tremor models additionally had positive weights associated with θ frequencies (5–7 Hz, p=0.003, permutation test). Altogether, although cortical signals supported equivalent decoding performance for tremor or slowness, decoding features were nonetheless distinct. On the other hand, effective motor control decoding performance (r2=0.469±0.112) was lower than tremor (effective motor control v. tremor, LMM β=−0.067, Z=−3.975, p=7.04*10–5) and slowness (effective motor control v. slowness, LMM β=−0.057, Z=−3.397, p=0.001). Nevertheless, effective motor control was represented in cortical decoding models by γhigh frequencies. These γhigh features additionally appeared to differentiate effective motor control models from both tremor and slowness models (125–175 Hz, p<0.026, permutation test) (Figure 5B–C). In addition, α/β (8–30 Hz, p<0.001, permutation test) and γlow (45–75 Hz, p<0.010, permutation test) frequencies exhibited an opposing relationship between effective motor control and slowness, much like the interaction between tremor and slowness models (Figure 5C). Taken together, although at different γ frequencies, both STN (γlow) and sensorimotor cortex (γhigh) exhibited features specific to effective motor control. And similarly to STN recordings, ECoG full-spectrum decoding was superior for tremor and effective motor control (full vs. beta-only decoding, LMM β=0.017–0.024, Z=2.224–2.577, p<0.026), but equivalent for slowness (full vs. beta-only decoding, LMM β=−0.013, Z=−1.531, p=0.126).

Cortical tremor and slowness decoding models exhibited opposing weights for multiple frequency bands, and co-expressed specific features with subthalamic recordings.

(A) Average cortical tremor and slowness decoding model coefficients for every recording along sensorimotor cortex (n=85 electrocorticography [ECoG] recordings, 10 subjects). Colored lines indicate average weights, with positive/negative values reflecting a positive or negative relationship with the metric. Error bars indicate s.e.m. across subjects. Black lines (top) represented contiguous spectral features that significantly differed between tremor and slowness decoding models. (B) Average model coefficients for effective motor control and tremor. (C) Average model coefficients for effective motor control and slowness. (D) Average subthalamic nucleus [STN]-cortical coherence tremor and slowness decoding model coefficients for every pairwise recording along sensorimotor cortex and macro contacts within the STN (n=85 ECoG recordings, 10 subjects). (E) Average coherence model for effective motor control and tremor. (F) Average model coefficients for effective motor control and slowness.

Finally, we analyzed whether motor features were selectively represented in different regions of cortex. ECoG recording sites (and their associated metric decoding performance) were plotted along a standard cortical surface (Figure 6A). Comparing tremor and slowness decoding performance by cortical anatomy revealed that slowness decoding had peaks in medial motor (n=31 ECoG recording sites, x=-37.5,y=-15.8,z=+70.2; T=2.210, p=0.030) and somatosensory (x=-34.4,y=-37.4,z=+71.0; T=2.340,p=0.022) cortices. We also observed a trending tremor decoding peak in lateral somatosensory cortex (x=-50.2,y=-31.8,z=+64.0; T=1.700,p=0.093) (Figure 6B). Similar to the STN, effective motor control decoding performance (relative to tremor or slowness decoding) was found not to differ by cortical anatomy (p>0.05). In sum, cortical signals supporting optimal tremor cortical decoding were found relatively lateral to those supporting better slowness decoding.

Cortical tremor and slowness decoding models were distributed throughout cortex, and generally were superior to subthalamic nucleus (STN) decoding models.

(A) Recording density of electrocorticography (ECoG) contacts (n=31 ECoG sites, 10 subjects) on an MNI reference surface. (B) Difference in tremor vs. slowness decoding r2-values for all cortical recordings. Warmer colors indicate surface vertices where tremor decoding was superior, whereas cooler colors indicate where slowness decoding was superior. (C) Decoding performance across metrics and recording types (n=85 ECoG recordings, 81 macroelectrode recordings, 81 microelectrode recordings, 10 subjects). Box represents interquartile range (25th–75th percentile), while whiskers indicate 5th–95th percentile of data range. Brackets indicate significant (*, p<0.05 in linear mixed model comparisons) differences in metric decoding r2-values. (D) Examples of tremor (top) and slowness (bottom) decoding for top-performing ECoG contacts. Each row within the panel represents a different subject. (E) Examples of tremor (top) and slowness (bottom) decoding for top-performing macroelectrodes. (F) Examples of tremor (top) and slowness (bottom) decoding for top-performing microelectrodes.

Cortical decoding outperformed STN decoding

To understand whether motor (dys-)function was better represented in cortical signals, decoding performance was compared between patients with ECoG and STN recordings (n=10 subjects, 85 ECoG recordings, 81 MER, 81 macroelectrode recordings) (Figure 6C). ECoG tremor decoding models exhibited higher performance than micro-STN recordings (ECoG v. micro-STN, r2=0.323±0.153 v. 0.240±0.200, LMM β=0.047, Z=2.504, p=0.012), but only trended higher relative to macro-STN recordings (ECoG v. macro-STN, r2=0.323±0.153 v. 0.294±0.222, LMM β=0.019, Z=1.001, p=0.317). Slowness decoding performance on the other hand was higher in ECoG models relative to both macro-STN (ECoG v. macro-STN, r2=0.314±0.143 v. 0.195±0.105, LMM β=0.061, Z=4.445, p=8.78*10–6) and micro-STN (ECoG v. micro-STN, r2=0.314±0.143 v. 0.128±0.134, LMM β=0.096, Z=7.025, p=2.14*10–12) recordings. Like tremor, effective motor control ECoG decoding models exhibited superior decoding performance to micro-STN (ECoG v. micro-STN, r2=0.233±0.105 v. 0.133±0.096, LMM β=−0.038, Z=−3.140, p=0.002), but not macro-STN (ECoG v. macro-STN, r2=0.233±0.105 v. 0.219±0.098, LMM β=0.004, Z=0.312, p=0.755) recordings. In summary, recordings from sensorimotor cortex were superior to micro- and macro-STN recordings for decoding slowness, while cortical recordings were also superior to micro-STN recordings for decoding tremor and effective motor control. However, both cortical and deep recordings were able to meaningfully track patient symptoms (Figure 6D–G).

Because cortical and subthalamic decoding models used similar features for metric decoding, we investigated whether combining cortical-subthalamic recordings would be superior to using local signals. Pairwise coupling between cortical and macro-STN contacts was calculated via coherence within the same frequency range and sub-band features for SVR decoding. Across all metrics, local decoding models outperformed pairwise decoding models (tremor: ECoG v. EcoG/macro-STN, r2=0.323±0.153 v. 0.075±0.029, LMM β=0.133, Z=15.073, p=2.43*10–51; slowness: ECoG v. EcoG/macro-STN, r2=0.314±0.143 v. 0.098±0.057, LMM β=0.123, Z=14.420, p=3.89*10–47; effective motor control: ECoG v. EcoG/macro-STN, r2=0.233±0.105 v. 0.096±0.050, LMM β=0.065, Z=10.213, p=1.73*10–24). When examining decoding coherence features, slowness models positively weighed frequencies within the β range (12–26 Hz, p<0.0005, permutation test) when compared to tremor and effective motor control decoding models. Effective motor control on the other hand positively weighed frequencies in the γlow range (34–45 Hz, p<0.010, permutation test) when compared to tremor and slowness decoding models. While tremor decoding models trended toward positive weights for θ (5–7 Hz, p<0.235, permutation test), these weights did not survive multiple comparisons. Overall, while cortical-subthalamic pairwise models relied on co-occurring physiological motifs, local neuronal information appeared more relevant to short-timescale decoding of motor dysfunction.

Discussion

In this study we quantified PD tremor and movement speed in a structured motor task as surrogates for coarser clinical measurements of tremor and bradykinesia. We decoded these metrics using linear models and mapped these results to basal ganglia and cortical anatomy. We demonstrated that tremor and bradykinesia were represented by different functional motifs with distinct localization in the STN and sensorimotor cortex. We also contrasted these pathological states with periods of effective motor control, revealing unique markers of relatively symptom-free states. To our knowledge this is the first study to not only characterize the behavioral interaction between tremor and slowness within a single behavioral context but also to compare directly each motor sign’s corresponding expression in neural activity, and further to compare this with relatively asymptomatic states. These results provide a holistic description of dynamic, spontaneous alternations in PD symptoms which reveal specific neurophysiological biomarkers of non-pathological and distinct pathological states.

We focused our neural decoding approach on two cardinal motor features of PD to isolate spectral features that reflected the expression of each. In the STN, tremor was characterized by lower-frequency (θ,α) oscillations in microelectrodes, whereas slowness was characterized by the presence of β oscillations and the absence of γ oscillations in macroelectrodes. Because γ frequency oscillations are commonly associated with hyperkinetic states (Lofredi et al., 2018; Swann et al., 2016), our slowness decoding results may be understood in part as an ‘anti-speed’ neural model. Indeed, effective motor control was distinguished by γlow frequency activity, highlighting the importance of γ frequency oscillations in effective movements. Some of these frequency bands in isolation (θ,β) have been found to correlate with clinical measures of tremor and bradykinesia (Asch et al., 2020; Neumann et al., 2016; Nie et al., 2021). However, here we show directly the contrasting nature of distinct PD motor states both behaviorally and neurophysiologically, complement an evolving literature of dynamic STN states in PD (Khawaldeh et al., 2022), and highlight the dependence of these neurophysiological ‘fingerprints’ on the particular neural recording technique.

We also identified where tremor and slowness were optimally decoded (i.e., where metric-specific spectral information was greatest). Within our STN MER, optimal tremor decoding sites were found to be located within dorsolateral STN whereas optimal slowness decoding sites were more centrally located within the STN. Optimal tremor decoding may have included activity from zona incerta, a stimulation site commonly thought to be critical for alleviating tremor (Plaha et al., 2008; Reck et al., 2009). Indeed, optimal stimulation sites to alleviate tremor and bradykinesia correspond to our dorsolateral-tremor/ventromedial-slowness topography (Akram et al., 2017). While several groups have localized β frequency activity to dorsolateral STN, this has been observed to be located inferiorly to tremor-related higher frequency oscillations (Tamir et al., 2020; Telkes et al., 2018; van Wijk et al., 2017). Our tremor/slowness-dorsal/ventral STN results correspond to prior work suggesting subdomains within the STN. In macaques, motor cortex projects to the dorsal portion of the STN, and the ventral STN receives projections from prefrontal cortex (Haynes and Haber, 2013). In both macaques and humans, the ventral STN has been associated with stopping movement, while the dorsal STN is more associated with motor initiation and selection (Chen et al., 2020; Mosher et al., 2021; Pasquereau and Turner, 2017). In addition, the co-occurrence of tremor-frequency (4–8 Hz; θ) activity in both STN and cortical recordings is consistent with previous work finding tremor-frequency oscillations (θ) originating in the STN and propagating/synchronizing to motor cortex during tremor (Hirschmann et al., 2013; Lauro et al., 2021). Our anatomical results suggest that this propagation may be specific to dorsal STN. Slowness decoding models alternatively relied upon β activity in macroelectrode recordings, perhaps reflecting anti-kinetic β bursts relayed to ventral STN from inferior frontal or supplementary motor cortex (Hannah et al., 2020; Oswal et al., 2021). But while prior work did not directly compare the neuroanatomical substrate of distinct PD features within or across subjects, this work demonstrates how alternating motor features of PD may manifest along these anatomical subdivisions.

In general, cortical recordings were equally capable of decoding tremor or slowness and were generally superior to STN-based decoding as previously reported (Merk et al., 2022b). When comparing the feature weights of these decoding models, we observed opposing relationships in both β and γ bands. As previous studies have shown that tremor decreases β oscillations across cortex (Qasim et al., 2016), and others have shown increased narrowband γ activity during hyperkinetic/dyskinetic states (Swann et al., 2016), here we demonstrate a ‘push-pull’ relationship between these frequency bands in the alternating expression of tremor and slowness, and when comparing slowness and effective motor control models. While cortical β oscillations (and their desynchronization with movement) are well characterized in PD (Rowland et al., 2015), the functional role of broadband γhigh /hfo activity is less clear. Although higher-frequency activity overlaps with phase-amplitude coupling peaks observed in cortex in unmedicated patients with PD (de Hemptinne et al., 2013), our models for effective motor control states suggested that γhigh is specifically associated with more normal movement in line with previous research on human cortical sensorimotor mapping (Crone et al., 1998; Fischer et al., 2017).

With advances in technology, DBS aspires toward incorporating chronic neurophysiological recordings to help guide therapeutic stimulation (Gilron et al., 2021). While current closed-loop DBS paradigms trigger stimulation based on one or two frequency bands representing PD symptoms (Kehnemouyi et al., 2021), our results argue for the potential utility of a more targeted neurophysiological approach to PD state identification: more dorsal STN contacts may better sense signals reflecting tremor, while more ventral STN contacts may better identify signals corresponding to bradykinesia. Precise patient- and symptom-specific models could not only inform where to stimulate, but also when and how to stimulate (i.e., identifying stimulation settings to best treat tremor vs. bradykinesia). Most importantly, future neuromodulation paradigms could be derived not simply to disrupt pathological activity but actually to sustain the neurophysiological γ frequency ‘targets’ of effective motor control. Looking ahead, chronic cortical recordings could work in concert with STN recordings to help identify precise motor states associated with specific aspects of disease expression much like prior work with essential tremor (He et al., 2021; Opri et al., 2020).

These results naturally suggest future prospective intra- or extra-operative studies applying these decoding models toward DBS control policies. While other studies decoding PD motor dysfunction from STN recordings have investigated the use of several decoding techniques (hidden Markov models, logistic regression, Kalman filters), their complexity/computational requirements will require balancing potential therapeutic benefits with the computational/power limitations of implanted pulse generators (Hirschmann et al., 2017; Merk et al., 2022a; Shah et al., 2018; Yao et al., 2020). Although relatively simple, our linear SVR approach was robustly able to fit patient- and symptom-specific decoding models in an interpretable and iterable fashion. Near-future studies using implanted hardware will require adaptation to device constraints, such as reducing the neural feature spectra based on potentially reduced sampling rates (e.g., 3–100 Hz based on a device-based 250 Hz sampling rate), and to conform to potentially fewer available frequency band estimates, though future devices are expected to allow for greater flexibility in signal processing and algorithmic complexity.

While our motor metrics correlated with UPDRS subscores, we recognize that our single intraoperative behavioral task does not capture all aspects of PD motor dysfunction, and is not as naturalistic as other studies using chronic recordings (Hirschmann et al., 2017); nevertheless, our approach had the advantages of being objective, quantitative, and consistent, while providing the unique opportunity to compare different modalities of recordings (e.g., micro- and macroelectrodes). Further, because patients were withdrawn from dopaminergic medication for at least 12 hr prior to surgery, exacerbated patient symptoms in the operating room may produce improved model decoding performance relative to models trained on patients taking dopaminergic medication in the natural setting. Future studies implementing these models in the extra-operative and clinical setting may help bridge these gaps.

While we focused on spectral power measurements in each structure and coherence across structures, future work may examine whether instantaneous phase-based measures of synchrony across structures may potentially better decode motor states. Given the nonuniform spatial sampling of imaging-based reconstructions, our imaging-based analyses may have lacked sufficient power to reveal smaller-scale or additional neurophysiological-anatomical relationships. Because of our method of identifying intraoperative recording sites suitable for the task (identifying motor-responsive single units through clinical somatotopic testing), our recording site distribution and subsequent decoding sub-regions may have been biased toward dorsolateral recordings. As tremor decoding exhibited higher decoding performance than slowness, our observed slowness peak in ventral STN may reflect overall decreased decoding performance. In addition, as macro-STN recordings were necessarily 3 mm above MER, some recordings may have been collected outside of the STN. High-density STN MER may have additionally been impacted by larger macro-STN contacts deforming tissue in the initial downward recording trajectory. Although our decoding performance may be improved by larger datasets and/or more advanced machine learning approaches, our linear approach robustly achieved patient-specific decoding while revealing metric-specific neurophysiological signatures. Nevertheless, with our parametric measurements of PD motor behavior in the intraoperative setting, we were able to delineate the contrasting, push-pull relationship between neural states underlying tremor, bradykinesia, and effective motor control in both the STN and sensorimotor cortex.

Materials and methods

Study and experimental design

Request a detailed protocol

All patients undergoing routine, awake placement of deep brain stimulating electrodes for intractable, idiopathic PD between June 2014 and December 2018 were invited to participate in this study. Patients were selected and offered DBS by a multi-disciplinary team based solely upon clinical criteria (Akbar and Asaad, 2017). In this report, 32 subjects (n=27 subjects in stationary-recording experiments, n=5 subjects in high-density recording experiments) (4F, 28M; aged 47.5–78.5 years) underwent STN DBS. Subjects were off all anti-Parkinsonian medications for at least 12 hr in advance of the surgical procedure (UPDRS Part III: 48.3 ± 12.6). Twelve subjects were considered tremor-dominant, and 13 subjects had average tremor UPDRS III scores >2 in their dominant hand (Jankovic et al., 1990). Seventeen age-matched controls (14F, 3M; often patients’ partners; aged 48.5–79.2 years) also participated in this study (patient v. control subjects ages, 65.2±7.9 vs 63.2±9.3, Mann-Whitney U-test, p=0.450). Controls were required simply to be free of any diagnosed or suspected movement disorder and to have no physical limitation preventing them from seeing the display or using the manipulandum. There was a strong male bias in the patient population (2F, 25M) and a female preponderance in the control population (14F, 3M), reflecting weaker overall biases in the prevalence of PD and the clinical utilization of DBS therapy (Accolla et al., 2007; Hariz et al., 2011; Rumalla et al., 2018). Both patients and control subjects were predominantly right-handed (patients: 30 right-handed, 2 left-handed; control subjects: 16 right-handed, 1 left-handed). Patients and control subjects agreeing to participate in this study signed informed consent, and experimental procedures were undertaken in accordance with an approved Rhode Island Hospital human research protocol (Lifespan IRB protocol #263157) and the Declaration of Helsinki.

Surgical procedure

Request a detailed protocol

MER from the region of the STN of awake patients are routinely obtained in order to map the target area and guide DBS electrode implantation. Microdrives (Alpha Omega Inc, Nazareth, Israel) were attached to a patient-customized stereotactic platform (STarFix micro-targeting system, FHC Inc) and then loaded with three parallel microelectrodes (Konrad et al., 2011). For 10 patients, ECoG strips were placed posteriorly along sensorimotor cortices through the same burr hole used for MER insertion to conduct intraoperative cortical recordings. The STN was identified electrophysiologically as a hyperactive region typically first encountered about 3–6 mm above estimated target (Gross et al., 2006). When at least one electrode was judged to be within the STN, electrode movement was paused and recordings were obtained in conjunction with patient performance of the visual-motor task.

High-density STN survey

Request a detailed protocol

In five subjects, once the bottom of the STN was identified using typical electrophysiological procedures, custom-built routines using an FDA-approved software development kit (Alpha Omega, Inc) were used to automatically raise electrodes by a pre-specified distance between trials to conduct a high-density STN survey. To start, clinical MER was conducted in typical fashion. Once the electrodes were judged to have exited the STN, the length of the STN recording span was calculated based on intraoperative neurophysiology. Based upon this length, the electrodes were automatically raised by the microdrives in pre-calculated steps in coordination with the visual-motor task, during the inter-trial intervals. During this task, a separate control computer was used to coordinate the behavioral task with robotic control of the Alpha Omega neurophysiology and microdrive systems. Specifically, the FDA-approved C++ Neuro Omega software development kit was compiled into a custom Python library that could communicate with the Neuro Omega systems with a ∼2 ms round-trip latency. From there, task-specific Python code enabled coordination with the behavioral control system. To acquire MER that spanned the STN, the length of the STN was estimated based upon standard neurophysiological assessment, and this length was divided by the number of task trials (typically 36). As the task was performed, the start of each inter-trial interval was detected by the control computer, and every few trials (typically 3), a command was issued to raise the electrodes by the appropriate distance. The task re-commenced once drive movement was complete (typically ∼10 s later). This process continued until the subject completed the task and the microelectrodes had reached the top of the STN.

Behavioral task and metrics

Request a detailed protocol

We employed a visual-motor target tracking task to estimate motor dysfunction in a quantitative and continuous fashion using MonkeyLogic (Asaad et al., 2013; Asaad and Eskandar, 2008a; Asaad and Eskandar, 2008b; Hwang et al., 2019). As per standard surgical procedure, patients were positioned in a reclined ‘lawn-chair’ position, supine on the operating table to maintain comfort while allowing patients to engage with the task and clinical assessment. For the task, a boom-mounted display was adjusted to the patient’s line of sight, and patients were asked to verbally confirm their ability to see on-screen task objects. Patients using a joystick had it placed in their lap, while patients using a tablet had it placed on a stand in their lap (angle adjusted to comfort) and a stylus placed in their dominant hand. Healthy control subjects performed the task by sitting in a chair at a table, with the manipulandum secured on the table. Similar to the patients, the adjustable arm-mounted task display was adjusted to patient’s line of sight. Subjects were instructed to follow a green target circle that moved smoothly around the screen by manipulating the joystick or stylus with their dominant hand with the goal of keeping the white cursor within the circle (Figure 1A). All subjects were instructed to not rest their dominant hand on their lap/the table while performing the task. The target circle followed one of several possible paths (invisible to the subject), with each trial lasting 10–30 s. Each session consisted of up to 36 trials (∼13 min of tracking data).

Tremor amplitude was calculated from 3 to 10 Hz bandpass filtered cursor traces, TMxt=Cax . Movement speed was calculated from cursor traces lowpass filtered at 3 Hz to remove the influence of tremor, Speedx=ΔxΔt . Both metrics were averaged into 7 s non-overlapping epochs to maintain consistency with our previous decoding approach (Ahn et al., 2020). To standardize movement speed within subjects, movement speed epochs within a session were min-max normalized into a measure of ‘slowness,’ where 0=highest speed and 1=lowest speed.

Effective motor control was quantified as the absence of tremor and slowness measures, relative to the entire session. Each epoch’s ‘effective motor control’ measure was then calculated as 1-Tremor+1-Slowness2 , where values of 0 indicated symptomatic states (tremor, slowness) whereas values of 1 indicated optimal motor performance.

Tremor and slowness were compared across control and PD populations using the following LMM: ymetric=Xpopulationβ+Zu+ϵ, where ymetric represented each epoch’s metric amplitude and Xpopulation represented categorical labels of populations. LMMs were used to calculate the correlation between tremor and slowness across the entirety of each subject’s behavioral data using the following model: ytremor=Xslownessβ+Zu+ϵ, where ytremor represented each epoch’s tremor amplitude and Xslowness represented the epoch’s simultaneous measurement of slowness.

To determine the timescale of metric fluctuation, autocorrelograms were calculated across each PD subject’s behavioral data using 100 ms epochs. The average FWHM of the autocorrelograms were considered the minimum time necessary to label motor metric data as a ‘symptomatic’ period. Tremor or slowness were considered ‘symptomatic’ if they exceeded the 95th percentile of aggregate control data, and sustained symptomatic periods were defined as those persisting beyond the population metric FWHM continuously. For effective motor control, epochs were labeled ‘symptomatic’ if they were above the median of the PD subject’s session distribution.

LMM design

Request a detailed protocol

Behavioral metrics (tremor and slowness) were compared across control and PD populations using an LMM to account for each subject’s asymmetric contribution of epochs: ymetric=Xpopulationβ+Zu+ϵ, where ymetric represented each epoch’s metric amplitude, Xpopulation represented categorical labels of populations and associated fixed-effect regression coefficients (β), Z represented the subject-specific random intercepts and their associated random effect coefficients (u), and ϵ represented the residuals. To understand the interactions and optimal timescales of tremor and slowness, each metric was calculated within smaller, 100 ms epochs. LMMs were used to calculate the correlation between tremor and slowness across the entirety of each subject’s behavioral data using the following model: ytremor=Xslownessβ+Zu+ϵ, where ytremor represented each epoch’s tremor amplitude, Xslowness represented the epoch’s simultaneous measurement of slowness and the fixed-effect regression coefficient (β), and Zu represented the same subject-specific random intercept design as above.

When assessing whether one type of metric (e.g., tremor) was preferentially decoded within a single type of recording (e.g., microelectrodes), SVR r2-values were compared using the following LMM: yr-value=Xmetricβ+Zu+ϵ, where yr-value represented SVR decoding r2-values from a single recording and metric, Xmetric represented categorical labels of metrics and associated fixed-effect regression coefficients (β), Z represented the subject-specific random intercepts and their associated random effect coefficients (u), and ϵ represented the residuals. When investigating whether one type of recording was superior at decoding a single metric, r2-values were compared using the following LMM: yr=Xrecordingβ+Zu+ϵ, where yr represented SVR decoding r2-values from a single recording and metric, Xrecording represented categorical labels of recording types and associated fixed-effect regression coefficients (β), and Zu represented the same subject-specific random intercepts model as above.

When comparing cross-metric model performance (i.e., determining the ability of a model trained on tremor to decode slowness), performance was assessed by linear regression between the model’s predicted metric (tremor) distribution and the co-occurring alternate metric (slowness) distribution. To compare the relative performance of tremor-trained models on decoding slowness, r2-values were compared within recording type using the following LMM: yr=Xmetricβ+Zu+ϵ, where yr represented the SVR decoding r2-value, Xmetric represented the categorical labels of either the model’s trained metric (tremor) or the alternate metric (slowness) and their associated fixed-effect regression coefficients (β), Z represented the subject-specific random intercepts and their associated random effect coefficients (u), and ϵ represented the residuals.

For datasets collected using the within-subject, high-density STN survey, SVR models were trained using recordings throughout the STN. Specifically, recordings at each depth were split into 2:1 train:test sets and aggregated for whole-STN SVR model fitting. r2-Values from SVR models trained on high-density microelectrode data were compared to r2-values from stationary microelectrode data using the following LMM: yr=Xexperimentβ+Zu+ϵ, where yr represented the SVR decoding r2-value, Xexperiment represented the categorical labels of experiment type and their associated fixed-effect regression coefficients (β), Z represented the subject-specific random intercepts and their associated random effect coefficients (u), and ϵ represented the residuals.

Neurophysiological signals and analysis

Request a detailed protocol

Microelectrode signals were recorded using ‘NeuroProbe’ tungsten electrodes (Alpha Omega), and macroelectrode signals were recorded from circumferential contacts 3 mm above the microelectrode tips. ECoG signals were acquired using Ad-Tech 8-contact subdural strips with 10 mm contact-to-contact spacing (Ad-Tech Medical, Racine, WI, USA). All signals were acquired at 22–44 kHz and synchronized using Neuro Omega data acquisition systems (Alpha Omega). Patients performed up to four sessions of the task, with microelectrodes positioned at different depths for each session. As microelectrodes were not independently positionable, some signals may have necessarily been acquired outside of the STN. All recorded signals were nevertheless considered and analyzed.

Neural data from the hemisphere contralateral to the patient’s dominant hand were analyzed using the ‘numpy/scipy’ Python 3 environment (Harris et al., 2020; Virtanen et al., 2020). Offline, ECoG signals were re-referenced to a common median reference within a strip (Liu et al., 2015). All resulting signals were bandpass filtered between 2 and 600 Hz, and notch filtered at 60 Hz and its harmonics. These resulting timeseries were then downsampled to 1 kHz. Timeseries were bandpass filtered using a Morlet wavelet convolution (wave number 7) at 1 Hz intervals, covering 3–400 Hz. The instantaneous power and phase at each frequency was then determined by the Hilbert transform. To analyze broad frequency bands, we grouped frequencies into canonical ranges: θ/α : 3–12 Hz, βlow : 12–20 Hz, βhigh : 20–30 Hz, γlow : 30–60 Hz, γmid : 60–100 Hz, γhigh : 100–200 Hz, and hfo (high-frequency oscillations): 200–400 Hz. Power within the hfo band was interpreted as multiunit spiking activity, rather than discrete oscillations/ripples.

For cortical-subthalamic pairwise decoding models, neural synchrony was quantified using Welch’s magnitude-squared coherence between timeseries within 7 s epochs (Hann windows, 1024 samples per segment with 512 sample shifts). The resulting coherence spectra were averaged into the same frequency bands as above.

Imaging-based reconstruction of recording sites

Request a detailed protocol

Preoperatively, magnetic resonance (MR) images were obtained that included T1- and T2-weighted sequences (T1: MPRAGE, T2: SPACE; Siemens Vario 3.0T scanner). Pre-, intra-, and postoperative (in some cases) computed tomography (CT) scans were also acquired (Extra-Op CT: GE Lightspeed VCT Scanner; Intra-Op CT: Mobius Airo scanner). Postoperative T1-weighted MR images were typically obtained 1–2 days after the operation. To reconstruct recording locations, MR and CT images were co-registered using the FHC Waypoint Planner software and AFNI (Cox, 1996; Li et al., 2016). Microelectrode depths were calculated by combining intraoperative recording depth information with electrode reconstructions obtained from postoperative images using methods described previously (Lauro et al., 2018; Lauro et al., 2016). To determine the anatomical distribution of MER sites across patients, preoperative T1-weighted MR images were registered to a T1-weighted MNI reference volume (MNI152 T1 2009c) (Fonov et al., 2009). The resulting patient-specific transformation was then applied to recording site coordinates, which were then assessed for proximity to the STN as delineated on the MNI PD25 atlas (Xiao et al., 2017; Xiao et al., 2015; Xiao et al., 2012). ECoG contacts were segmented from intraoperative CT volumes, and were then projected onto individual cortical surface reconstructions generated from preoperative T1 volumes (Dale et al., 1999; Fischl et al., 2002; Saad and Reynolds, 2012; Trotta et al., 2018). Individual cortical surface reconstructions were co-registered to a standard Desikan-Destrieux surface parcellation (Argall et al., 2006; Desikan et al., 2006; Destrieux et al., 2010). Contacts within sensorimotor cortex (labeled as motor or somatosensory cortex by parcellation label – ‘ctx_G_precentral’ and ‘ctx_G_postcentral’, respectively) were considered for the present study.

Neural decoding of behavioral metrics

Request a detailed protocol

To investigate whether STN or cortical activity could be used to estimate co-occuring behavioral metrics, SVR with a linear kernel using ‘scikit-learn’ was applied toward multi-spectral decoding of tremor or slowness (Pedregosa et al., 2011). Spectral power estimates for each canonical band (θ/α , βlow , βhigh , γlow , γmid , γhigh , hfo) were further subdivided into 7 sub-bands for a total of 42 spectral features across 3–400 Hz (Ahn et al., 2020). SVR models trained on a single electrode’s spectral features were fit using 100-fold Monte Carlo cross-validation with a 2:1 train/test split of temporal epochs within a task session. Model performance was assessed by linear regression (specifically, the squared Pearson r-value – r2) between the observed and predicted metric distributions. To verify that these decoded results were not spurious, a separate set of SVR models were fit with a shuffled correspondence between behavioral metric data and neurophysiological signals in the training set.

When assessing whether one type of metric (e.g., tremor) was preferentially decoded within a single type of recording (e.g., microelectrodes), SVR r2-values were compared using the following LMM: yr-value=Xmetricβ+Zu+ϵ, where yr-value represented SVR decoding r2-values from a single recording and metric and Xmetric represented categorical labels of metrics. When investigating whether one type of recording was superior at decoding a single metric, r2-values were compared using the following LMM: yr=Xrecordingβ+Zu+ϵ, where yr represented SVR decoding r2-values from a single recording and metric and Xrecording represented categorical labels of recording types.

When comparing the relative ability of a model trained on tremor to decode tremor or slowness, r2-values were compared within recording type using the following LMM: yr=Xmetricβ+Zu+ϵ, where yr represented the SVR decoding r2-value and Xmetric represented the categorical labels of either the model’s trained metric (tremor) or the alternate metric (slowness).

Because these SVR models used a linear kernel, we extracted SVR model coefficients (‘weights’) to understand which spectral features were used to decode behavioral metrics. As linear SVR estimates of behavioral metrics (Ybehavior) are a combination of neural weights (Wneural) and power estimate (Xneural) inputs (Ybehavior=WneuralXneural+intercept), positive weights described the association between the presence of a specific frequency band with higher metric output values. Conversely, negative weights described the absence of a neural feature when metric output values were high.

To test whether specific clusters of features (≥3 contiguous spectral features) were consistently weighted across recordings, the distribution of each feature’s SVR model weights (averaged over three adjacent features) across recordings were compared to the distribution of metric-shuffled SVR model weights using a contiguity-sensitive permutation test (Ahn et al., 2020). Over 10,000 iterations, each recording’s SVR weight values were shuffled across the two models (empirical vs. shuffled), and the difference between individual feature distributions across electrodes was assessed using a paired t-test. The empirical T-statistic for each feature distribution (e.g., microelectrode β1-3 empirical v. shuffled) was then compared to the null T-statistic distribution. The probability of the observed difference in T-statistics was calculated empirically from the null T-statistic distribution. This procedure was also used for understanding whether specific features were weighted differently for specific metric decoding models (e.g., microelectrode β1-3 tremor v. slowness).

For datasets collected using the high-density STN survey, SVR models were trained using recordings throughout the STN. Specifically, data from each depth were split into 2:1 train:test sets and aggregated for whole-STN SVR model fitting. r2-Values from SVR models trained on high-density microelectrode data were compared to r2-values from stationary microelectrode data using the following LMM: yr=Xexperimentβ+Zu+ϵ, where yr represented the SVR decoding r2-value and Xexperiment represented the categorical labels of experiment type.

To determine if whole-STN models could decode metrics above chance, r2-value distributions from empiric and metric-shuffled decoding models were compared using the Wilcoxon test. To understand if specific recording sites contained information specific to individual metrics, metrics were estimated at each depth by applying the whole-STN SVR model weights to spectral features recorded at that depth. From there, r2-values for each recording/depth were calculated between estimated and observed metrics.

Anatomical analysis of metric decoding

Request a detailed protocol

To compare whether specific motor features were better decoded in different regions of the STN, tremor and slowness r2-values were plotted in MNI coordinate space. All voxels and their associated r2-values with MER recordings were then compared with a voxel-wise paired t-test with AFNI’s ‘3dttest++’ function. Each resulting voxel had an associated Z-statistic that was generated from 10,000 permutations of shuffling tremor and dataset r2-values across voxels using AFNI’s equitable thresholding and clustering (ETAC) algorithm (Cox, 2019). Briefly, ETAC was used to estimate dataset-specific empirical statistical (e.g., T-statistic) and cluster-size (number of adjacent voxels) thresholds for significant results by running several (10,000) permutations of testing with shuffled data. However, due to the relatively low number of voxels per recording dataset (1 mm3 per recording site compared to whole-brain coverage typically acquired with fMRI scans), no distinct clusters were isolated by the ETAC algorithm. For cortical surface-based comparisons between tremor and slowness r2-values, ‘3dttest++’ was also used, although without the ETAC algorithm as it is not currently implemented for surface-based datasets. Thus, instead of Z-scores computed by ETAC, we examined the T-statistic from ‘3dttest++’ and its associated p-value. Colormaps for cortical figures were obtained from https://github.com/snastase/suma-colormaps ( Nastase, 2018).

Statistical analysis

Request a detailed protocol

Data in text are represented as mean±standard deviation. All statistical tests, unless otherwise specified, were carried out in the ‘scipy’ environment. p-Values were adjusted for multiple comparisons wherever appropriate using the Benjamini-Hochberg procedure with q=0.05 (Benjamini and Hochberg, 1995). Data points (epochs) were aggregated across trials within a single session. When comparing data aggregated across multiple subjects, LMMs were performed using the ‘statsmodels’ toolbox to disentangle the effect of interest (continuous or categorical) from the random effects/unequal contributions of each subject’s dataset (Lindstrom and Bates, 1988; Seabold and Perktold, 2010). In other words, individual epochs/trials/sessions/recordings were considered repeated measurements within individual subjects, and individual subject variability was considered as fixed factors during all statistical comparisons to avoid one subject unduly influencing results. All LMMs were random intercepts models, where each random intercept corresponded to a subject’s dataset, and generally followed the formula of y=Xβ+Zu+ϵ , where y represented the outcome variable, X represented the continuous or categorical predictor variables, β represented the fixed-effect regression coefficients, Z represented the subject-specific random intercepts, u represented the random-effects regression coefficients, and ϵ represented the residuals of the model fit. Once a model was fit, fixed-effect p-values were calculated from Z-scored parameter estimates (fixed-effect coefficients divided by their standard errors) against the normal distribution. When reporting LMM results, categorical comparisons were delineated as ‘A v. B’ whereas continuous regressions were noted as ‘A × B.’ When comparing r2-value distributions with LMMs, descriptive statistics (mean±standard deviation) corresponded to distributions containing only each subject’s highest-performing recording, while LMM statistics (β, Z, p) corresponded to comparisons of all recordings.

Data availability

The raw datasets supporting the current study contain patient information and are unique datasets under continued investigation for additional projects, including those of junior trainees. Deidentified neural/behavioral estimates and related code to reproduce all analyses in the manuscript will be made available in a public repository (Dryad; https://doi.org/10.5061/dryad.h9w0vt4n4). To request raw datasets, please contact the corresponding authors (me@peterlauro.me, wael_asaad@brown.edu) with a project proposal. Based upon the granularity of the data requested and potential for patient information exposure, data sharing would be granted in consultation with the Lifespan IRB. There are no commercial restrictions for these data currently.

The following data sets were generated
    1. Lauro PM
    (2023) Dryad Digital Repository
    Data from: Concurrent decoding of distinct neurophysiological fingerprints of tremor and bradykinesia in Parkinson's disease.
    https://doi.org/10.5061/dryad.h9w0vt4n4

References

    1. Akbar U
    2. Asaad WF
    (2017)
    A comprehensive approach to deep brain stimulation for movement disorders
    Rhode Island Medical Journal 100:30–33.
    1. Parkinson J
    (2002) An essay on the shaking palsy. 1817
    The Journal of Neuropsychiatry and Clinical Neurosciences 14:223–236.
    https://doi.org/10.1176/jnp.14.2.223
    1. Pedregosa F
    2. Varoquaux G
    3. Gramfort A
    4. Michel V
    5. Thirion B
    6. Grisel O
    7. Blondel M
    8. Prettenhofer P
    9. Weiss R
    10. Dubourg V
    11. Vanderplas J
    12. Passos A
    13. Cournapeau D
    14. Brucher M
    15. Perrot M
    16. Duchesnay É
    (2011)
    Scikit-learn: machine learning in python
    Journal of Machine Learning Research 12:2825–2830.
  1. Conference
    1. Shah SA
    2. Tinkhauser G
    3. Chen CC
    4. Little S
    5. Brown P
    (2018) Parkinsonian Tremor Detection from Subthalamic Nucleus Local Field Potentials for Closed-Loop Deep Brain Stimulation
    2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC. pp. 2320–2324.
    https://doi.org/10.1109/EMBC.2018.8512741

Decision letter

  1. Hayriye Cagnan
    Reviewing Editor; University of Oxford, United Kingdom
  2. Tamar R Makin
    Senior Editor; University of Cambridge, United Kingdom
  3. Wolf-Julian Neumann
    Reviewer; Charité - Universitätsmedizin Berlin, Germany

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Concurrent Decoding of Distinct Neurophysiological Fingerprints of Tremor and Bradykinesia in Parkinson's Disease" for consideration by eLife. Your article has been reviewed by 3 peer , one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Tamar Makin as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Wolf-Julian Neumann (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

Reviewers all agree that the work is interesting and beneficial, however, there are several aspects that the authors should address:

1) Clarifying how behavioural measures (tremor, slowness) were derived, and how these relate to clinical scores. The authors should also support their analysis choices such as normalisation of cursor speed and derivation of effective motor control;

2) Analysing sub-thalamic and cortical recordings together to explore connectivity and coupling measures;

3) Clarifying how multiple recordings from each subject have been dealt with for decoding and using R2 as a performance metric.

Reviewer #1 (Recommendations for the authors):

I enjoyed reading this manuscript and think it will be a valuable contribution to a number of research fields following a revision.

1) The manuscript would significantly improve if the authors provided additional information regarding their methods and reorganised their results.

A) How did the authors compute tremor amplitude probability densities in age-matched controls and patients who did not exhibit tremors? Panels 1B and C linked to tremor analysis are difficult to read – did most patients and controls have no tremor and therefore the density functions are decaying from 0 or is there a low amplitude peak?

B) What was the motivation for normalising cursor speed to its minimum and maximum and what are the implications of this normalisation when comparing speed within and across participants?

C) I am not entirely sure how useful age-matched control behaviours are to understanding the main results of the paper – the authors could consider removing these to streamline the Results section.

D) Could authors further clarify the analysis using FWHM to delineate periods of time where metrics were sustained above control levels?

E) the authors start referring to effective motor control in figure 2 but the description appears later in the paper; re-organising figures 2 and 3 would improve readability.

2) If ECOG and STN recordings were acquired simultaneously, how did signals in both structures co-vary? Would considering envelope – envelope, phase – phase, envelope – phase information improve decoding beyond what can be achieved from a single recording site?

3) In the discussion, the authors state that their model supports "tremor related oscillations originating in the STN and propagating to cortex" – what is the evidence for this in the manuscript?

4) Could the authors further discuss how their full spectrum decoder may be implemented in the future for DBS control taking into account device and real-time processing constraints?

5) Could the authors further discuss how tremor/slowness/effective movement decoding from micro-electrode recordings reflect overall activity levels of units (in particular those linked to higher frequencies (γ high and hfo)) (Figures 3 and 4)?

Reviewer #2 (Recommendations for the authors):

Thank you for inviting me to review this interesting study. This is an amazing paper, that I read with enthusiasm. I am not sure why the authors have chosen to neglect all brain signal decoding papers in the field of deep brain stimulation, perhaps they were afraid that this would diminish the novelty. I personally think that the paper and results are sufficiently novel and the paper would have further gained from a thorough discussion of the current field (for review from our group see Merk et al., Exp Neurol 2022; https://doi.org/10.1016/j.expneurol.2022.113993 ; I do not aim to get this paper cited, just want to provide some inspiration).

Reviewer #3 (Recommendations for the authors):

There are a few questions and suggestions that would strengthen the overall conclusions of the manuscript.

The approach they use to obtain high-density recordings of the STN involves first driving the microelectrode to the bottom of the STN, and then in an automated fashion and based on the length of the recording track, incrementally and automatically moving the microelectrode dorsally, interleaving experimental sessions at each increment, until the top of the STN is reached. This is a nice approach for mapping the entire STN. However, there are two questions. First, as far as I understand, the microelectrodes used will have a larger macro electrode contact 3mm dorsal to the microelectrode tip. This means that any tissue that lies 3mm to the most ventral aspect of the recording will be damaged by the larger macro contact. If the recording span of the STN is larger than 3mm, then some microelectrode recordings will be in this damaged region. How do the authors account for this, and would it make sense to discard data that is obtained >3mm dorsal to the dental STN border? Second, different recording tracks may have different lengths and spans. Which recording track is used for this process? As they note, not all the electrodes, therefore, recorded data from the STN, and so the question is whether it would make sense to discard these non-STN recordings.

Tremor amplitude is quantified as the magnitude of the 3-10 Hz filtered signal. How is this converted to a tremor score that is then used for effective motor control? In addition, effective motor control appears to be simply the average of the (effective) tremor and slowness scores. However, the simple average may be misleading as a very good score in one domain could potentially compensate for a poor score in the other. Had the authors instead considered using an effective motor control score that is the minimum of either the tremor or slowness effective scores (1 – their value)? On a related note, the authors convert cursor speed to slowness by normalizing within the session. But should they instead normalize within subjects?

Figure 2 could be presented more clearly. For example, there are no scale bars in 2A, and panels B, C, and D are all different sizes for example.

For figure 3, the authors use decoding models to compare the decoding of tremor to the decoding of slowness to determine which spectral features can distinguish between the two categories and between each category and effective motor control. This is a nice analysis. Interestingly, there is a difference in these features between micro and macro recordings, and the macro electrode features appear somewhat similar to the motor cortex ECoG recordings. Is there a reason this should differ from the micro features? More interestingly, is there a possible link between the macro STN recordings and the ECoG recordings? Have the authors investigated any measures of coupling or connectivity between the two regions?

There are multiple experiments performed at different depths throughout the STN in each subject. It is not clear to me, and I am sure the authors have addressed this point, but when constructing the LLMs, there is included a factor for the subject. However, are all experiments or trials recorded within the same subject considered as independent samples? Could one subject's data be driving these significant regression coefficients?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Concurrent Decoding of Distinct Neurophysiological Fingerprints of Tremor and Bradykinesia in Parkinson's Disease" for further consideration by eLife. Your revised article has been evaluated by Tamar Makin (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1) Could the authors please (a) clarify patients' and control subjects' limb position/posture during the behavioural task; (b) indicate that Louis et al. 2001 observed a relation between Parkinson's rest and action tremor when the UPDRS rest tremor sub-scores and the Washington Heights-Inwood Genetic Study of Essential Tremor Rating Scale were correlated; and (c) discuss why in this study there is a deviation between UPDRS rest tremor sub-score and tremor severity on task vs UPDRS kinetic tremor sub-score despite the previously reported relationship between the two (Louis et al. 2001).

https://doi.org/10.7554/eLife.84135.sa1

Author response

Essential revisions:

Reviewers all agree that the work is interesting and beneficial, however, there are several aspects that the authors should address:

1) Clarifying how behavioural measures (tremor, slowness) were derived, and how these relate to clinical scores. The authors should also support their analysis choices such as normalisation of cursor speed and derivation of effective motor control;

We have sought to clarify the derivation of our motor metrics, and their relationship to UPDRS clinical assessments. In particular, we discussed the relationship of our tremor measurement with resting v. action/postural tremor – while our tremor measurement positively correlates with resting tremor, it trended towards positive correlation with action tremor as well. Because resting tremor is typically thought to be more specific to Parkinson’s disease (and our task uses continuous motion instead of static postures), we felt our task-based tremor measure was appropriate and relevant to PD.

Within-session speed normalization was used primarily to control for differing target speeds across different task trials and versions, while also accounting for variability in patient speed. To ensure this decision did not significantly alter our results, we now have confirmed that decoding slowness using within-subject or no speed normalization produced the same ultimate results (see details later in this response).

In addition, we have described the derivation of effective motor control earlier in the manuscript. When comparing our metric-average approach vs. the reviewer’s metric-minimum approach we saw equivalent decoding performance. We thus continue to report the same metric-averaged results as before.

2) Analysing sub-thalamic and cortical recordings together to explore connectivity and coupling measures;

We apologize for not investigating the (in retrospect) obvious question of combining subthalamic and cortical recordings for metric decoding. We quantified connectivity/coupling between electrodes using magnitude-squared coherence, and used the resulting spectra as neural features for decoding. Perhaps surprisingly, local decoding models out-performed paired-recording decoding models. However, when analyzing the relevant coherence decoding features, macro-STN and cortical coherence-based decoding relied on shared frequency bands (lines 402–418).

3) Clarifying how multiple recordings from each subject have been dealt with for decoding and using R2 as a performance metric.

We have clarified that for datasets compared across subjects using linear mixed models, individual epochs/trials/sessions/recordings are considered as repeated measurements for individual subjects, and individual subject variability was considered during all statistical comparisons to avoid one subject unduly influencing results. We describe our statistical approach further in the “Linear Mixed Model Design” sub-section of Materials and methods (including individual model/analysis design) in lines 646–690, and provide a brief overview of our approach in the “Statistical Analysis” sub-section of Materials and methods (lines 819–822).

All decoding-based results have been converted to r2, and the relevant statistics have been updated with our overall findings remaining intact.

Reviewer #1 (Recommendations for the authors):

I enjoyed reading this manuscript and think it will be a valuable contribution to a number of research fields following a revision.

1) The manuscript would significantly improve if the authors provided additional information regarding their methods and reorganised their results.

A) How did the authors compute tremor amplitude probability densities in age-matched controls and patients who did not exhibit tremors? Panels 1B and C linked to tremor analysis are difficult to read – did most patients and controls have no tremor and therefore the density functions are decaying from 0 or is there a low amplitude peak?

Probability density plots were calculated in a similar fashion for all age-matched controls and patients regardless of the presence of tremor. Specifically, a normalized histogram of tremor amplitudes across all subjects and epochs was calculated. As patients with tremor were often not continuously exhibiting tremor throughout the task, the low-amplitude peaks for patient and control distributions likely represents a noise baseline within the 3–10 Hz band. The right-sided tail in the patient distribution, however, represents the wide range of tremor amplitude expressed by patients.

B) What was the motivation for normalising cursor speed to its minimum and maximum and what are the implications of this normalisation when comparing speed within and across participants?

Because the task’s moving target has a certain speed, this produces a sort of normalization in that the range of movement speed is at least somewhat related to that target speed in most instances. The practical effect of this is that normalization does not introduce a major change, and normalized speed produces the same ultimate results. When comparing speed across subjects and populations, speed was not normalized in order to understand if our task measurements of speed corresponded to bradykinesia in patients when compared to controls. Within individual subjects, speed was normalized in order to better identify when patients were slowing down relative to their own baseline.

Nevertheless we compared decoding r2 values across session-normalized and raw speed values, and found no difference in microelectrode recordings (session-normalized slowness vs. raw speed; r2 = 0.125 ± 0.108 v. 0.113 ± 0.094, LMM β = -0.004, Z = -0.741, p = 0.459), macroelectrode recordings (session-normalized slowness vs. raw speed; r2 = 0.198 ± 0.147 v. 0.181 ± 0.148, LMM β = 0.014, Z = 1.879, p = 0.060), and ECoG recordings (session-normalized slowness vs. raw speed; r2 = 0.314 ± 0.143 v. 0.317 ± 0.144, LMM β = -0.003, Z = -0.523, p = 0.601).

C) I am not entirely sure how useful age-matched control behaviours are to understanding the main results of the paper – the authors could consider removing these to streamline the Results section.

We included age-matched control subjects in the study to help validate our motor metrics as reflecting motor symptoms of PD (i.e. control behavior was used to establish cut-offs for pathologic behavior), across both fixed- and random-pattern versions of our task. In addition, the control data allowed us to characterize the distinct timescales of symptomatic tremor and slowness expression.

D) Could authors further clarify the analysis using FWHM to delineate periods of time where metrics were sustained above control levels?

We now realize this may have been unclear. As described in the Materials and methods section, to determine the timescale of metric fluctuation, autocorrelograms were calculated across each PD subject’s behavioral data using 100 ms epochs. The average full-width half-maximum (FWHM) of the autocorrelograms were considered the minimum time necessary to label motor metric data as a “symptomatic" period. Tremor or slowness were considered “symptomatic" if they exceeded the 95th percentile of aggregate control data, and sustained symptomatic periods were defined as those persisting beyond the population metric FWHM continuously. For effective motor control, epochs were labeled “symptomatic” if they were above the median of the PD subject’s session distribution.

E) the authors start referring to effective motor control in figure 2 but the description appears later in the paper; re-organising figures 2 and 3 would improve readability.

We apologize for the lack of clarity, in which effective motor control was introduced in Figure 2 before being introduced in the text. We have moved the text from the “Effective motor control had characteristic neural signatures” Results sub-section to the “Tremor and slowness were distinct and opposing symptomatic states” (lines 150–156).

2) If ECOG and STN recordings were acquired simultaneously, how did signals in both structures co-vary? Would considering envelope – envelope, phase – phase, envelope – phase information improve decoding beyond what can be achieved from a single recording site?

We thank the reviewer for suggesting this analysis. We analyzed synchrony between cortical and STN recordings by calculating magnitude-squared coherence between these timeseries within 7-second epochs. We then averaged the resulting coherence spectra into the same frequency bins as spectral power, and used these features for metric decoding. Perhaps surprisingly, local decoding models outperformed coherence decoding models for all metrics (lines 402–418). Nevertheless, we were able to observe shared neural features across ECoG and macro-STN recordings.

3) In the discussion, the authors state that their model supports "tremor related oscillations originating in the STN and propagating to cortex" – what is the evidence for this in the manuscript?

We interpreted the co-occurrence of tremor-frequency (4-8 Hz; θ) activity in both STN and cortical recordings as being consistent with tremor-frequency synchronization during tremor episodes, as described in work from our group (Lauro et al., 2021) and others (Hirschmann et al., 2013). We have since clarified the text in the discussion (lines 483–485).

4) Could the authors further discuss how their full spectrum decoder may be implemented in the future for DBS control taking into account device and real-time processing constraints?

Given the translational implications of our results, we appreciate the reviewer for the opportunity to discuss future considerations for implementing our full-spectrum decoder for future prospective experiments. We have added a paragraph to the discussion describing these considerations alongside other studies/techniques for decoding motor dysfunction in PD (lines 522–533).

Reviewer #2 (Recommendations for the authors):

Thank you for inviting me to review this interesting study. This is an amazing paper, that I read with enthusiasm. I am not sure why the authors have chosen to neglect all brain signal decoding papers in the field of deep brain stimulation, perhaps they were afraid that this would diminish the novelty. I personally think that the paper and results are sufficiently novel and the paper would have further gained from a thorough discussion of the current field (for review from our group see Merk et al., Exp Neurol 2022; https://doi.org/10.1016/j.expneurol.2022.113993 ; I do not aim to get this paper cited, just want to provide some inspiration).

Reviewer #3 (Recommendations for the authors):

There are a few questions and suggestions that would strengthen the overall conclusions of the manuscript.

The approach they use to obtain high-density recordings of the STN involves first driving the microelectrode to the bottom of the STN, and then in an automated fashion and based on the length of the recording track, incrementally and automatically moving the microelectrode dorsally, interleaving experimental sessions at each increment, until the top of the STN is reached. This is a nice approach for mapping the entire STN. However, there are two questions. First, as far as I understand, the microelectrodes used will have a larger macro electrode contact 3mm dorsal to the microelectrode tip. This means that any tissue that lies 3mm to the most ventral aspect of the recording will be damaged by the larger macro contact. If the recording span of the STN is larger than 3mm, then some microelectrode recordings will be in this damaged region. How do the authors account for this, and would it make sense to discard data that is obtained >3mm dorsal to the dental STN border? Second, different recording tracks may have different lengths and spans. Which recording track is used for this process? As they note, not all the electrodes, therefore, recorded data from the STN, and so the question is whether it would make sense to discard these non-STN recordings.

We agree that, in our experience, recordings “on the way up” where the microelectrode traverses a region already traversed by macroelectrodes, recording quality can be slightly diminished. However, it is still possible to observe single units which suggests that, to a useful extent, surrounding tissue may “collapse” around the smaller probe as it moves upward. Nonetheless, it is true that this tissue may be altered somewhat due to the earlier traversal, and so we now make mention of this as a potential limitation (lines 543–549).

Given the relative scarcity of data in the intraoperative setting, all recording tracks were considered. In the case of high-density STN survey data, we based drive movement on the “best” trajectory (i.e. longest span of representative STN neurophysiology and somatotopic responses). As you point out however, this means that some recordings may have been collected outside the STN. However, despite this, fixed-location recordings were able to decode well above chance, with macro-STN recordings (which would be more likely to be outside STN) demonstrating superior performance. Because of these shared limitations, we decided to include all high-density STN recordings.

Tremor amplitude is quantified as the magnitude of the 3-10 Hz filtered signal. How is this converted to a tremor score that is then used for effective motor control? In addition, effective motor control appears to be simply the average of the (effective) tremor and slowness scores. However, the simple average may be misleading as a very good score in one domain could potentially compensate for a poor score in the other. Had the authors instead considered using an effective motor control score that is the minimum of either the tremor or slowness effective scores (1 – their value)? On a related note, the authors convert cursor speed to slowness by normalizing within the session. But should they instead normalize within subjects?

Tremor amplitude was min-max normalized within each subjects from 0-1, which was subsequently used in the effective motor control equation.

We thank the reviewer for raising this important point of potentially obscuring motor dysfunction by the use of motor metric averaging. We re-calculated EMC using the method described by the reviewer, which was shown to have lower r2 values across micro-electrode recordings (average vs. minimum EMC; r2 = 0.140 ± 0.014 v. 0.092 ± 0.069, LMM β = -0.026, Z = -4.650, p = 3.32 * 10-6), but was equivocal for macro-electrode recordings (average vs. minimum EMC; r2 = 0.204 ± 0.097 v. 0.170 ± 0.116, LMM β = -0.014, Z = -1.837, p = 0.066) and ECoG recordings (average vs. minimum EMC; r2 = 0.233 ± 0.105 v. 0.245 ± 0.099, LMM β = 0.012, Z = 0.321, p = 0.748).

When comparing session- and subject-normalized slowness measurements, there was no difference in r2 values for micro-electrode recordings (session- v. subject-normalized slowness; r2 = 0.125 ± 0.108 v. 0.126 ± 0.105, LMM β = 0.001, Z = 0.195, p = 0.845), macro-electrode recordings (session- v. subject-normalized slowness; r2 = 0.198 ± 0.147 v. 0.197 ± 0.151, LMM β = 0.000, Z = 0.043, p = 0.966) and ECoG recordings (session- v. subject-normalized slowness; r2 = 0.139 ± 0.135 v. 0.140 ± 0.135, LMM β = -0.001, Z = -0.093, p = 0.926).

Because overall these metric modifications were mostly equivocal, we will continue reporting results with our original metrics.

Figure 2 could be presented more clearly. For example, there are no scale bars in 2A, and panels B, C, and D are all different sizes for example.

We have since updated all subpanels within Figure 2, and hope the results are more clear.

For figure 3, the authors use decoding models to compare the decoding of tremor to the decoding of slowness to determine which spectral features can distinguish between the two categories and between each category and effective motor control. This is a nice analysis. Interestingly, there is a difference in these features between micro and macro recordings, and the macro electrode features appear somewhat similar to the motor cortex ECoG recordings. Is there a reason this should differ from the micro features? More interestingly, is there a possible link between the macro STN recordings and the ECoG recordings? Have the authors investigated any measures of coupling or connectivity between the two regions?

We agree with the reviewer that micro-electrode and macro-electrode STN decoding models revealed intriguing differences in relevant neural features. We think the micro v. macro difference is due in part to the former recording type being less sensitive to local field potentials than the latter (Marmor et al., 2017). As for the similarities between macro-STN and ECoG recordings, this may be due to synchrony or volume conduction. As we re-referenced macro-STN and ECoG recordings in order to reduce volume conduction, shared neural features more likely represent functional coupling or shared representations of a motor state. It should be noted that the ranges of relevant frequencies in the β and γ bands differed between macro-STN and ECoG models, suggesting that there may be specific local features (particularly in sensorimotor cortex) in addition to shared motifs (Dubey and Ray, 2019).

However, to address whether these features represented metric-specific coupling between the STN and cortex, we calculated connectivity via magnitude squared coherence between these timeseries within 7-second epochs. We then averaged the resulting coherence spectra into the same frequency bins as spectral power, and used these features for metric decoding. Perhaps surprisingly, local decoding models outperformed coherence decoding models for all metrics (lines 402–418). Nevertheless, we were able to observe shared neural features across ECoG and macro-STN recordings.

There are multiple experiments performed at different depths throughout the STN in each subject. It is not clear to me, and I am sure the authors have addressed this point, but when constructing the LLMs, there is included a factor for the subject. However, are all experiments or trials recorded within the same subject considered as independent samples? Could one subject's data be driving these significant regression coefficients?

When comparing decoding accuracies across subjects with linear mixed models, individual decoding models from one session-recording (and their associated r or r2 values) are treated as one sample. As you mentioned however, subject is treated as a fixed factor to prevent one subject from skewing overall results. We describe our statistical approach further in the “Linear Mixed Model Design” sub-section of Materials and methods (including individual model/analysis design) in lines 646–690, and provide a brief overview of our approach in the “Statistical Analysis” sub-section of Materials and methods (lines 819–822).

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1) Could the authors please (a) clarify patients' and control subjects' limb position/posture during the behavioural task; (b) indicate that Louis et al. 2001 observed a relation between Parkinson's rest and action tremor when the UPDRS rest tremor sub-scores and the Washington Heights-Inwood Genetic Study of Essential Tremor Rating Scale were correlated; and (c) discuss why in this study there is a deviation between UPDRS rest tremor sub-score and tremor severity on task vs UPDRS kinetic tremor sub-score despite the previously reported relationship between the two (Louis et al. 2001).

We apologize for not describing the subject/limb positioning previously, this had been eliminated to fit formatting requirements. We now describe these details in the “Behavioral Task and Metrics” sub-section of “Materials and methods” on lines 644–656:

“As per standard surgical procedure, patients were positioned in a reclined “lawn-chair” position, supine on the operating table to maintain comfort while allowing patients to engage with the task and clinical assessment. For the task, a boom-mounted display was adjusted to the patient’s line of sight, and patients were asked to verbally confirm their ability to see onscreen task objects. Patients using a joystick had it placed in their lap, while patients using a tablet had it placed on a stand in their lap (angle adjusted to comfort) and a stylus placed in their dominant hand. Healthy control subjects performed the task by sitting in a chair at a table, with the manipulandum secured on the table. Similar to the patients, the adjustable arm-mounted task display was adjusted to patient’s line of sight. Subjects were instructed to follow a green target circle that moved smoothly around the screen by manipulating the joystick or stylus with their dominant hand with the goal of keeping the white cursor within the circle (Figure 1A). All subjects were instructed to not rest their dominant hand on their lap/the table while performing the task.”

We now acknowledge the rating scales used to correlate resting and action tremor in the Louis et al. 2001 in text (lines 118–122):

“Although postural tremor can correlate with resting tremor when patients with PD are measured by the UPDRS for the former and the Washington Heights-Inwood Genetic Study of Essential Tremor (WHIGET) Rating Scale for the latter, resting tremor is thought to be more specific to the PD pathophysiology (Louis et al., 2001).”

Finally, there are several reasons to expect a complex relationship between UPDRS tremor scores (resting and kinetic) and performance on the target tracking task. First, we do observe (and now report) a relationship between resting and action/postural tremor UPDRS scores within our patient sample (Spearman ρ = 0.602, p = 0.002) in the text (lines 115–118), consistent with prior observations. Second, our task entails patient positioning that is different from that which is typically used to assess resting, postural, and action tremors, and is undertaken in the operating room environment under circumstances distinct from the standard in-clinic UPDRS evaluations. Third, there are likely aspects of the movements in our task that make it less directly comparable to resting, postural, or movement tremors at least as those are typically assessed. Importantly, we do observe a relationship between movement velocity and tremor on our task (Figure 2B, left), such that tremor is more strongly expressed at lower movement velocities, as one would expect for PD. Ultimately, the paramount consideration in our experiment was that neural activity should be directly compared against ground-truth behavioral metrics at identical, short timescales. One might argue that because our task did not fully recapitulate UPDRS findings in the clinic, the degree and type of tremor we observed during our experiment in the operating room may not fully reflect an individual’s level of motor dysfunction in a natural environment. Conversely, one might argue that because our task was naturalistic and required fine-motor goal-directed activity (in a way that is different from key aspects of the UPDRS exam), that we might have leveraged a behavior that is perhaps in some way more relevant to everyday function. In either case, it is highly unlikely our measure of tremor is irrelevant to the ground truth of an individual’s disease state. Most important, however these questions might be resolved, they do not diminish the value of an immediate, objective measure of tremor in order to understand its neurophysiological basis.

https://doi.org/10.7554/eLife.84135.sa2

Article and author information

Author details

  1. Peter M Lauro

    1. Department of Neuroscience, Brown University, Providence, United States
    2. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, United States
    3. The Warren Alpert Medical School, Brown University, Providence, United States
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing
    For correspondence
    me@peterlauro.me
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work)
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8569-6427
  2. Shane Lee

    1. Department of Neuroscience, Brown University, Providence, United States
    2. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, United States
    3. Norman Prince Neurosciences Institute, Rhode Island Hospital, Providence, United States
    4. Department of Neurosurgery, Rhode Island Hospital, Providence, United States
    Contribution
    Conceptualization, Supervision, Investigation, Methodology, Writing – review and editing
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work)
  3. Daniel E Amaya

    1. Department of Neuroscience, Brown University, Providence, United States
    2. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, United States
    Contribution
    Investigation, Methodology, Writing – review and editing
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work)
  4. David D Liu

    Department of Neurosurgery, Brigham and Women’s Hospital, Boston, United States
    Contribution
    Methodology, Writing – review and editing
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work)
  5. Umer Akbar

    1. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, United States
    2. The Warren Alpert Medical School, Brown University, Providence, United States
    3. Norman Prince Neurosciences Institute, Rhode Island Hospital, Providence, United States
    4. Department of Neurology, Rhode Island Hospital, Providence, United States
    Contribution
    Supervision, Investigation, Writing – review and editing
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work)
  6. Wael F Asaad

    1. Department of Neuroscience, Brown University, Providence, United States
    2. Robert J. and Nancy D. Carney Institute for Brain Science, Brown University, Providence, United States
    3. The Warren Alpert Medical School, Brown University, Providence, United States
    4. Norman Prince Neurosciences Institute, Rhode Island Hospital, Providence, United States
    5. Department of Neurosurgery, Rhode Island Hospital, Providence, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Writing – review and editing
    For correspondence
    Wael_Asaad@brown.edu
    Competing interests
    The authors have patents and patent applications (US patent #: 17312155, 17470710) broadly relevant to Parkinson's disease (but not directly based upon this work). WFA has received proprietary equipment and technical support for unrelated research through the Medtronic external research program
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4406-9096

Funding

National Institute of Neurological Disorders and Stroke (T32MH020068)

  • Peter M Lauro

Doris Duke Charitable Foundation (Clinical Scientist Development Award#2014101)

  • Wael F Asaad

National Institute of General Medical Sciences (P20 GM103645)

  • Wael F Asaad

Neurosurgery Research and Education Foundation

  • Wael F Asaad

Lifespan Norman Prince Neurosciences Institute

  • Shane Lee

Brown University Robert J. and Nancy D. Carney Institute for Brain Science

  • Peter M Lauro

NIH Office of the Director (S10OD025181)

  • Wael F Asaad

Medtronic

  • Wael F Asaad

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful for the generous participation of our patients in this study. We thank Kelsea Laubenstein-Parker for technical assistance, Karina Bertsch for administrative support, and Ann Duggan-Winkle for clinical support. We also thank James Yu, Minkyu Ahn, David Segar, Tina Sankhla, and Daniel Shiebler for helping develop the motor task experiment. Finally, we thank Menem Andria and Alpha Omega for adding computer-controlled microdrive capabilities. National Institutes of Health Training Grant NINDS T32MH020068 (PML), Doris Duke Clinical Scientist Development Award #2014101 (WFA), National Institutes of Health COBRE Award: NIGMS P20 GM103645 (PI: Jerome Sanes) supporting WFA, Neurosurgery Research and Education Foundation grant (WFA), Lifespan Norman Prince Neurosciences Institute Brown University Robert J and Nancy D Carney Institute for Brain Science. Part of this research was conducted using computational resources and services at the Center for Computation and Visualization at Brown University, with funding provided by an NIH Office of the Director grant S10OD025181. WFA has received proprietary equipment and technical support for unrelated research through the Medtronic external research program.

Ethics

Human subjects: Patients and control subjects agreeing to participate in this study signed informed consent, and experimental procedures were undertaken in accordance with an approved Rhode Island Hospital human research protocol (Lifespan IRB protocol #263157) and the Declaration of Helsinki.

Senior Editor

  1. Tamar R Makin, University of Cambridge, United Kingdom

Reviewing Editor

  1. Hayriye Cagnan, University of Oxford, United Kingdom

Reviewer

  1. Wolf-Julian Neumann, Charité - Universitätsmedizin Berlin, Germany

Version history

  1. Preprint posted: February 14, 2022 (view preprint)
  2. Received: October 12, 2022
  3. Accepted: May 26, 2023
  4. Accepted Manuscript published: May 30, 2023 (version 1)
  5. Version of Record published: June 13, 2023 (version 2)

Copyright

© 2023, Lauro et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 434
    Page views
  • 67
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Peter M Lauro
  2. Shane Lee
  3. Daniel E Amaya
  4. David D Liu
  5. Umer Akbar
  6. Wael F Asaad
(2023)
Concurrent decoding of distinct neurophysiological fingerprints of tremor and bradykinesia in Parkinson’s disease
eLife 12:e84135.
https://doi.org/10.7554/eLife.84135

Further reading

    1. Neuroscience
    Amanda J González Segarra, Gina Pontes ... Kristin Scott
    Research Article

    Consumption of food and water is tightly regulated by the nervous system to maintain internal nutrient homeostasis. Although generally considered independently, interactions between hunger and thirst drives are important to coordinate competing needs. In Drosophila, four neurons called the interoceptive subesophageal zone neurons (ISNs) respond to intrinsic hunger and thirst signals to oppositely regulate sucrose and water ingestion. Here, we investigate the neural circuit downstream of the ISNs to examine how ingestion is regulated based on internal needs. Utilizing the recently available fly brain connectome, we find that the ISNs synapse with a novel cell-type bilateral T-shaped neuron (BiT) that projects to neuroendocrine centers. In vivo neural manipulations revealed that BiT oppositely regulates sugar and water ingestion. Neuroendocrine cells downstream of ISNs include several peptide-releasing and peptide-sensing neurons, including insulin producing cells (IPCs), crustacean cardioactive peptide (CCAP) neurons, and CCHamide-2 receptor isoform RA (CCHa2R-RA) neurons. These neurons contribute differentially to ingestion of sugar and water, with IPCs and CCAP neurons oppositely regulating sugar and water ingestion, and CCHa2R-RA neurons modulating only water ingestion. Thus, the decision to consume sugar or water occurs via regulation of a broad peptidergic network that integrates internal signals of nutritional state to generate nutrient-specific ingestion.

    1. Neuroscience
    Lucas Y Tian, Timothy L Warren ... Michael S Brainard
    Research Article

    Complex behaviors depend on the coordinated activity of neural ensembles in interconnected brain areas. The behavioral function of such coordination, often measured as co-fluctuations in neural activity across areas, is poorly understood. One hypothesis is that rapidly varying co-fluctuations may be a signature of moment-by-moment task-relevant influences of one area on another. We tested this possibility for error-corrective adaptation of birdsong, a form of motor learning which has been hypothesized to depend on the top-down influence of a higher-order area, LMAN (lateral magnocellular nucleus of the anterior nidopallium), in shaping moment-by-moment output from a primary motor area, RA (robust nucleus of the arcopallium). In paired recordings of LMAN and RA in singing birds, we discovered a neural signature of a top-down influence of LMAN on RA, quantified as an LMAN-leading co-fluctuation in activity between these areas. During learning, this co-fluctuation strengthened in a premotor temporal window linked to the specific movement, sequential context, and acoustic modification associated with learning. Moreover, transient perturbation of LMAN activity specifically within this premotor window caused rapid occlusion of pitch modifications, consistent with LMAN conveying a temporally localized motor-biasing signal. Combined, our results reveal a dynamic top-down influence of LMAN on RA that varies on the rapid timescale of individual movements and is flexibly linked to contexts associated with learning. This finding indicates that inter-area co-fluctuations can be a signature of dynamic top-down influences that support complex behavior and its adaptation.