1. Neuroscience
Download icon

Distinct population code for movement kinematics and changes of ongoing movements in human subthalamic nucleus

  1. Dennis London  Is a corresponding author
  2. Arash Fazl
  3. Kalman Katlowitz
  4. Marisol Soula
  5. Michael H Pourfar
  6. Alon Y Mogilner
  7. Roozbeh Kiani  Is a corresponding author
  1. Center for Neural Science, New York University, United States
  2. Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, United States
  3. Neuroscience Institute, NYU Langone Health, United States
  4. Department of Psychology, New York University, United States
Research Article
  • Cited 0
  • Views 891
  • Annotations
Cite this article as: eLife 2021;10:e64893 doi: 10.7554/eLife.64893

Abstract

The subthalamic nucleus (STN) is theorized to globally suppress movement through connections with downstream basal ganglia structures. Current theories are supported by increased STN activity when subjects withhold an uninitiated action plan, but a critical test of these theories requires studying STN responses when an ongoing action is replaced with an alternative. We perform this test in subjects with Parkinson’s disease using an extended reaching task where the movement trajectory changes mid-action. We show that STN activity decreases during action switches, contrary to prevalent theories. Furthermore, beta oscillations in the STN local field potential, which are associated with movement inhibition, do not show increased power or spiking entrainment during switches. We report an inhomogeneous population neural code in STN, with one sub-population encoding movement kinematics and direction and another encoding unexpected action switches. We suggest an elaborate neural code in STN that contributes to planning actions and changing the plans.

Introduction

The ability to change an ongoing plan of action is essential to success in a dynamic world. How do we switch from executing one action to another? The gating of movements is thought to differentially engage cortico-striatal-thalamic pathways. Executing a planned movement involves activation of the direct pathway (Kravitz et al., 2010; Sippy et al., 2015), while preventing or stopping a movement involves activation of the indirect (Kravitz et al., 2010) and hyperdirect pathways (Aron and Poldrack, 2006; Chen et al., 2020). The subthalamic nucleus (STN), located at the intersection of the indirect and hyperdirect pathways, is theorized to globally suppress movement through excitation of the output structures of the basal ganglia (Frank, 2006). These theories are supported by the experimental observations that withholding a planned but uninitiated action is associated with increased firing rates in STN neurons (Bastin et al., 2014; Isoda and Hikosaka, 2008; Pasquereau and Turner, 2017; Schmidt et al., 2013) and elevated local field potential (LFP) power in the beta frequency range – a hallmark of inhibition effects of STN on action planning circuitry (Alegre et al., 2013; Benis et al., 2014; Kuhn et al., 2004; Ray et al., 2012). A prediction of the current theories is that STN contributes to switching an action plan by increasing its activity to halt the ongoing action.

However, it is possible that STN plays a more elaborate role in governing actions than just cancelling existing plans. This possibility is supported by past studies that showed STN neurons increase their activity during movement initiation (Bastin et al., 2014; Pasquereau and Turner, 2017) and distinct LFP activity occurs during behavioral conflicts (Zavala et al., 2013; Zavala et al., 2014), suggesting that the STN may do more than globally suppress movement.

To investigate STN neural responses during action switches, we developed a novel behavioral task in which human subjects with Parkinson’s disease undergoing implantation of deep brain stimulation (DBS) electrodes moved their hands on a straight or turning trajectory from a home position to target regions on either side of the home position (Reach or Planned Turn trials). On a random subset of trials, intermixed with others, subjects were cued to initiate a straight reach but then cued to switch to a turning trajectory mid-movement (Impromptu Turn trials). We recorded spiking activity of STN units and LFPs from the same sites during the task. Our new task contrasts with prior studies of STN which focused largely on the starting or withholding of simple movements, such as button presses or joystick movements in humans (Alegre et al., 2013; Bastin et al., 2014; Benis et al., 2014; Ray et al., 2012), saccades or reaches in non-human primates (Isoda and Hikosaka, 2008; Pasquereau and Turner, 2017), or nose-pokes in rats (Schmidt et al., 2013). Comparing STN neural responses during Impromptu Turn, Planned Turn, and Reach trials provides a framework to examine current theories on the role of STN in movement planning and execution.

We found that interruption of a planned movement with a cue to switch movement trajectory was associated with decreased firing rates in STN, contrary to predictions of prevailing theories. Careful investigation of the neural responses revealed an elaborate code, in which the population of neurons encoded the movement plan, hand kinematics during the movement, and changes in the movement plan. Unsupervised clustering of the population response dynamics revealed functional clusters of neurons, one specialized for encoding movement kinematics and the other for encoding changes in movement plan. The neurons also exhibited a wide range of response lags, some predicting upcoming changes in kinematics and some reflecting recently executed changes, creating a time reservoir for planning and tracking of actions. Finally, we show that switching the movement trajectory was associated with constant beta power in the LFP and similar entrainment of action potentials with beta frequency oscillations, further supporting the idea that switching ongoing actions has a distinct neural mechanism from withholding uninitiated actions. Together, our results suggest an elaborate population code in STN and an intricate role in planning and execution of complex actions, beyond what is predicted by existing theories.

Results

Diverse responses of STN units

We recorded single and multiunit spiking activity from the STN while human subjects with Parkinson’s disease performed single-step or multistep reaching movements during surgeries to implant DBS electrodes. Subjects controlled a cursor on a screen with free hand movements in 3D space as recorded by a stereoscopic camera (Figure 1A). They used the hand contralateral to the recorded STN to move the cursor into a fixation area. After a variable-duration fixation period, they received an instruction cue to make a horizontal reaching movement to a target region (Reach trials) or a horizontal and then upward reaching movement to another target region (Planned Turn trials). Movement was allowed to begin only after a go signal (disappearance of fixation point), separating movement planning and execution phases of the trial. On a fraction of trials initially cued as Reach trials, the target changed mid-movement, instructing the subject to switch the ongoing horizontal movement to a vertical movement (Impromptu Turn trials) to reach the same target region as in the Planned Turn trials. Subjects completed all trial types with similarly high accuracies and comparable movement onset latencies (Figure 1B,C), with their hand trajectories following the instructed paths (Figure 1D–F). However, the latencies of turns in Impromptu Turn trials were systematically longer than movement onset latencies (Figure 1C) likely due to the overhead of switching an ongoing movement into a new movement trajectory. This latency difference suggests that the computations involved in changing an ongoing plan might be distinct from those for starting a planned movement.

Subjects completed cued single-step and multistep reaching movements during intra-operative recordings from STN.

(A) Subjects completed trials involving planned linear movements (Reach), planned two-step movements (Planned Turn), and trials initially cued as Reach trials with a later cue mid-movement to alter the trajectory (Impromptu Turn). Cursor sequence on the screen indicates movement trajectory. (B, C) Trials that met the accuracy criteria as a fraction of the total trials (B) and movement latency (C) were similar on all trial types. However, the time from the turn cue to the change in movement trajectory on Impromptu Turn trials (turn latency) was larger than the movement latency. (D) Example movement trajectories. Each line illustrates a sample trial. The fixation point is indicated by the most leftward point on each trajectory. Start positions are shifted horizontally and leftward trials are flipped horizontally to improve clarity and facilitate comparison of trials in the same group. (E) Mean movement trajectories across all subjects show successful completion of distinct movement plans across different trial types. Shading is SE. (F) Hierarchical clustering of movement trajectories reveals segregation of Reach trajectories from the Planned and Impromptu Turn trajectories. Ns and Nt indicate the number of subjects and the number of trials, respectively, across all panels.

We recorded from 39 units in 9 STNs of 8 subjects (one subject completed separate recording sessions from STNs of both hemispheres). Subject demographics and the number of units recorded in each subject are shown in Supplementary file 1. The recorded neurons showed a variety of response dynamics (Figure 2, Figure 2—figure supplement 1), suggesting an inhomogeneous neural population in the STN. To demonstrate this response diversity, we calculated the peri-stimulus time histograms (PSTHs) of each unit relative to the fixation, instruction, movement onset, turn onset, and feedback events for each trial type. We analyzed the population dynamics using a principal component analysis on unit PSTHs across the recorded population, focusing on the first four principal components (PCs), which explained 61% of the variance. The population responses occupied distinct locations in this state-space at different times throughout the trial (Figure 2E–H), suggesting encoding of task events in STN. Furthermore, the population response trajectories for different trial types diverged 60 ms after movement onset and remained separated until 610 ms after the feedback (Video 1; Holm–Bonferroni corrected pairwise permutation tests p<0.001), indicating that the STN population distinguishes different movements both during and after their execution.

Figure 2 with 2 supplements see all
Firing rates of STN units represent different task events and movement types.

(A–D) Raster plots (A, B) and corresponding PSTHs (C, D) of two example units. Neural responses are aligned to the fixation (Fix), instruction (Inst), movement onset (Move), turn onset (Turn), and feedback (Fdbk). Each point in (A, B) represents an action potential. Significant differences (p<0.1) between pairs of PSTHs within each unit are indicated by corresponding pairs of colored stars. S and U indicate subject and unit identities, respectively, and Nt indicates the total number of trials per unit. (E–H) Projections of population responses on the top four principal components reveals distinct response patterns associated with different task events and trial types. Ns, Nu, and Nt indicate the number of subjects, units, and trials, respectively. Stars indicate Holm–Bonferroni corrected p<0.001 for pairwise permutation tests. Shadings in CH are SE. Because there was no turn onset on reach trials, we sampled from the distribution of turn onset times on turn trials to create the turn-aligned PSTHs for reach trials (see Materials and methods).

Video 1
State space trajectories for firing rate displayed separately for Reach, Planned Turn, and Impromptu Turn trials.

Shading indicates standard error. PC2 distinguishes different trial types the least and is omitted from these plots for clarity. The first four principal components (including PC2) were included in statistical analysis of principal components data.

Although response patterns varied considerably across units, many units elicited similar motifs across trials. We, therefore, hypothesized that the STN neural population consists of functional clusters with a mixture of selectivities (Figure 3A). To better characterize the population diversity, we performed a hierarchical clustering analysis on the unit responses using the contribution of units to the first four population response PCs (see Methods). The analysis revealed two main classes that together comprised 87% of the units (Figure 3): units of one group modulated their firing rates with ongoing movement kinematics (Movement Units; 41%, present in 6/8 subjects, Figure 3B and Figure 2—figure supplement 1A,C,E), and did not distinguish planned and unplanned turns, whereas the units of the other group had a more attenuated representation of movement kinematics but encoded unexpected changes of the ongoing movement plan, distinguishing different trial types (Turn Units; 46%, present in 6/8 subjects, Figure 3C and Figure 2—figure supplement 1B,D,F). Since we used unsupervised clustering of the PSTH dynamics to identify these groups, units in the same group were not perfectly homogenous and demonstrated diversity that could suggest sub-clusters within each group. To preserve our statistical power, however, we focus only on the two top clusters as a basis for deeper exploration of STN responses in the following sections.

Figure 3 with 1 supplement see all
Two functional groups of neurons with distinct response dynamics.

(A) Normalized responses of all recorded units. Hierarchical clustering of responses revealed two dominant functional groups: Movement Units (41%) and Turn Units (46%). White asterisks indicate single units, which are not clustered together. Example units from each group are shown in Figure 2—figure supplement 1. (B, C) The mean PSTHs across the two functional groups. Horizontal pairs of colored stars indicate times with significantly different firing rates for the corresponding pairs of trial types (p<0.05, Holm–Bonferroni corrected linear mixed effect models). Ns, Nu, and Nt indicate the number of subjects, units, and trials represented.

All subjects had recordings of Movement and/or Turn Units. Our dataset includes both single and multi-units, in which most spikes originate from a single unit (Supplementary file 1). Single and multi-units had similar proportions in the Turn and Movement Unit groups. We did not find distinct distributions of the locations of Turn and Movement Units within STN (Figure 2—figure supplement 2). There was no significant difference in the median baseline firing rates of Turn and Movement Units (36.6 Hz, IQR: [22.2–56.6 Hz] vs. 35.5 Hz IQR: [15.5–88.0], respectively, p=0.82 Wilcoxon rank sum test).

Movement Units respond to kinematic properties of movements

The firing rates of Movement Units were modulated by task events and trial types. They increased at movement onset and decreased after feedback. Firing rates were significantly lower on Planned Turn trials beginning 140 ms after movement onset and lasting until 130 ms before turn onset (Figure 3, p=0.035 and p=0.0070 for Planned Turn vs Reach and Planned Turn vs Impromptu Turn, respectively, Holm-Bonferroni corrected linear mixed effect models with random effects for subjects and units nested within subjects). These changes were also seen when examining individual unit PSTHs in this time interval. Despite limited trial counts for individual units, 38% of Movement Units showed decreased firing on Planned Turn compared to Impromptu Turn trials, significantly above chance (p=5.0 × 10–4, binomial test). However, only 19% of Movement Units showed decreased firing on Planned Turn compared to Reach trials (p=0.21, binomial test). The difference of the firing rates at the population level afforded accurate classification of trials. A decoder trained to classify trials based on population responses achieves a cross-validated accuracy of 67% and 71% for decoding Reach from Planned Turn trials and Planned Turn from Impromptu Turn trials, respectively.

While Movement Units increased their firing rates at the onset of movement, this was not a direct response to the appearance of the go cue itself. As shown in Figure 3—figure supplement 1A, an increase in Movement Unit firing rate occurred ~300 ms after the go cue, a significantly longer delay than would be expected for visual processing of the cue. In contrast, firing rates began to increase ~200 ms before the movement onset and continued to rise until after movement onset, suggesting that the increased activity is associated with movement planning and execution.

Because STN neurons are known to change their activity with the kinematics of movement (Georgopoulos et al., 1983; Tankus et al., 2017), we sought to determine if movement kinematics differed across trial types and if the pre-turn differences of firing rates could be explained by contrasting movement kinematics.

The kinematics of movement varied by trial type, as shown in Figure 4—figure supplement 1. During the initial fixation period, hands were stationary, but after the go signal, horizontal acceleration and horizontal speed increased on all trial types (Figure 4—figure supplement 1A-B). Horizontal acceleration then decreased and became negative, compatible with deceleration necessary to stop the hand at the target location. On Planned and Impromptu Turn trials, vertical acceleration and vertical speed increased at turn onset (Figure 4—figure supplement 1C-D). On Impromptu Turn trials, horizontal speed decreased before the vertical speed increased, and the magnitude of these changes was such that the total speed (calculated as the Euclidean norm of the vertical and horizontal speeds) decreased before turn onset. (Figure 4—figure supplement 1E). To facilitate pooling data across subjects and comparison of movement kinematics with neural responses, the total speed and acceleration were normalized (z-score) within each subject (Figure 4A–D). Mirroring the kinematic profile of subjects’ hands, Movement Units increased their firing rates at movement onset on all trial types and decreased them prior to turn onset on Impromptu Turn trials (Figure 4E).

Figure 4 with 4 supplements see all
Firing rates of Movement Units reflect kinematic parameters of past, ongoing, and future movements.

(A, B) Average movement speed during recordings from Movement and Turn Units. Movement speed was z-scored across trials of each session before averaging across sessions. (C, D) Average movement acceleration. (E–F) Average firing rates of Movement and Turn Units aligned to task events. These panels are identical to those in Figure 3B–C and replicated here to facilitate comparison of firing rate dynamics with changes of movement speed and acceleration. Shading indicates SE in AF. (G–H) The firing rate of each unit was modeled as a linear function of time-lagged speed and acceleration (Equation 12). The distributions of best-fitting time lags are shown across units. Positive lags indicate that firing rates predict future movement kinematics, and negative lags indicate that firing rates reflect past movement kinematics. Insets display the actual and predicted firing rates and the regression coefficients for speed and acceleration. Error bars are 95% confidence intervals.

If the responses of Movement Units reflect upcoming movement kinematics, we would expect their firing rates to be a function of the speed and acceleration that would occur in the future. To test this, we modeled the firing rates of Movement Units as a function of lagged speed and acceleration (Equation 12). We determined the lag of each unit as the time shift that best explained the firing rate changes as a function of speed and acceleration. These kinematic models of firing rate generated predicted PSTHs that closely matched the data (inset of Figure 4G and Figure 4—figure supplement 2A,C,E; in-sample R2 = 0.59, 95% CI: [0.57 0.61], out-of-sample R2 = 0.48, 95% CI: [0.45 0.51] across units) and explained the firing rate dynamics (speed mean β = 0.19, 95% CI: [0.01 0.38]; acceleration mean β = 0.25, 95% CI: [0.07 0.43]). However, the best-fitting models had lags that varied considerably in the population of Movement Units (Figure 4G and Figure 4—figure supplement 2G), ranging from +190 ms to -110 ms. Therefore, some units predicted future movement kinematics (positive lags; 56% of units, mean lag± s.e., 110 ms ± 24 ms), whereas others tracked past movement kinematics (negative lags; 44% of units, mean lag± s.e., –89 ms ± 13 ms). Overall, the population of Movement Units represented the kinematics of future movements as well as those of past and ongoing movements.

For the results above, we used a model that explained neural responses based on total movement speed and acceleration, regardless of movement direction. To test whether acceleration and speed along particular movement directions were more effective at modulating neural responses, we explored a multitude of alternative models in which the firing rate was a function of acceleration and speed along the horizontal and/or vertical directions on the screen. These models were fit separately for each unit and the model with the highest R2 for the held-out data (out-of-sample R2) across all units was chosen (see Methods for description of alternative models and process for holding out data). The best model was identical to the main model in the paragraphs above and contained two predictors: the total speed and total acceleration. The alternative models underperformed in predicting the held-out data. But among them, models that included total speed or total acceleration tended to explain response modulations better than the models that included only directional speed and acceleration (Figure 4—figure supplement 3). Overall, our results suggest that kinematic-dependent response modulations of Movement Units are not strongly dependent on the direction of movement.

Movement Units increase firing rate for medial movements

The models above demonstrate that the modulation of Movement Unit responses with movement speed and acceleration is not directionally tuned. However, that does not necessarily imply that these units lack other forms of directional tuning. In particular, Movement Units have higher mean firing rates for the medial movement direction but similar kinematic-dependent response modulations around the mean for different movement directions. In our task, the instructed movement direction was either medial (towards the recorded STN) or lateral (away from the recorded STN) to the recorded STN. For successful completion of trials, subjects had to select not only the correct movement type (Reach vs. Turn) but also the correct direction (medial vs. lateral). The mean firing rates across trial types from 10 ms before to 270 ms after movement onset were higher on trials where the movement direction was medial, as shown in Figure 5A–D,H (p=0.019 linear mixed effect model). This effect was seen in 8/12 individually examined Movement Units (p=1.7 × 10–7, binomial test), and a population response decoder trained on population responses can separate medial from lateral trials with an accuracy of 80%.

Figure 5 with 2 supplements see all
Movement Units have higher firing rates at the onset of medial movements.

(A–C) PSTHs separated by trial type and medial (lighter shades) vs lateral movements (darker shades). Stars denote times with statistically significant differences between medial and lateral movements. (D) Average PSTHs for medial and lateral movements across all trial types. (E–G) Residual firing rate PSTHs for different trial types after regressing out kinematic effects. (H) Average firing rates in a window starting 10 ms before and ending 270 ms after movement onset. Stars denote p<0.05 (linear mixed effect model). Error bars in all panels are SE.

Since the firing rate of these units is also associated with movement kinematics, we tested whether the mean firing rate difference could be ascribed to different speeds and kinematics on medial compared to lateral movements. However, the speed and acceleration were comparable on medial and lateral trials (Figure 5—figure supplement 1). Furthermore, we used the kinematic model of each unit to predict its firing rate for medial and lateral movements. Residual firing rates around these predictions offer a test for the hypothesis that kinematic differences account for the mean firing rate differences on medial and lateral trials. This hypothesis would predict no systematic discrepancy in the residuals as the model would have already captured the effects of kinematic differences on firing rates. Contrary to this prediction, the model residuals reflected systematic firing rate differences, similar to those explained in the previous paragraph (Figure 5E–H) with the mean residual firing rate across this interval greater for medial trials than lateral trials (p=0.034, linear mixed effect model). A population response decoder is still able to separate trial types based on residual firing rates; its accuracy is 70%. Therefore, Movement Units showed a sizeable difference in firing rates around movement onset of medial and lateral trials. However, their firing rate modulations with movement kinematics were not noticeably directional.

If movement direction affected firing rate, why did models that incorporated the x- and y- components of speed and acceleration perform poorly predicting the firing rate? Those models could accommodate differential modulation of firing rates with the kinematics of medial and lateral movements. For example, if firing rates changed by different amounts with changes of speed on medial compared to lateral movements, models with directional components would outperform. This did not occur in the data; firing rate modulated positively and by the same amounts with the magnitude of speed and acceleration on both medial and lateral trials, but the mean firing rate was overall lower on lateral trials. These differences between medial and lateral trials were only present after movement onset and disappeared prior to Turn Onset. Thus, there appears to be a firing rate offset determined by movement direction that occurs after movement onset.

Turn Units predict whether an ongoing movement plan will change

Movement Unit PSTHs exhibited their greatest modulation at movement onset and at feedback (Figure 3B), which coincided with movement cessation. In contrast, Turn Units changed their firing rates prior to turn onset, as shown in Figure 3C. Between 240 ms before to 80 ms after turn onset, the firing rate on Impromptu Turn trials was significantly smaller than on Planned Turn trials (p=0.0031, Holm–Bonferroni corrected linear mixed effect model). Despite limited trial data, 33% (6/18) of individual Turn Units had significantly smaller firing rates on Impromptu Turn than Planned Turn trials in this interval, significantly more than expected by chance (p=0.0012, binomial test). The difference in firing rates between Planned and Impromptu Turns is 0.2 z-scores, which represents half of the range of modulation of firing rates of Turn Units prior to feedback. This difference is sufficient to decode Planned and Impromptu Turn trials from firing rates; a decoder trained to separate Planned Turn and Impromptu Turn trials based on population responses achieves an accuracy of 74%.

Can this difference be explained by movement kinematics? Like Movement Units, Turn Unit firing rates were also correlated with the lagged kinematics (Figure 4H). The best-fit lags were similarly widely distributed with firing rate lags ranging between 190 ms before to 110 ms after changes in movement kinematics (Figure 4H and Figure 4—figure supplement 2H). However, the strength of this relationship was considerably weaker than for Movement Units (out-of-sample R2, Turn Units 0.13, 95% CI: [0.085 0.18], Movements Units 0.48, 95% CI: [0.45 0.51]; in-sample R2, Turn Units, 0.47, 95% CI: [0.44 0.49], Movement Units 0.59, 95% CI [0.57 0.61]). Furthermore, movement kinematics did not accurately predict the PSTHs of these units (Figure 4—figure supplement 2B,D,F). As with Movement Units, models that included separate predictors for the horizontal and vertical components of velocity and acceleration underperformed as shown in Figure 4—figure supplement 4. We used the residual firing rates to determine whether movement kinematics explained away the modulation of firing rates prior to Turn Onset. As in previous sections, residual firing rates close to zero would imply that firing rate modulations before turn onset were explained by the changes in speed and acceleration. In contrast, the mean residual firing rates of Turn Units across the same interval (240 ms before to 80 ms after Turn Onset) were significantly greater than zero on Planned Turn trials (p=0.028, Holm–Bonferroni-corrected Wilcoxon-signed rank test). Furthermore, the mean residual firing rate over this interval was smaller on Impromptu Turn trials than on Planned Turn trials (p=0.0069 and Holm–Bonferroni corrected linear mixed effect model). Similar results were obtained with a variety of kinematic models. For example, an alternate model with total acceleration as the only predictor, which provided slightly better out-of-sample R2 for the Turn Units alone (Figure 4—figure supplement 4C), also had smaller residual firing rates on Impromptu Turn trials than Planned Turn trials over this interval (p=0.0031, Holm–Bonferroni corrected linear mixed effect model). In fact, using the residual firing rates does not appreciably change the decoding accuracy (72% vs 74% for the full responses), suggesting that movement kinematics minimally contributed to the decoder’s success. Therefore, independent of changes in movement kinematics, Turn Unit firing rates were significantly smaller before turn onset on Impromptu Turn trials compared to Planned Turn trials. Unlike Movement Units, Turn Units did not display any difference when comparing medial and lateral trials (Figure 5—figure supplement 2), and they did not respond to the go or turn cues (Figure 3—figure supplement 1B). Overall, Turn Units provide a clear code for changes in ongoing action plans.

Population activity in STN encoding movement direction and movement type precedes the movement

The response properties of the Movement and Turn Units suggest that the population of STN neurons simultaneously encode the kinematics of movements, as well as movement directions and changes of ongoing motor plans. To gain insight into the dynamics of the representation of movement directions and plan changes, we used our kinematic models to remove response fluctuations related to movement speed and acceleration. We then performed a principal components analysis on the residual firing rates of all units to identify dominant response modulations unrelated to movement kinematics. This allowed us to capture population dynamics hidden in the analyses of mean residual firing rates (Figure 5 and Figure 5—figure supplement 2).

Figure 6A–I shows distinct population activity for medial and lateral trials from 30 ms after the instruction to 80 ms after movement onset. The mean squared distance between trajectories over this time period was higher than expected by chance (Figure 6K, p<0.001, permutation test; also see Video 2 for state space trajectories). Thus, medial and lateral trials were associated with distinct population responses well before movement onset and for tens of milliseconds following the movement onset.

Distinct population response dynamics for different movement directions and trial types.

(A–I) Residual firing rate trajectories in 3D principal component space for different trial types and around different task events. Later times are indicated by darker colors. Pairs of trajectories in each panel correspond to medial and lateral trials. Triangles indicate alignment times: instruction (first column), movement onset (second column), and turn onset (third column). Shading indicates standard error. (J) Mean squared distance of trajectories of pairs of trial types around movement onset relative to that expected by chance. A distance equal to zero indicates that the difference is not distinct from chance. (K) Mean squared distance between medial and lateral trajectories around movement onset relative to that expected by chance. The distance is aggregated across all trial types. ★ and ★★ indicate p<0.05 and p<0.001, respectively (permutation test). Error bars are 95% confidence intervals.

Video 2
State space trajectories for residual firing rate after removing the contribution of kinematics separated by trial type and direction (medial vs lateral).

The population activity also distinguished different movement types prior to movement onset. Planned Turn trials showed distinct activity patterns from Reach trials between 130 ms before until 60 ms after movement onset (Figure 6J; p=0.024, Holm–Bonferroni corrected permutation test). A similar difference was also present between Planned Turn and Impromptu Turn trials (Figure 6J, p=0.002, Holm–Bonferroni corrected permutation test). There was no difference between activity patterns on Reach and Impromptu Turn trials (p=0.087, Holm–Bonferroni corrected permutation test). This is expected since Impromptu Turn trials begin as Reach trials.

Movement and turn units desynchronize from beta band oscillations at movement onset

Neuronal spiking in STN is typically phase-locked to beta oscillations (Kühn et al., 2005; Weinberger et al., 2006), which are themselves shaped by inhibition in the local circuit and related to movement suppression. These response properties are hypothesized to be key to the role of STN for inhibitory control of movement in the cortical-basal ganglia network. Since changes in firing rate of the two subpopulations of STN neurons in our study were associated with changes in movement, we examined how these changes related to spike-field locking. Existing theories predict that locking of action potentials to beta oscillations would decrease prior to the movement onset. Furthermore, they predict that beta spike-field locking might increase prior to the turn onset on Impromptu Turn trials, when subjects decrease their speed and halt horizontal movement before initiating a vertical movement toward the new target location (Figure 4—figure supplement 1).

Both Movement and Turn Units spiked in synchrony with beta oscillations as quantified by the Pairwise-Phase Consistency (PPC), shown in Figure 7. PPC averaged across all trial types and times was significantly different from zero for Movement and Turn Units (median PPC, 0.0042, IQR: [0.0014, 0.010], p=0.0023 for Movement Units, and 0.0056, IQR: [6.6 × 10–4, 0.026], p=0.0033 for Turn Units, Wilcoxon signed-rank tests). This corresponded to a firing rate peak-to-trough modulation of 30% and 35%, respectively, for Movement and Turn Units.

Figure 7 with 1 supplement see all
STN units fire in synchrony with beta oscillations except during the movement phase of the task.

(A–F) Pairwise-phase consistency as a function of time and frequency for Reach (A, B), Planned Turn (C, D), and Impromptu Turn (E, F) trials shows decreased synchronization between spikes and LFP in the beta range around movement onset and turn onset. (G–H) Changes of PPC within the beta band aligned to task events. Shading indicates SE. There was no significant difference between trial types. (I, J) Box plots show the median and interquartile range of the beta band PPC across all trial types in the analysis windows for each task event. Whiskers indicate the range of the data excluding outliers, which are defined as more extreme than 1.5× IQR outside of the IQR (stars indicate PPC significantly different from 0, p<0.05, Holm–Bonferroni corrected Wilcoxon signed-rank test).

However, contrary to past theoretical predictions, PPC was not significantly different between trial types (no significant clusters detected by cluster mass test, PPC averaged across movement and turn epochs not different between trial types, p=0.94, p=0.68, p=0.93, for Movement Units, Turn Units, and the combined neural population, respectively, Friedman’s test). Furthermore, the action potentials of Movement Units were uncoupled from the beta oscillation around the time of the turn onset (median PPC = 5.1 × 10–4, IQR: [−0.0045, 0.0043], Wilcoxon signed-rank test p=0.88). Notably, the uncoupling occurred even earlier for Turn Units, around the time of movement onset (median PPC = 0.0020 IQR: [–8.6 × 10–4, 0.012], p=0.064, Wilcoxon signed-rank test; compare Figure 7I,J). Compatible with these observations, the induced beta power did not differ across trial types (no significant clusters detected by cluster mass test, induced power averaged across movement and turn epochs not different between trial types, p=0.72, Friedman’s test), decreased below baseline shortly before movement onset, and returned to baseline after feedback (Figure 7—figure supplement 1). These negative results were not simple byproducts of lack of statistical power. Based on our sample size and measurement noise, our Friedman’s test could detect median PPC differences across trial types as small as 0.0054 for Movement Units, 8.9 × 10–4 for Turn Units, and 9.2 × 10–4 for the combined neural population (significance p-value threshold set to 0.05 for the power analyses). Also, our signed rank test could detect median action potential couplings with beta oscillations as small as 0.004 for Movement Units and 0.003 for Turn Units, well within the range observed in other task epochs (Figure 7G–J). Furthermore, the detection threshold for the beta power change was 0.056, again well within the plausible changes in experimental observations. Overall, there was an absence of noticeable increase in beta power and the phase locking of action potentials to beta oscillations prior to turn onset on Impromptu Turn trials. In addition, firing rates increased during the movement phases of the trials (Figures 2 and 3, and Figure 2—figure supplement 1), while beta power and beta spike-phase locking decreased, suggesting that spiking of STN neurons switches from beta entrainment to alternative patterns during the movement phases of the trial including around turn onset. These observations further suggest that the underlying mechanisms for changes in ongoing action plans are different from those involved in withholding an uninitiated plan (e.g., stop signal tasks), corroborating our earlier conclusion based on the reduction of firing rates prior to turn onset on Impromptu Turn trials.

Discussion

The subthalamic nucleus’ location at the junction of the indirect and hyperdirect basal ganglia pathways and its synchronous activity in Parkinson’s disease have led to the hypothesis that its main role is to inhibit movement (Frank, 2006). Here we use a novel temporally extended reaching task to show that neuronal activity in human STN increases with increasing movement speed and acceleration and decreases prior to changes in an ongoing movement plan. Changes of movement kinematics and motor plans are encoded distinctly in the STN population and represented by different subpopulations within STN. The firing rates of Movement Units increase linearly with speed and acceleration, and while they encode the kinematic changes associated with a change of plan, they do not represent the change of plan itself. In contrast, the firing rates of Turn Units primarily respond to an unpredictable need to switch movement trajectory independent of ongoing movement kinematics.

Our study is the first in humans to examine the differential role of STN neuronal activity during planned and unplanned changes of movements. Our task had three trial types, Reach, Planned Turn, and Impromptu Turn; the latter two consisted of similar two-step trajectories but were distinguished by the presence of an unpredictable cue on Impromptu Turn trials that instructed subjects to switch their ongoing reaching movement to a turning movement. We found that beginning 240 ms before the switch in movement trajectory, population activity was significantly decreased on Impromptu Turn trials compared to Planned Turn trials. Furthermore, we found no changes in induced beta power or phase locking of spike times to beta oscillations during this period.

Our results contrast with STN activity during suppression of uninitiated movement plans, which have commonly been associated with increased spiking activity in the STN population (Bastin et al., 2014; Isoda and Hikosaka, 2008; Pasquereau and Turner, 2017; Schmidt et al., 2013) and increased beta band activity (Bastin et al., 2014; Ray et al., 2012). Prior studies in STN focused on Go/No-Go (Isoda and Hikosaka, 2008; Pasquereau and Turner, 2017) or Stop Signal tasks (Bastin et al., 2014; Pasquereau and Turner, 2017; Ray et al., 2012; Schmidt et al., 2013). In Go/No-Go tasks, subjects choose from one of two actions: movement or withholding of movement. Likewise, in Stop Signal tasks, subjects are cued to make a movement and subsequently cued to withhold it after a short gap from the initial cue. In both cases, successful No-Go or Stop trials involve no movement at all. In contrast, on Impromptu Turn trials in our task, subjects are already executing a movement when they are cued to change it; they then convert their ongoing movement into a new one and change trajectory. We observed a decrease in the firing rate of Turn Units independent of ongoing kinematics in contrast to the increase in firing rate on No-Go Units observed by Isoda and Hikosaka, 2008 and Pasquereau and Turner, 2017. In addition, we did not see an increase in beta power or beta range spike-phase locking around the unplanned turns though beta power was overall decreased after movement onset. It is possible that the decrease in activity we observed prior to changes in ongoing actions is reflective of the STN’s role in interrupting action sequences. However, this would not be explained by existing theories of STN function, which predict an increase in activity. The decreased firing rates may reflect an alternative mechanism or additional role for STN in action control.

In our task, Turn Unit firing rates were highest prior to the turn on trials where subjects had planned to make this movement. They were lowest prior to the turn when subjects had to change their ongoing movement. And they were between these extremes on trials where subjects executed a planned movement but knew that they could be cued to switch movements. An intriguing possibility is that Turn Units may reflect the moment-by-moment confidence or expectation to complete an ongoing or recently planned movement. Consistent with this, at movement onset, the firing rates of Turn Units were at their baseline and not distinct between trials; at this point in the trial, subjects have to make the same movement regardless of cues that may occur in the future. They should thus have the same level of confidence in their initial movement.

While it is possible that our subjects’ underlying Parkinson’s disease influenced our data, this is not sufficient to explain our novel results. Though the STN in Parkinson’s disease contains hypersynchronous beta activity (Brittain and Brown, 2014), beta activity increases during stopping in both Parkinson’s disease (Ray et al., 2012) and obsessive compulsive disorder (OCD) (Bastin et al., 2014). Furthermore, the involvement of beta oscillations and spiking activity in response inhibition has been consistent in patients with Parkinson’s disease (Alegre et al., 2013; Benis et al., 2014; Ray et al., 2012; Zaghloul et al., 2012) and OCD (Bastin et al., 2014) as well as in healthy non-human primates (Isoda and Hikosaka, 2008; Pasquereau and Turner, 2017) and rodents (Schmidt et al., 2013). Our observation that beta activity and beta-locked spiking decrease during movement is also consistent with this framework. However, beta activity and beta-locked spiking did not increase with unplanned movements, suggesting that different mechanisms underlie altering ongoing actions and preventing uninitiated ones. Future studies will determine if this finding can be replicated in non-Parkinsonian subjects.

While Turn Units may encode subjects’ confidence about a selected action, Movement Units seem to more directly encode the kinematics of movement, consistent with studies going back several decades (Georgopoulos et al., 1983). Recently, Tankus et al., 2017 showed that the correlation of firing rates and movement kinematics of uncued movements is lagged by different time periods. We replicated this finding in cued, task-based movements, showing that the firing rate of Movement Units predicts future speed and acceleration or reflects past speed and acceleration. By examining alternative firing rate models, we show that the modulation of Movement Unit activity with movement kinematics is not directionally tuned and is best modeled as functions of total speed and acceleration, quantified as the vector norm of speed or acceleration in the movement plane. In contrast, movement direction is represented in the overall activity of the units and seems to be distinctly coded from speed and acceleration (Figures 5 and 6). Thus, the neural codes for movement direction and kinematics are linearly separable at each moment.

The varying time lags of STN unit responses relative to kinematics could indicate a functional division, in which some neurons are involved in prospective planning of movements, whereas others monitor past movements, implementing a time reservoir of kinematic information for recent movements. The neurons that represent a recent movement could reflect past responses of the STN neurons that lead movements, represent proprioceptive inputs to STN (Abosch et al., 2002; Rodriguez-Oroz et al., 2001), or reflect inputs from the motor cortices that project to STN, either directly (Aron and Poldrack, 2006; Chen et al., 2020) or indirectly through other nodes of the basal ganglia. While proprioceptive feedback cannot shape responses of STN units that encode future movements, cortical inputs could contribute to the activity of those neurons. The same cortical projections could also play a role in shaping the responses of Turn Units.

Prior studies have found differences in neural activity during preparation of movements in different directions (Isoda and Hikosaka, 2008; Zaghloul et al., 2012). Our results complement past studies by showing a similar difference during the movement itself. Examining population activity after regressing out the kinematic effects, we found that medial and lateral movements are associated with distinct population responses, starting 30 ms after the instruction (at least 190 ms before movement onset) and lasting until 80 ms after movement onset (Figure 6). Within the STN population, the mean residual activity of Movement Units after removal of kinematic modulations distinguished medial and lateral movement directions until 270 ms after Movement Onset (Figure 5). In addition to the movement direction, the population activity of different sub-populations of STN neurons distinguished different trial types around key task events (Figures 3 and 6). Thus the population activity in STN encodes diverse properties of the movement plan and is not limited to encoding only the kinematics.

The division of the STN population into two major sub-populations – Movement Units and Turn Units – through unsupervised clustering indicates substantial differences in response properties of the two sub-populations. However, the two sub-populations also share some response properties. For example, Turn Units are moderately responsive to movement kinematics, and the residual firing rates of Movement Units show a trend toward reduction before unplanned turns. Further, the diversity of responses within each sub-population suggests the possibility that each may comprise multiple sub-groups. These sub-groups could be better characterized in a larger dataset, and especially by recording from the same neurons in diverse tasks and movements. Because of partial overlap between response properties of Movement and Turn Units, it is also reasonable to treat the STN population as a mixture of neurons with diverse selectivity. In fact, adopting such a perspective and analyzing the STN population responses was quite effective at revealing that the movement plans for different trial types as well as movement directions are encoded hundreds of milliseconds before movement onset. However, mechanistic understanding of the population response requires accurate characterization of individual neurons in the population and discovering distinct sub-populations within the circuit.

An intriguing hypothesis is that Movement Units encode kinematic information to contribute to execution and possibly planning of actions, whereas Turn Units contribute to maintaining or changing ongoing action plans. Since Movement Unit activity is predictive of future kinematics, this information is likely rapidly transmitted directly from cortex via the glutamatergic hyperdirect pathway. STN neurons that receive these glutamatergic inputs project to GABAergic SNr neurons (Xiao et al., 2015). Compatible with these anatomical connections, stop signals emerge in the STN neural responses earlier than in substantia nigra pars reticulata (SNr) responses (Schmidt et al., 2013). Furthermore, stopping has been linked to activity in the projections from inferior frontal gyrus to STN (Chen et al., 2020; Swann et al., 2009) and low-frequency stimulation of glutamatergic afferents in STN (i.e., the hyperdirect pathway) worsens bradykinesia though stimulation of the motor cortex-STN projection improves it (Gradinaru et al., 2009). We suggest that Movement Units in STN may exert their influence on movement control through this inhibitory pathway. Because the representation of movement kinematics is much weaker in the activity of Turn Units, they likely have different inputs compared to Movement Units, possibly including inputs from the indirect pathway. These insights offer testable hypotheses for future studies on the neural mechanisms that underlie changing ongoing plans of action.

Our results support a versatile and diverse neural code in STN capable of supporting a variety of functions in planning, executing, monitoring, and changing plans of action. As Pasquereau and Turner, 2017 note, this richness in neural code appears incompatible with a narrowly defined role for STN that includes only movement inhibition. We therefore suggest a broader role for STN in motor control. Future studies should explore the rich neural code in STN and clarify its causal role in various aspects of motor control.

Materials and methods

Subjects

Subjects were patients with Parkinson’s disease undergoing implantation of deep brain stimulation electrodes in the subthalamic nucleus (see Supplementary file 1 for demographic information on subjects). They were recruited at their regular appointments, signed informed consent for participation in the study, and completed a training session on the behavioral task at this time. They received monetary compensation for participation in the training session and intra-operative session. All experimental procedures were approved by the Institutional Review Board at NYU Langone Medical Center.

Surgery and electrophysiological recordings

Request a detailed protocol

Subjects did not take their anti-Parkinsonian medications on the day of surgery. Surgical implantation was staged; the DBS lead in each hemisphere and the pulse generator were implanted in separate surgeries. Subjects were given moderate sedation for application of the Leksell stereotactic frame (Elekta, Stockholm, Sweden) and underwent a CT scan for frame localization, which was then aligned with a preoperative structural MRI with gadolinium.

Microelectrodes were targeted to the ventral border of the STN based on structural MRI and the AC-PC coordinates (11 mm lateral to midline, 5 mm inferior to the AC-PC plane, and 4 mm posterior to the AC-PC midpoint). The trajectory to the target was chosen such that the skull entry point was in the vicinity of the coronal suture and modified to avoid structures including vessels, ventricles, and the caudate nucleus. Sedation was typically discontinued after making the incision, and subjects were awake on minimal or no sedation for the electrophysiological recordings. Recordings were performed from a single microelectrode or concurrently from two microelectrodes spaced 2 mm apart at their insertion. For each electrode, an outer cannula was advanced to 15 mm above the target. A microelectrode in an inner cannula (Neuroprobe Sonus Shielded tungsten microelectrode, 1 MΩ impedance, Alpha Omega, Alpharetta, GA) was advanced to the end of the outer cannula and further advanced with a microdrive (Alpha Omega, Alpharetta, GA). When two microelectrodes were used, they were driven concurrently on the same microdrive. The microelectrode tip protruded 3 mm distal to the end of the inner cannula. Recordings were referenced to the outer cannula. Online analysis of spiking activity was performed with the Alpha Omega recording system.

The recorded units spanned the space between the dorsal and ventral borders of STN. Entrance to the STN was noted based on two criteria: (1) real-time alignment of a CT scan performed at the time of surgery with a high-resolution pre-surgical structural MRI, and (2) an increase in background spiking (Novak et al., 2007) and identification of rhythmically firing cells at high firing rates (Hutchison et al., 1998). A decrease in background spiking marked the ventral boundary of STN. In many of our surgeries deeper cells along this trajectory fired rapidly with a relatively constant firing rate, characteristic of SNr cells. This procedure has previously been described in detail by our group (Sterio et al., 2002).

We used the pre-operative high resolution structural MRI with gadolinium, pre-operative high resolution CT scan with the Leksell frame in place, planned trajectory, and intra-operative or post-operative CT scan with DBS lead in place to localize our recordings. Our procedure is as follows. The pre-operatively planned target in AC-PC coordinates was identified on MRI. This MRI was merged with the pre-operative frame-localizing CT scan. Using the frame localization, the target was converted into frame coordinates. The pre-operatively planned trajectory was identified using this target and the ring and arc angle of the Leksell frame. Based on the alignment of the post-operative CT scan, which shows the DBS lead, with the intra-operative CT scan, the recording trajectory was confirmed and corrected whenever necessary (corrections were always <1 mm). The distance of the microelectrode from the ventral STN target was measured using our microdrive. The location of each recording was identified as the point along the trajectory with the measured distance from target. In the case of multiple simultaneously used microelectrodes, the additional trajectories were offset by the distance of electrodes (2 mm) in the appropriate direction. These calculations were performed using the Brainlab Elements stereotaxy software.

Behavioral task

Request a detailed protocol

After identification of a stable single or multiunit, we proceeded with our behavioral task. A monitor on a jointed arm was positioned in front of the patient. A stereoscopic infrared camera (Leap Motion Controller, San Francisco, CA) was positioned under the subject’s contralateral hand with respect to the recorded STN. The subject was instructed to rest his or her elbow on the armrest and make forearm movements to control the cursor. The subject’s hand, therefore, was assumed to move on the surface of a sphere with elbow at its center. The built-in software of the Leap Motion Controller identified the absolute position of the subject’s palm relative to the controller (sampling rate ~100 Hz). The device latency was ~10 ms. The palm position was mapped to a spherical surface with center located at approximately the subject’s elbow and radius equal to the forearm length. The horizontal and vertical angles of this position were scaled to generate the x- and y-coordinates of the cursor on the screen.

A central white circle (fixation point) appeared on the screen at the beginning of each trial. The subject moved the cursor into the fixation point to indicate readiness, beginning the fixation period. After a variable duration (truncated exponential distribution, range, 300–600 ms, mean 400 ms), an instruction cue appeared on the right or left of the screen, which indicated both the type and the direction of the trial. A green circle indicated a Reach trial (60–70% of trials), instructing subjects to move the cursor to a target area defined by the location of the cue. A blue arrow indicated a Planned Turn trial (30–40% of trials), instructing subjects to move toward the cue and then upward to a target area at the top of the screen. A minimum horizontal distance from the fixation point had to be attained before the turn (i.e., subjects could not make a straight, diagonal movement to the target region). Subjects kept their hands still to maintain the cursor in the fixation area until the fixation point disappeared (go cue; truncated exponential delay, range 400–1000 ms, mean 600 ms). They were instructed to initiate their movements as soon as they detected the go cue. Trials on which subjects left the fixation area before the go cue were halted and treated as fixation breaks. On a random half of the trials that were initially cued as Reach trials, the cue changed to a blue arrow during the movement toward the target (turn cue; 43–67% of Reach trials, 30–40% of total trials). The cue changed when subjects entered a rectangular window whose center was chosen from a uniform distribution that spanned 35–65% of distance between the fixation and target points. Subjects had to change their movement paths, making a vertical movement to a target region at the top of the screen. These were designated as Impromptu Turn trials. After completion of each trial, distinct auditory feedback indicated correct responses, incorrect responses (trajectories that did not meet the criteria above), or fixation breaks. The inter-trial interval was 1.5 s.

A session was terminated when the subject showed signs of fatigue or after approximately 10 min. If subjects wished to perform a second block of the task, the microelectrode was moved to a new location with stable units and a new block started. Surgical implantation of DBS electrodes then proceeded per our usual protocol.

Spiking data analysis

Request a detailed protocol

Spike detection and sorting were performed using an offline sorter (Plexon, Dallas, TX) on raw electrophysiological recordings sampled at 44 kHz. Spike-sorting was performed manually using principal components of spike waveforms. Single units, multi-units, and unsorted background spiking (hash) were considered as separate units. The characteristics of all units are shown in Supplementary file 2. Single units were those with the clearest isolation from background and most consistent action potential waveforms (8/39 units). They had a baseline median firing rate of 32.9 Hz (IQR: 18.3–53.9), consistent with single units identified in past studies (Sharott et al., 2014; Weinberger et al., 2006). Furthermore, the median firing rate of units that were not well-isolated was 40.4 Hz (IQR: 20.7–82.6), consistent with the results of Weinberger et al., suggesting that most multi-units are largely contributed to by a single neuron. The baseline firing rate was calculated from spikes in the intervals between trials. Similar conclusions were reached using well-isolated units only.

All analyses were performed using Matlab R2017a. Spikes were aligned to behavioral events. In addition to events defined by the stimuli on the screen, we also identified two events based on hand movements on each trial. Movement onset was defined as the period of maximum acceleration that most closely preceded the peak velocity of cursor. Turn onset was found as the time of maximum angular acceleration of the cursor. The calculated movement and turn onsets were examined visually relative to each movement trajectory and manually corrected as needed by author D.L. Manual corrections of these events, if needed, were small and blinded to the electrophysiological data and the trial type.

To compare spiking activity around the turn cue on Impromptu Turn trials with spiking activity on Reach and Planned Turn trials, we sampled randomly from the distribution of turn cue onsets on Impromptu Turn trials to define corresponding analysis windows on Reach and Planned Turn trials. On Reach trials, subjects did not make a turn movement, and therefore, there was no turn onset event. To compare responses around the onset of the turn movement on Planned and Impromptu Turn trials with Reach trials, we sampled randomly from turn onset times on Planned Turn and Impromptu Turn trials to define corresponding analysis windows on Reach trials.

We sought to determine spiking activity associated with the variety of task events. Inter-event intervals were variable because they were selected from random distributions or depended on the subject’s responses. In calculating PSTHs aligned to an event, we censored spikes that occurred before the preceding event or after the following event. Spikes occurring within 500 ms of the end of the preceding trial were also censored. For smoothing the PSTHs, we used a causal kernel (alpha function, α2τeατ for τ>0 and 0 for τ≤0, with α=20s−1). PSTH error bars were calculated as the standard error over trials. Units were excluded from analyses if there were fewer than four correct trials of each type or if the subject was unable to perform the task (mean ± s.d., 21.8 ± 8.1 trials per trial type per subject). There was a mean ± s.d. of 18.8 ± 11.9 trials per subject split across three conditions that did not meet the accuracy criteria. Adding these trials to the analyses did not critically change our results. The number of correct trials of each type for each unit is shown in Supplementary file 2. For analysis of PSTHs on medial and lateral trials, we only examined units with at least four correct medial and lateral trials of each type (mean ± s.d., 13.5 ± 3.3 trials per trial type per subject). We had planned to collect data from approximately 25 units in this study as this would allow us to detect a difference of 0.5 standard deviations between conditions assuming pairwise comparisons and a Holm-Bonferroni corrected p-value of 0.05 needed for significance. We had expected to record 3–4 units per STN resulting in a sample size of 9 STNs. We recorded from an average of 1.95 units per recording location. 62.5% of subjects (5/8) had multiple recording sites.

The analysis window aligned to each task event was chosen such that at least 2/3 of trials of each subject contributed to the analysis at each point in time. Because fixation, instruction, and feedback events were separated from their neighboring events by large intervals, their analysis windows were wide (fixation, [–500 ms, + 300 ms]; instruction, [–200 ms, + 300 ms]; feedback [–230 ms, + 1000 ms]; windows are with respect to the corresponding event times). These large windows were adopted to maximize our precision for estimation of firing rates, but similar results were obtained with shorter windows too. Since movement and turn onsets had short latencies relative to the Go and Turn cues, analysis windows aligned to these events had to be narrower (movement onset, [–300 ms, + 270 ms]; turn onset, [−270 ms, + 220 ms]).

To explore the similarity of responses across recorded units, we created a response profile for each unit by concatenating the unit’s PSTHs for each trial type in the analysis windows listed in the previous paragraph. To account for different ranges of firing rates across units, the response profiles of each unit were z-scored based on the mean and standard deviations obtained across all trials. Then, we performed PCA on the response profiles across recorded units to find the key patterns that shaped the profiles. As there was no turn onset event on Reach trials, the principal components were determined based on Planned and Impromptu Turn trials only. Reach trial response profiles were then projected onto this space. The top four principal components explained 61% of the total variance on Planned and Impromptu Turn trials, indicating that a small number of patterns captured the diversity of neural responses. The coefficients of the top four components were used to perform a hierarchical clustering analysis on the recorded units using the Chebyshev distance metric and a complete linkage function. The analysis revealed two main clusters that accounted for 87% of units. The remaining units were part of a third cluster with strong responses around planned turns and activity profiles that resembled those of the movement-selective cluster of units. However, we could not extensively explore this third cluster due to low sample size.

Statistical analysis

Request a detailed protocol

Significant differences between firing rates across different trial types were assessed by averaging firing rates within time intervals of interest and making pairwise comparisons using a linear mixed effects model with a fixed effect for the difference of firing rates and random effects for subjects and for units nested within subjects. When reporting pairwise comparisons in cases where there are more than two groups, the p-values were corrected using the Holm–Bonferroni method.

Choosing analysis intervals in time series data is susceptible to Type I error as one can often define a large number of potential analysis intervals within an epoch of interest. It would therefore be ideal to have a systematic approach that guards against such errors while also maintaining adequate temporal resolution to capture transient effects in the data. We developed such an approach (available at https://github.com/dlondon12/InterruptedReachClustMass); copy archived at SWH swh:1:rev:c597f70832fa3bc86e53aaa6030cba643e60f92b ;London, 2021),, for our dataset by adapting the between-subjects cluster mass test (Maris and Oostenveld, 2007). The cluster mass test allows detection of significant differences over contiguous time points while controlling for the family-wise error rate, addressing the multiple comparisons and low temporal resolution problems in conventional tests.

Traditionally, two groups are compared in the cluster mass test. However, since most of our analyses concerned differences among three groups (i.e., the different trial types), we extended this method. The cluster mass test consists of two basic operations: selecting the time clusters and testing their significance. In the original method, two groups are compared at each point in time to generate a t- or z-statistic that is thresholded to select time clusters at a desired false discovery rate. Instead of calculating a t-statistic for two groups, we calculated an F-statistic for three groups at each point in time, and then thresholded it to select our clusters of interest.

In our analyses, we consider units to be random effects and trial types to be fixed effects (analogous to a repeated measures ANOVA). The F-statistic was calculated as follows. For each unit j, consider a firing rate response profile, rj,m(t) , where m represents the trial type. Let us denote the mean response across units with r¯m(t) , across trial types with r¯j(t) , and across units and trial types with r¯(t) . The total sum of squares, SStot , the sum of squares across the M trial types, SSM , and across the the N units, SSN , are:

(1) SStot(t)=m=1Mj=1N(rj,m(t)r¯(t))2SSM(t)=Nm=1M(r¯m(t)r¯(t))2SSN(t)=Mj=1N(r¯j(t)r¯(t))2

Considering each unit as a repeated measure, the residual sum of squares is

(2) SSerror=SStotSSMSSN

and the statistic for the difference across trial types is

(3) FM1,(N1)(M1)(t)=SSM(t)M1SSerror(t)(N1)(M1)

This statistic was thresholded at the 90th percentile of the F-distribution with parameters M-1 and (N1)(M1) . Contiguous blocks of test statistics above threshold were considered to be time clusters. We then calculated a statistic for each cluster, termed the cluster statistic. Among the commonly used statistics, we chose the sum of SSerror(t) within each cluster in this paper, but we also obtained similar results with other relevant statistics including the sum of SSM(t) or of FM1,(N1)(M1) . The null hypothesis is that the M trial types are no different and thus interchangeable. We therefore created our null distribution by permuting the labels on each trial type within each unit, calculating the F-statistic for each permutation to select clusters, and determining the maximum cluster level test statistic. The cluster level statistics for the real data were compared with this null distribution to calculate p-values.

Because our framework is analogous to a repeated measures ANOVA with a single factor (e.g. trial type) per unit, it can also be generalized seamlessly to a repeated measures ANOVA with two factors (e.g. trial type and movement direction). For this extension, consider the additional factor s, with S levels. Our goal is to determine significant time clusters of each main effect, while controlling for interaction effects. For each combination of the two factors, each unit has a different firing rate profile, rj,m,s(t) . For factor m, the sum of squares is

(4) SSM(t)=NSm=1M(r¯m(t)r¯(t))2

For factors m and s, we must also consider the interaction sum of squares:

(5) SSM,S(t)=Nm=1Ms=1S(r¯m,s(t)r¯m(t)r¯s(t)+r¯(t))2

where r¯m,s(t) is the mean across trials for the factors m and s. Analogous expressions can be constructed for each individual factor and each pair of factors. F-statistics are calculated for each main effect:

(6) FM1,(N1)(M1)(t)=SSM(t)M1SSM,N(t)(N1)(M1)FS1,(N1)(S1)(t)=SSS(t)S1SSS,N(t)(N1)(S1)

These F-statistics are thresholded to generate clusters. The cluster-level statistic is the sum of the sum of squares in the denominator of the respective F-statistic calculation. For each tested effect, there is a separate permutation distribution depending on the ‘exchangeable units’ (Anderson, 2001). For testing factor m, the levels of m are permuted within each unit and without permutation of the other factor. The remainder of the calculation continues as for the one factor case, resulting in separate p-values for each main effect, as in a multi-way ANOVA.

We further extended this framework to the case of multidimensional data, as in the case of principal components analysis. Up until this point we have considered the case of N units which represent repeated measures. Principal components analysis of these units results in dimensionality reduction to D dimensions. These dimensions are not repeated measures; they represent orthogonal components of the data. Therefore, instead of considering one-dimensional repeated measures data of the form rj,m,s(t) , we consider D -dimensional data of the form Qm,s(t)=[q1,m,s(t),,qD,m,s(t)] , where Qm,s(t) is the projection of rj,m,s(t) onto the first D principal components, and qd,m,s(t) is the projection onto the d-th principal component.

Our goal is to find the times at which the trajectories in the PC space are significantly different. As there are no repeated measures, we cannot calculate F-statistics, but we can use the sum of squares directly. Again, we analyze each factor separately: for each trial type m, we consider Q¯m(t), where each time point is the centroid across the levels of s, and Q¯(t) , where each time point is the centroid across the levels of m and s. We calculate the sum of squares by adding the squared Euclidean distances from the centroid by analogy with Equation 4:

(7) SSM=Sm=1MQ¯m(t)Q¯(t)2

where . is the Euclidean norm.

We threshold these sums of squares above the 90th percentile as was done with the F-statistics. To determine the value of the 90th percentile, we generate a permutation distribution as follows. For each unit, we randomly permute the trials across the appropriate factor without permutation of the other factors. The mean firing rates are calculated for each unit across the permuted trials, and the population activity patterns projected on the same principal components as the original dataset. These projections are then used to calculate the sums of squares. After creation of the clusters, the cluster-level statistics are calculated by adding the sum of squares values within each cluster. The maximum cluster-level statistic on each permutation is used to generate the distribution of cluster level statistics across permutations. The p-values for clusters of the original dataset are calculated as the probability of values in the permutation distribution being higher than the cluster level statistics.

Once significant clusters are located with this multidimensional test, we perform pairwise comparisons of the relevant conditions to determine if the distances in state space between conditions in these time intervals are larger than those expected by chance. At each time point, we calculate the mean squared Euclidean distance between the trajectories of levels 1 and 2 of factor m, 1ΔttQ¯1(t)Q¯2(t)2. To generate the null hypothesis distribution for comparison, we calculate this same statistic based on the permutation procedure explained above. If factor m has only two levels, the projections of permutated data on principal components are the same as those generated above for the cluster mass calculation. If factor m has more than two levels, we generate new permutation projections on the principal components by permuting trials as above but only for the levels being compared. The mean squared Euclidean distance is compared to this distribution to calculate p-values. The mean squared Euclidean distance is subtracted from each value in this distribution to generate a new distribution, which is then used to calculate the mean and 95% confidence interval of the mean squared Euclidean distance above chance-level distance (Figure 6J–K). A sufficient number of permutations was used for these tests to allow for a minimum Holm-Bonferroni corrected p-value of 0.001 (1000 iterations for a single pairwise comparison and 3000 iterations for three comparisons).

Finally, we used this same framework to analyze firing rate differences between trial types within individual units. Here, we focus on the firing rate on each individual trial instead of the PSTH across all trials for the unit. Each trial contains only a single condition (e.g., Planned Turn). The single-trial firing rate at time t on trial i of condition m is ri,m(t) , resulting in the following sums of squares:

(8) SStot(t)=m=1Mi=1Nm(ri,m(t)r¯(t))2SSM(t)=m=1Nm(r¯m(t)r¯(t))2SSerror=SStotSSM

where Nm is the number of trials on condition m. The F-statisic is

(9) FM1,NM(t)=SSM(t)M1SSerror(t)NM

where N is the total number of trials. The 90th percentile of the F-distribution is used to select clusters whose cluster-level statistic is the sum of SSerror as in the repeated measures case. These cluster statistics are compared to those generated using a null distribution created by randomly partitioning the trials among the conditions keeping the trial count in each condition constant. This comparison results in a p-value for each cluster.

We used this method to quantify the number of units that reflect the firing rate differences identified on the population level. If a time interval contained a significant difference between a pair of conditions across the population of units, we determined if this time interval contained a significant difference between that pair of trial types within each unit. A significant difference was considered to be a cluster with p<0.1. Under this criterion, we would expect 10% of units to contain significant differences by chance alone; this is the null hypothesis. We calculated a p-value for this null hypothesis as the probability of obtaining at least as many significantly different units under the binomial distribution with success rate of 10% (this is the binomial test).

Mixed effect models

Request a detailed protocol

We tested if the firing rate of the recorded units changes as a function of the movement speed and acceleration. During trial , the hand position was used to update the cursor position (xi(t),yi(t)) on the screen. The cursor position was filtered using a Gaussian kernel (10 ms standard deviation) and used to calculate instantaneous speed and acceleration at time t:

(10) (vx,i(t),vy,i(t))=(xi(t)xi(tdt)dt,yi(t)yi(tdt)dt)(ax,i(t),ay,i(t))=(vx,i(t)vx,i(tdt)dt,vy,i(t)vy,i(tdt)dt)

where dt is the time interval between consecutive measurements. We then calculated the magnitude of the instantaneous speed and acceleration vectors on the screen:

(11) vi(t)=(vx,i(t),vy,i(t))ai(t)=(ax,i(t),ay,i(t))

where . is the Euclidean norm. To combine data across sessions, we normalized vi(t) and ai(t) by z-scoring across all trials in each session. Figure 4 shows averages of normalized speed and acceleration for each trial type, m, aligned to different task events.

To evaluate the relationship between the firing rate of each unit and movement kinematics, we used a lagged regression analysis:

(12) r¯m(t)=β0+β1v¯m(t+Δt)+β2a¯m(t+Δt)

where r¯m(t) is the mean firing rate of the unit at time t on trial type m, v¯m and a¯m are the mean movement speed and acceleration, respectively, and Δt is the time lag. The β parameters were fit as separate fixed effects for each unit. We systematically varied Δt, searching for the model with the highest R2 as the best fitting model. Across units and across Δt that were modeled, the firing rate time points that were fitted were kept constant. The maximum positive and negative Δt determined which data points were fitted and which were held out for model testing. Consider an r¯m(t) on the range [tstart,tend] and modeled on the range [Δtmin,Δtmax]. There is a subset of points in [tstart,tend] , such that there is data for v¯m(t+Δt) and a¯m(t+Δt) (i.e. t+Δt is also on the interval tstart,tend). Therefore, r¯m(t) was only fit on the range [tstartΔtmin,tendΔtmax] . The optimal value of Δt determined which time points were on this range and were, therefore, used to fit the model. The time points that were not fit but were on the range [tstartΔt,tendΔt] (note this range contains the best fit ) were used as a testing set for the models. Equation 12 was chosen as a good descriptive model of firing rates across all units after extensive search of possible linear models of speed and acceleration. Confidence intervals for R2 were calculated using bootstrapping (1,000 iterations).

Overall, we compared 14 alternative models that contained a combination of speed, acceleration, and their interaction (Figure 4—figure supplements 34). In some models, the x- and y-components of speed and acceleration in the direction of target or their absolute values were included instead of the total magnitude of speed and accleration. Our critical results were replicated in the alternative models that had the top out-of-sample fits to the data.

Ideal observer analysis

Request a detailed protocol

To determine whether population firing rate differences were sufficient for decoding of behavior, we used a cross-validated ideal observer analysis to determine if the firing rate changes across the population of units were sufficient to correctly classify the trial type. Our procedure is as follows. On each iteration, we randomly held out one trial of the same type from each unit (test set). Next, we calculated the mean response of each unit on the rest of the trials (training set) marginalized by the trial type. We then calculated the mean Euclidean distance of the firing rate patterns of the held out trials across the population from each of the patterns of marginalized mean responses for the training trials. The held out trials were classified according to the smallest Euclidean distance. We performed 10,000 iterations, randomly dividing the trials for each unit into training and test sets. The decoding accuracy was the fraction of iterations correctly classified.

LFP and joint LFP-spiking analysis

Request a detailed protocol

For LFP analyses, the raw 44 kHz recorded voltage signals from each electrode were downsampled to 2 kHz. We then locally detrended and removed line noise using the Chronux toolbox (Mitra and Bokil, 2008), augmented with manual artifact removal where needed. The LFP from 2.5 s before to 2.5 s after an artifact was censored from further analysis. We then calculated the continuous wavelet transform using a Morse wavelet with symmetry parameter of 3 and time-bandwidth product of 20. We aligned the wavelet-transformed data across trials for each session, averaged across all trials of the same type (in the complex plane), and calculated the squared magnitude to determine the evoked power. To calculate the induced power, we computed the squared magnitude of the aligned wavelet-transformed data, averaged across trials, and subtracted the evoked power.

Evoked and induced power calculations are sensitive to varying trial counts. In our study, each condition had 4–47 trials depending on the session. Furthermore, because of variability in the timing of events in our task, the number of trials available for analysis varied by time. To control for potential noise induced by this variability, we used the time period between 500 ms to 250 ms before the fixation onset to calculate a baseline for the average evoked and induced power for each trial. We then used a permutation procedure where we sampled without replacement from trials in each session, without regard to trial type, in order to estimate the average powers for a variety of trial counts. Repeating this process one hundred times generated a baseline distribution of evoked and induced power for any trial count at each time point. These baselines were separately calculated for each recording session. We normalized our calculated evoked and induced power values at each point in time by dividing by the mean of the appropriate baseline distribution for the trial count of that particular trial type at that particular point in time in each session.

We also calculated evoked and induced power of the beta frequency range. To do this we applied a notch filter between 13 and 30 Hz to the downsampled, detrended voltage data. We then used the Hilbert transform to calculate the analytic signal. Evoked and induced power were calculated from this signal and normalized by the baseline for the trial counts in the session.

The coupling of spiking activity to the LFP was quantified using the pairwise-phase consistency metric (PPC). We use the P^2 metric described by Vinck et al., 2012 which is not vulnerable to biases caused by varying firing rates or dependencies between the firing rate and the distribution of spike phases. Briefly, we used the continuous wavelet transfrom to calculate the phase of each spike at every frequency of the LFP, termed the spike-triggered phase. The spikes were then aligned to task events and segmented into bins. For each time point, PPC was calculated using spikes that preceded that point in time by no more than 10 cycles; bin size, therefore, varied with different frequencies. We capped the bin size at 1.5 s, and spaced time points by 50 ms. We randomly selected a pair of trials of the same type (e.g. two Reach trials). For each pair of spikes in the same bin across the two trials, we calculated the dot product of spike-triggered phases and averaged the results across the two trials. These averaged dot products were averaged over all pairs of trials of the same type to calculate PPC for a particular bin. This was repeated for all bins allowing us to calculate a spectrogram for PPC. The PPC for the entire beta range was calculated in the same way using the the spike-triggered phase from the analytic signal of the beta-filtered LFP instead of the continous wavelet transform.

We examined the data for differences between spectrograms using the cluster mass test in the same fashion as explained above for firing rates. No significant clusters were located. As explained in the section on Statistical Analysis, we avoided analyzing arbitrarily selected intervals to control for Type I error.

Because the shortest of our event-aligned analysis intervals was approximately 500 ms, we limited our frequency-based analysis to above 4 Hz.

Movement trajectory analysis

Request a detailed protocol

To confirm that subjects performed the task as instructed, we generated average movement trajectories across trial types. Since movements were self-paced and could take different times, averaging the hand position across trials at a particular time would not accurately reflect the variability of trajectories across space. To calculate the mean and standard deviation of movement trajectories in a particular trial type, we devised the following analysis using pairs of trajectories. The mean of an arbitrary vector {xk} is:

μ=1Nk=1Nxk

which can be rewritten as the mean of the mean of independently-sampled pairs of data points, i and j:

(13) μ=12Nk=1Nxk+12Nk=1Nxk=12N(i=1Nxi+j=1Nxj)=1N2i=1Nj=1Nxi+xj2

To calculate the mean of an arbitrary number of trajectories, we used dynamic time warping (Sakoe and Chiba, 1978) to align each pair of trajectories, and took the mean of the aligned positions. We then selected 200 data points from this mean over equally spaced time intervals. This is the pairwise mean trajectory. Finally, we calculated the mean across pairwise mean trajectories to calculate the mean trajectory.

The variance can be calculated from pairs of trajectories in a similar way:

(14) σ2=1Nk=1N(xk1Nl=1Nxl)2=1Nk=1N(1Nl=1Nxk1Nl=1Nxl)2=1N3k=1N(l=1Nxkxl)2

The method above was used only to calculate the mean and standard deviation of movement trajectories shown in Figure 1E.

To classify single-trial movement trajectories, we performed hierarchical clustering of the trajectories using the discrete Frechet distance (Eiter et al., 1994) as a dissimilarity metric between pairs of trajectories. We then used Ward’s method (Ward, 1963) to generate the cluster tree.

Data availability

All raw and processed data has been uploaded to Dryad. DOI:https://doi.org/10.5061/dryad.2jm63xsq2.

The following data sets were generated
    1. London D
    2. Fazl A
    3. Katlowitz K
    4. Soula M
    5. Pourfar M
    6. Mogilner A
    7. Kiani R
    (2021) Dryad Digital Repository
    Distinct population code for movement kinematics and changes of ongoing movements in human subthalamic nucleus.
    https://doi.org/10.5061/dryad.2jm63xsq2

References

    1. Anderson MJ
    (2001) Permutation tests for univariate or multivariate analysis of variance and regression
    Canadian Journal of Fisheries and Aquatic Sciences. Journal Canadien Des Sciences Halieutiques et Aquatiques 58:626–639.
    https://doi.org/10.1139/f01-004
  1. Software
    1. Eiter T
    2. Eiter T
    3. Mannila H
    (1994)
    Computing discrete frÉchet distance
    Computing Discrete FrÉchet Distance.

Decision letter

  1. Nicole C Swann
    Reviewing Editor; University of Oregon, United States
  2. Richard B Ivry
    Senior Editor; University of California, Berkeley, United States
  3. Andrew Sharott
    Reviewer; University of Oxford, United Kingdom
  4. Harrison Walker
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This manuscript reports on a unique data set of human subthalamic nucleus electrophysiological recordings during a behavioral task. The results (which include single unit, multi-unit, and local field potential activity) point to the complex nature of the neural code in this structure and challenge some existing models related to the role of the subthalamic nucleus.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Distinct population code for movement kinematics and changes of ongoing movements in human subthalamic nucleus" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Andrew Sharott (Reviewer #1); Harrison Walker (Reviewer #3).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

While the reviewers were overall positive about many aspects of the manuscript, there were a number of points they brought up. Upon discussion, we felt that the changes required were substantial enough that we were not confident they could be fully addressed in a revision. In particular, the reviewers noted missing information about the type of patient(s) studied and noted that this could have a substantial impact on the interpretations of the results. Additionally, we felt that important information about number of units contributed by each patient (and how those may have influenced different STN sub-populations was not clear). Finally, the reviewers felt additional statistical analyses and justification were necessary. While the reviewers agreed that it might be possible for all these points to be addressed in a revision, they also noted it may substantially change the interpretation (depending on the updated results). The full reviewers are included below.

Reviewer #1:

1. There are important details missing regarding the data and data acquisition. What type of patients were these recordings made in? (a) Are they from different patient populations? This is an important consideration for several aspects of the study including the clustering of neurons into different response types and the analysis of β oscillations. (b) What is the baseline firing rates of the STN neurons and how many were considered single units? The homogeneity of STN unit spiking means that single units should be in consistently firing at a relatively high rate (10-35 spikes/s). Defining the entrance, the STN alone is not sufficient to establish that the units in the STN (particularly at the ventral border). Further evidence that the units were in the STN and were single units is important to substantiate the authors claims as to the relevance of their results to STN function.

2. On a related issue, do the authors have information as to the spatial parameters of the STN recordings? If so, do these parameters correlate with the function clustering of the response types?

3. The authors use custom use a range of custom designed data analyses to tease out STN neurons with tendencies to respond to differently to parts of the task in the different conditions. In some cases, this is a strength of the study and many of the analyses elegantly pull out features of the data that may not be immediately clear (e.g Figures3 and 4). However, in it's current form it is difficult to establish the effect sizes for other analyses – particularly in the case of the "Turn Units." While the authors go to impressive lengths to show that the differences at the Turn are different, from Figure 3C they seem to be very small and a fraction of those seen in response to kinematic parameters. Currently, many of the novel claims are dependent on this difference (Abstract: "We show that STN activity decreases during action switches, contrary to prevalent theories".). Can the authors provide evidence of the effect size of the responses to task switching as well as being able to say they are significantly different? How many individual units are able to distinguish between specific aspects of the task (such as the example in Figure 2C).

4. Given the size of the sample – the authors should provide some evidence that their main findings are scene across the majority of patients. For example, how may units in each cluster are derived from different patients?

5. The analysis of β oscillations/synchrony is potentially very interesting, but I cannot interpret it's relevance without knowing the disease-status of the patients.

Reviewer #2:

The authors report on subthalamic nucleus neural population dynamics and its role in movement planning and execution, especially as it relates to changes in the motor plans (they compared trials of straight reach, planned turn or unexpected turn while recording 39 single units from 8 patients undergoing DBS surgery). They used principal component analysis of firing rates and unsupervised clustering to identify two types of cells – movement and turning units, each contributing specific information regarding movement dynamics encoded in the STN (the third identified cluster was too small to analyze). The manuscript is generally well written and the complex analysis is logically presented although sometimes conclusions are difficult to appreciate from the figures. This type of analysis is increasingly necessary to try and understand the complex dynamics in the nervous system, but the challenge will be to relate to behavior in a meaningful way. The authors should relate their findings to (presumably) parkinsonian pathology especially when discussing synchronization of units with β band oscillations which is the hallmark of the disease.

– The analysis is presented only for move and turn onsets. Was there a difference in firing rates at Go and Impromptu Turn cues? Also, when was the impromptu turn cue given and did this differ by movement speed (within and across subjects)?

– "This result was not due to lack of statistical power of the cluster mass test, which showed great sensitivity to systematic trends in our datasets." – It is difficult to assess such statements which occur several times. Could you explain what trends?

– "Neuronal spiking in STN is typically phase-locked to β oscillations (Kühn et al., 2005; Weinberger et al., 2006)" – is this because of parkinsonian state? The authors need to acknowledge that subjects had Parkinson's disease (presumably, this is not stated anywhere) so findings may be influenced by pathological parkinsonian state, especially in regards to single cell firing being phase locked to β oscillations which is known to be exaggerated in PD.

– "Interestingly, these sub-populations are partially spatially segregated within STN, with Movement Units located more laterally…" – provide a depiction of unit locations if commenting on this in the Discussion.

– "Further, stimulation of glutamatergic cortical inputs to STN (i.e., the hyperdirect pathway) worsens bradykinesia (Gradinaru et al., 2009)" – that paper showed quite the opposite – stimulation of hyperdirect pathway improved bradykinesia. Please address in your conclusions.

– "To compare spiking activity around the turn on the Impromptu Turn trials with spiking activity on Reach and Planned Turn trials, we sampled randomly from the distribution of turn cue onsets on Impromptu Turn Trials to define corresponding analysis windows on Reach and Planned Turn trials."

– I understand why you had to do this for Reach trials, but why did you not use the actual turn on Planned turn trials as the analysis window?

– Methods: "Subjects were patients undergoing implantation of deep brain stimulation…" Please specific that these were presumably patients with Parkinson's disease, and were they off dopaminergic medications?

– Authors report "13.5{plus minus}3.3 trials per trial type per subject" – How may trials of each type were done for each unit recorded? How many units were recorded per patient?

Reviewer #3:

Conceptually, even if 'impromptu' units are an actual subpopulation of neurons, a decrease in STN activity just before the execution of the unexpected turn seems consistent, not inconsistent, with previous hypotheses on STN as a "brake" or "global stop" nucleus. I think they could clarify their thoughts on this more in their writing.

Would like to see more specific examples of sorted units and examples that clearly demonstrate the phenomena they are claiming.

All group level plots and panels should include the number of unique participants, the number of individual units, and the number of trials. For example, in Figure 2e PC1 appears to just be the PTH from the unit in Figure 2C. This unit appears to overwhelmingly dominate the population variance. Was PCA performed across all subjects/trials/units?

Given 39 units / 9 STNs = 4 units per STN. How many subjects had multiple recording sites, what was the average number of units per recording site?

Plots should generally be aligned to impromptu turn cue onset. Without this time alignment, it’s unclear whether the observed changes are related to the cue, the prior movement, or the subsequent movement. Even with this alignment there are some questions about prior events spilling into the 'impromptu' turn condition, depending on the latency of the second command and the nature of the kinematic response of the units in a given experiment.

Figure 2B/D: How do these trajectories differ? Following feedback, the traces look identical. Were stats performed across all principal components?

Figure 3C: Turn versus impromptu turn changes are < 0.2 Z-scores – these are not convincing changes.

Use of terms ipsilateral and contralateral is confusing. I think they mean medial/lateral or flexor/extensor as these tasks appear to be occurring in the same hand.

Figure 7E-F – description says that turn trials show a decreased spike synchronization, but this is because β disappears. Also it appears that this happens regardless of trial type (not just during turns). They might consider also including high frequency LFP activity for comparison with β.

The authors should explicitly state what disease the diagnoses of the patients (assuming Parkinson's disease in all, since these are STN implants).

Supplemental Video with 3D PCA visualization, why are the plotted dimensions PC1, PC3, and PC4? Where is PC2?

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

Thank you for resubmitting your work entitled "Distinct population code for movement kinematics and changes of ongoing movements in human subthalamic nucleus" for further consideration by eLife. Your revised article has been evaluated by Richard Ivry (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:

Primarily we had some concerns about the interpretation of the model especially with response to the "Turn units" and what they may actually be encoding. We have further elaborated on this concern, as well as highlighting a few additional concerns below.

Essential revisions:

1) Please comment more on what the 'turning' units are actually encoding, and if they could actually be encoding the acceleration and velocity in the vertical axis? Further, do the results still support the statement that 'one sub-population encoding movement kinematics and direction and another encoding unexpected action switches'. Specifically, Figure 4 indicates that both 'Movement Unit' and 'Turn Unit' firing rate can reflect kinematic parameters with different time lags. The authors argued that (Line 276-278): '… the strength of this relationship (for Turn Units) was considerably weaker than for Movement Units (mean R2 = 0.45 compared to 0.61 for Movement Units). Furthermore, movement kinematics did not accurately predict the PSTHs of these units (Figure 4 —figure supplement 2B, D, and F).' However, it is quite subjective to say the R2 values are 'considerably weaker'. Comparing the plots of A, C and E against B, D, and F in Figure 4 —figure supplement 2, can the author's confidently say that movement kinematics can predict Movement Units (shown in plots of A, C and E), but cannot predict Turn Units (shown in plots of B, D and F)? Furthermore, even if the model fitting is a bit worse for Turn Units compare to Movement Units, did the authors used the same model that fits best for movement units? Or to be more fair, did the authors do model comparison for best fitting models for the two types of Units separately?

2) The authors mentioned that 'Turn Units decreased their firing rates before turn onset on Impromptu Turn trials but increased them before turn onset on Planned Turn trials.' (Line 289-291). However, we don't see this statement being supported by data/figure. Can the authors provide evidence/data for this? Or did we miss anything?

3) Reviewer 3 from the previous round of review raised that 'Conceptually, even if 'impromptu' units are an actual subpopulation of neurons, a decrease in STN

activity just before the execution of the unexpected turn seems consistent, not inconsistent, with the previous hypotheses on STN as a "brake" or "global stop" nucleus. We think the authors should clarify their thoughts on this more in their writing.' We feel this was not sufficiently addressed.

4) For the statistical analysis Equation. 1 – Equation. 9, it is not clear whether the authors have used self-programmed scripts for these tests or existing toolbox/functions in Matlab to do this. As far as I understand, Equation. 1 – Equation. 9 are similar to what is implemented in the linear mixed effects models available in Matlab (https://www.mathworks.com/help/stats/linear-mixed-effects-models). This can be applied to each time point, and then using cluster mass test to identify the time window with significant effect. Or alternatively, can the authors specify what is difference between the method used in this study and the linear mixed effects model? The authors mentioned that 'For testing factor ;, the levels of ; are permuted within each unit and without permutation of the other factor. The remainder of the calculation continues as for the one factor case, resulting in separate p-values for each main effect, as in a multi-way ANOVA.' Is it equivalent to repeating the multilevel modeling permuting the levels of m? If the method is based on the framework of multi-level modeling, it might not be necessary to present all the Equations.

5) The accurate rate of different conditions was around 0.7. Did the authors only focus on the 'accurate trials' which satisfy the 'accuracy criteria' or considered all trials? Is there any difference between the neural activities for 'accurate trials' vs 'inaccurate trials'?

7) Can the authors compare the baseline spiking rate of the 'movement related' and 'turning related' units identified by clustering analysis applied on the PSTHs.

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

Author response

[Editors’ note: The authors appealed the original decision. What follows is the authors’ response to the first round of review.]

Reviewer #1:

1. There are important details missing regarding the data and data acquisition. What type of patients were these recordings made in? (a) Are they from different patient populations? This is an important consideration for several aspects of the study including the clustering of neurons into different response types and the analysis of β oscillations. (b) What is the baseline firing rates of the STN neurons and how many were considered single units? The homogeneity of STN unit spiking means that single units should be in consistently firing at a relatively high rate (10-35 spikes/s). Defining the entrance, the STN alone is not sufficient to establish that the units in the STN (particularly at the ventral border). Further evidence that the units were in the STN and were single units is important to substantiate the authors claims as to the relevance of their results to STN function.

Many thanks for your thoughtful comments, which helped us close the gaps in the paper and improve it in the process. All recordings are from a homogeneous population – patients with medically refractory Parkinson’s disease undergoing surgical implantation of DBS electrodes in the STN. They did not take their anti-Parkinsonian medications on the day of surgery. We regret omitting the description of the disease status. The disease status has now been explicitly noted in the abstract, introduction, results, methods, and discussion. Patient age, pre-operative UPDRS motor score, side of recording, and number of single and multi-units recorded per patient are listed in a new table (Supplementary File 1) which we reference in the results. Our conclusions are not changed if any of subjects 1, 2, or 7, who contributed relatively more units than others, are removed from the data (page 43).

The characteristics of all units are now listed in a new table (Supplementary File 2). We were conservative in classifying single-units and chose only units with clearest isolation from background and consistency of action potential waveforms (8/39). The single units had a baseline median firing rate of 32.9 Hz (IQR: 18.3-53.9), consistent with single units identified in past studies (Sharott et al., 2014, Weinberger et al., 2006). Furthermore, the median firing rate of multi-units was 40.4 Hz (IQR: 20.7-82.6), consistent with the results of Weinberger et al., suggesting that most multi-units are largely contributed to by a single neuron. The baseline firing rate was calculated from spikes in the intervals between trials. We discuss this on page 20.

We analyzed both single- and multi-units as is commonly done (Sharott et al., 2014, Zavala et al., 2017, Lipski et al., 2018), and our analyses do not depend on the units being well-isolated. We have largely consistent tuning across single and multi-units, and because all analyses were performed on Z-scored firing rates, they did not depend on differences in the range of firing rates. Critically, single units individually display characteristics seen in the population of all units. We now show in the revised Figure 3 which neural responses are from single units (white asterisks in Figure 3A). Units in this figure are sorted by proximity in the hierarchical clustering algorithm. If our results do not depend on whether a unit is well-isolated we would expect single units to not be clustered together by the algorithm, as is the case in the figure (discussed on page 6).

Microelectrodes were targeted to the ventral border of the STN based on structural MRI and the AC-PC coordinates (11 mm lateral to midline, 5 mm inferior to the AC-PC plane, and 4 mm posterior to the AC-PC midpoint). However, the recorded units spanned the space between the dorsal and ventral borders of STN. Entrance to the STN was noted based on two criteria: (1) real-time alignment of a CT scan performed at the time of surgery with a high-resolution pre-surgical structural MRI (for a detailed explanation, please see our response to the next comment), and (2) an increase in background spiking and identification of rhythmically firing cells at high firing rates. A decrease in background spiking marked the ventral boundary of STN. In many of our surgeries deeper cells along this trajectory fired rapidly with a relatively constant firing rate, characteristic of SNr cells. This procedure has been previously described in detail by our group (Sterio et al., 2002). We have added these points to the Methods (page 19).

2. On a related issue, do the authors have information as to the spatial parameters of the STN recordings? If so, do these parameters correlate with the function clustering of the response types?

We used the pre-operative high resolution structural MRI with gadolinium, pre-operative high resolution CT scan with the Leksell frame in place, planned trajectory, and intra-operative or post-operative CT scan with DBS lead in place to localize our recordings. Our procedure is as follows. The pre-operatively planned target in AC-PC coordinates was identified on MRI. This MRI was merged with the pre-operative frame-localizing CT scan. Using the frame localization, the target was converted into frame coordinates. The pre-operatively planned trajectory was identified using this target and the ring and arc angle used. Based on the alignment of the post-operative CT scan with the DBS lead with the intra-operative CT scan, the trajectory was confirmed and corrected whenever necessary (corrections were always <1mm). The distance of the microelectrode from the ventral STN target was measured using our microdrive. The location of each recording was identified as the point along the trajectory with the measured distance from target. In the case of multiple simultaneously used microelectrodes, the additional trajectories were offset by the distance of electrodes (2 mm) in the appropriate direction. These calculations were performed using the Brainlab Elements stereotaxy software. This description has been added to the Methods (page 19).

We have added Figure 2 —figure supplement 2 to show the locations of all recorded units in AC-PC coordinates (A) and in coordinates normalized to the electrophysiological entry point of STN (B, we followed the method of Sharott et al., 2014). There was a nonsignificant trend toward movement units being located more dorsolaterally in raw AC-PC coordinates (A) but this was not evident when examining the normalized coordinates. We state on page 6 that there was no clear association with unit classification and recording location.

3. The authors use custom use a range of custom designed data analyses to tease out STN neurons with tendencies to respond to differently to parts of the task in the different conditions. In some cases, this is a strength of the study and many of the analyses elegantly pull out features of the data that may not be immediately clear (e.g Figures3 and 4). However, in it's current form it is difficult to establish the effect sizes for other analyses – particularly in the case of the "Turn Units." While the authors go to impressive lengths to show that the differences at the Turn are different, from Figure 3C they seem to be very small and a fraction of those seen in response to kinematic parameters. Currently, many of the novel claims are dependent on this difference (Abstract: "We show that STN activity decreases during action switches, contrary to prevalent theories".). Can the authors provide evidence of the effect size of the responses to task switching as well as being able to say they are significantly different? How many individual units are able to distinguish between specific aspects of the task (such as the example in Figure 2C).

4. Given the size of the sample – the authors should provide some evidence that their main findings are scene across the majority of patients. For example, how may units in each cluster are derived from different patients?

While we focused our manuscript on population-level results, these results are indeed evident when examining individual units. Author response table 1 shows the individual unit responses for each population level result. For each unit, we used the cluster mass test to determine if any significant differences were present across the time interval over which population-level changes are seen. For example, Turn Units have decreased firing on Impromptu Turn trials compared to Planned Turn trials between 240 ms before and 80 ms after the turn. The firing rate across this interval contained significant differences between Planned and Impromptu Turn trials on 33% of Turn Units (6/18), and the rest of the units showed a trend in the same direction although they could not reach significance due to limited data for each unit. Our threshold for significance at the individual unit level was p=0.1, as is commonly done (e.g., Mallet et al., 2016). Under this criterion, we would expect 10% of units to show significance under the null hypothesis of no difference between trial types. The effect was seen on 33% of neurons allowing us to confidently reject this null hypothesis (p=0.0012, binomial test). The population level results for firing on medial and lateral trials (we have changed our terminology from ipsilateral and contralateral at the suggestion of Reviewer 3) on Movement Units similarly holds at the individual unit level (and has quite a large effect size, statistically significant in 67% of units). These results are now included in the revised manuscript (pages 7, 9, 11).

Author response table 1
Population Level ResultFraction of Individual Unitsp-value
Decreased firing rate of Movement Units on Planned Turn trials compared to Impromptu Turn trials, 140 ms after movement onset to 130 ms before turn onset (Figures 2B, 3B)0.38 (6/16)5.0 × 10-4(binomial test)
Decreased firing rate of Movement Units on Planned Turn trials compared to Reach trials 140 ms after movement onset to 130 ms before turn onset (Figure 2B, 3B)0.13(2/16)0.21(binomial test)
Decreased firing rate of Turn Units on Impromptu Turn trials compared to Planned Turn trials 240 ms before to 80 ms after turn onset (Figure 2D, 3C)0.33(6/18)0.0012(binomial test)
Decreased firing rate of Movement Units on lateral trials compared to medial trials 10 ms before to 270 ms after movement onset (Figure 5)0.67(8/12)1.7 × 10-7(binomial test)

While we acknowledge that our sample size of patients was small, the results were consistent across them. The table of unit properties (Supplementary File 2) shows the number of units contributed by each patient and the classification of the units. Movement units were seen in 75% of patients (6/8), and turn units were seen in 75% (6/8), as well. All patients had recordings of Movement and/or Turn Units (i.e., no patients contained zero Movement Units and Turn Units). This is reported on pages 4-5 of the revised manuscript.

5. The analysis of β oscillations/synchrony is potentially very interesting, but I cannot interpret it's relevance without knowing the disease-status of the patients.

All patients had Parkinson’s disease. Their age and pre-operative UPDRS scores are now reported in Supplementary File 1 of the revised manuscript. We now discuss these results in the context of Parkinson’s disease (page 16) , suggesting that our conclusions are more broadly applicable given concordant results between studies of STN in humans with Parkinson’s disease or OCD, and non-pathologic recordings in non-human primates and rodents.

Reviewer #2:The authors report on subthalamic nucleus neural population dynamics and its role in movement planning and execution, especially as it relates to changes in the motor plans (they compared trials of straight reach, planned turn or unexpected turn while recording 39 single units from 8 patients undergoing DBS surgery). They used principal component analysis of firing rates and unsupervised clustering to identify two types of cells – movement and turning units, each contributing specific information regarding movement dynamics encoded in the STN (the third identified cluster was too small to analyze). The manuscript is generally well written and the complex analysis is logically presented although sometimes conclusions are difficult to appreciate from the figures. This type of analysis is increasingly necessary to try and understand the complex dynamics in the nervous system, but the challenge will be to relate to behavior in a meaningful way. The authors should relate their findings to (presumably) parkinsonian pathology especially when discussing synchronization of units with β band oscillations which is the hallmark of the disease.

We very much appreciate the reviewer’s comments and have used them as an opportunity to improve the paper, as we explain below.

– The analysis is presented only for move and turn onsets. Was there a difference in firing rates at Go and Impromptu Turn cues? Also, when was the impromptu turn cue given and did this differ by movement speed (within and across subjects)?

Figure 3 —figure supplement 1 shows average neural activity around the go and turn cues (A and B) compared to the movement and turn onset (C and D) for both unit types. There was no significant difference between trial types around the go cue. However, a clear increase in firing rate is visible beginning about 200 ms before movement onset and about 300 ms after the go cue, suggesting that this response is likely related to movement or movement planning, not the processing of the go cue itself.

The turn cue was presented after the movement onset and before the hand reached the target position, prompting subjects to alter their movement trajectory mid-movement. Turn Units responded similarly in different trial types around the time when the turn cue appeared (B). Movement Units had significantly smaller firing rates on Planned Turn trials than Reach trials between 50 ms before and 170 ms after the expected time of the turn cue (A) because subjects had already begun their turn on Planned Turn trials and movement units track hand kinematics (note that subjects observed the actual turn cue only on the Impromptu Turn trials, but we can compare firing rates around that time across trials).

Because our goal in this study was to determine neural responses in STN to planned and unplanned movements, our key comparisons are with respect to turn onset not the turn cue. Past human studies in STN frequently used go/no-go or stop-signal tasks with button presses or joystick movements as the operant response. Reporting cue-aligned responses is reasonable in those cases. However, our task has the advantage of recording a full, extended movement, which could also be multi-segmented (planned and impromptu turn trials). We can thus parse when exactly movement onset and turn onset occur. Indeed, this is where we also found significant modulation of firing rates.

The turn cue on Impromptu Turn trials was presented when subjects were within a specific rectangular window on their path to the reach target. The center of this window was randomly chosen from a uniform distribution between 35% and 65% of the distance between the fixation and target points. Since the time the turn cue was presented was determined by the subject’s position, faster movements were more likely to result in a sooner presentation of the turn cue. However, as explained in the paper, Turn Units were minimally sensitive to changes in movement in kinematics, including speed and acceleration. The turn cue was presented 206 ± 71 ms after movement onset and 504 ± 154 ms after the go cue (mean ± SD across sessions).

To clarify these points, we have added the figure above to the manuscript (Figure 3 —figure supplement 1) and explain the findings and the conceptual advance on pages 7 and 11.

– "This result was not due to lack of statistical power of the cluster mass test, which showed great sensitivity to systematic trends in our datasets." – It is difficult to assess such statements which occur several times. Could you explain what trends?

We apologize for being vague and have now clarified our approach and results in the revised manuscript. A decoder trained to classify population neural responses performs with ∼75% accuracy (cross-validated) for differences of 0.2 z-scores and performs above chance for differences as low as 0.05 z-scores. These analyses give us confidence that absence of significance is not simply due to low power. We have now included the accuracy with which a decoder classifies population neural responses for individual results in the manuscript.

– "Neuronal spiking in STN is typically phase-locked to β oscillations (Kühn et al., 2005; Weinberger et al., 2006)" – is this because of parkinsonian state? The authors need to acknowledge that subjects had Parkinson's disease (presumably, this is not stated anywhere) so findings may be influenced by pathological parkinsonian state, especially in regards to single cell firing being phase locked to β oscillations which is known to be exaggerated in PD.

Thanks for bringing up this important point. We regret the omission of the fact that all of our patients had medically refractory Parkinson’s disease and were off medication on the day of surgery. While our study took place in these patients that are known to have underlying hypersynchronous β activity (Brittain et al., 2014), β oscillatory activity in STN is not a unique phenomenon in Parkinson’s disease; β activity in the STN has been found in patients with OCD (Rappel et al., 2018). Furthermore, β band activity increases during stopping in both Parkinson’s (Ray et al., 2012) and OCD (Bastin et al., 2014), suggesting that the functional role of β oscillations in action inhibition is not unique to Parkinson’s disease. It is, therefore, notable that we did not find an increase in β induced power or spike-field locking prior to subjects interrupting their planned behaviors. However, our results may still be influenced by a pathological Parkinsonsian state consisting of elevated β activity and generalizing them would require recordings in non-Parkinsonian patients or in animal models without underlying disease. We discuss these points on page 16 of the discussion.

– "Interestingly, these sub-populations are partially spatially segregated within STN, with Movement Units located more laterally…" – provide a depiction of unit locations if commenting on this in the Discussion.

We now provide our methodology for locating recorded units on page 19 of the revision, and we have included Figure 2 —figure supplement 2, which depicts unit locations. There was a nonsignificant trend toward movement units being located more dorsolaterally in raw AC-PC coordinates but this was not evident when examining the normalized coordinates (normalization method was adopted from Sharott 2014). While there may be some segregation of unit types, this is not a strong trend and we have removed this claim from our results and discussion.

– "Further, stimulation of glutamatergic cortical inputs to STN (i.e., the hyperdirect pathway) worsens bradykinesia (Gradinaru et al., 2009)" – that paper showed quite the opposite – stimulation of hyperdirect pathway improved bradykinesia. Please address in your conclusions.

Gradinaru et al.’s findings regarding the effects on bradykinesia caused by stimulation of the hyperdirect pathway were complex. Gradinaru et al. explain that stimulation of the hyperdirect pathway from motor cortex to STN improved bradykinesia, but stimulation of glutamatergic inputs in STN resulted in differing effects depending on the frequency of stimulation. We have now corrected our statement in the revised manuscript (page 18)

– "To compare spiking activity around the turn on the Impromptu Turn trials with spiking activity on Reach and Planned Turn trials, we sampled randomly from the distribution of turn cue onsets on Impromptu Turn Trials to define corresponding analysis windows on Reach and Planned Turn trials." – I understand why you had to do this for Reach trials, but why did you not use the actual turn on Planned turn trials as the analysis window?

Thank you for the great question. Our wording here was not clear enough in distinguishing analyses aligned to the turn cue and to the turn movement itself. The turn cue is not present on Reach or Planned Turn trials so we sampled from the turn cues on Impromptu Turn trials to generate analysis windows for other trial types. For analyses of the turn movement itself, we used the distribution of the times of onset of turn movement in planned turn and impromptu turn trials to define temporally comparable analysis windows for Reach trials on which there was no turn. We have corrected the confusion caused by our wording in the revision (page 20-21).

– Methods: "Subjects were patients undergoing implantation of deep brain stimulation…" Please specific that these were presumably patients with Parkinson's disease, and were they off dopaminergic medications?

All recordings are from a homogeneous population – patients with medically refractory Parkinson’s disease undergoing surgical implantation of DBS electrodes in the STN. They did not take their anti-Parkinsonian medications on the day of surgery. The description of the disease status is now specified in the abstract, introduction, results, discussion, and methods. Patient age, pre-operative UPDRS motor score, side of recording, and number of single and multi-units recorded per patient are now listed in Supplementary File 1. We also confirm that our conclusions do not critically depend on subjects 1, 2, or 7, who contributed relatively more units than others (page 43).

– Authors report "13.5{plus minus}3.3 trials per trial type per subject" – How may trials of each type were done for each unit recorded? How many units were recorded per patient?

The number of units recorded per patient is shown in Supplementary File 1. The number of trials of each type per unit is shown in Supplementary File 2. Please also see our response to the first comment of Reviewer #1, which includes copies of the tables in Supplementary Files 1 and 2.

– Could not play the second movie

Reviewer #3:

Conceptually, even if 'impromptu' units are an actual subpopulation of neurons, a decrease in STN activity just before the execution of the unexpected turn seems consistent, not inconsistent, with previous hypotheses on STN as a "brake" or "global stop" nucleus. I think they could clarify their thoughts on this more in their writing.

Many thanks for bringing up this key point, which we now explain more extensively in the revised manuscript. Briefly, theories positing STN as mediating stop processes predict that stopping is accompanied by increased firing rates. Schmidt et al. 2013 is an example of this, where STN spiking immediately follows a stop cue. In our task, movement speed decreases before the turn on both Planned and Impromptu Turn trials, though more so for the latter (Figure 4 —figure supplement 1). If the Turn Units contribute to stopping, we would expect them to have increased firing rates on Impromptu Turn and Planned Turn trials compared to Reach trials and also higher firing rates on Impromptu Turn than Planned Turn trials. Our data indicate the opposite result; Turn Unit firing rates decrease before impromptu turns and increase before planned turns.

Would like to see more specific examples of sorted units and examples that clearly demonstrate the phenomena they are claiming.

Absolutely. In our revision, we now provide Figure 2 —figure supplement 1, which shows 3 example Turn Units and 3 example Movement Units. S and U indicate the subject and unit identity, and Nt indicates trial counts. Asterisks above the time axis indicate periods with statistically significant differences between Planned and Impromptu Turn trials. Although the dynamics of responses vary across units, our main findings are apparent in single units. The Movement Units show increased firing beginning around the movement onset and returning to baseline after the feedback. Meanwhile, the Turn Units show decreased firing in the period around the onset of turn movement on Impromptu Turn trials. The units in panels B and D reach statistical significance, and the unit in F shows a trend in the same direction as the majority of the population (please see our answer to comment 3 of reviewer 1 for a count of the individual units with significant statistics).

All group level plots and panels should include the number of unique participants, the number of individual units, and the number of trials. For example, in Figure 2e PC1 appears to just be the PTH from the unit in Figure 2C. This unit appears to overwhelmingly dominate the population variance. Was PCA performed across all subjects/trials/units?

We have now added the information to the figures. Author response table 2 summarizes the subject, unit, and trial counts for in each figure.

Author response table 2
FigureUnique subjectsUnitsTrials
18N/A913
2A,C (example unit)11107
2B,D (example unit)1145
2E-H8392550
3A8392550
3B6161117
3C6181177
4 A,C,E,G6161117
4 B,D,F,G6181177
5312970
63191501
7 A,C,E,G,I6161117
7 B,D,F,H,J6181177
2 Supplement 1 (example units)46402
2 Supplement 2631N/A
3 Supplement 1 A,C6161117
3 Supplement 1 B,D6181177
4 Supplement 18392550
4 Supplement 2 A,C,E,G6161117
4 Supplement 2 B,D,F,H6181177
4 Supplement 36161117
5 Supplement 1 A,C,E,G,I,K312970
5 Supplement 2 B,D,F,H,J,L38537
7 Supplement 18N/A1318

Regarding the specific example of Figure 2, Figure 2C shows an individual unit while 2E-H show the firing rates of all units projected on the principal components. The first principal component shown in Figure 2E was recapitulated in many individual units, including the example unit shown in 2C, but the PCs were not based on any individual unit or overwhelmingly shaped by a single unit. The example unit does not dominate the variance; it is representative of it. In fact, removing this unit for the calculation of PCs would make little difference in the outcome. PCA was done on PSTHs of correct trials on all units across all subjects. The PCA procedure is described on page 21 of the manuscript.

Given 39 units / 9 STNs = 4 units per STN. How many subjects had multiple recording sites, what was the average number of units per recording site?

We obtained spiking data from 20 recording sites (on average 1.95 units per site). 62.5% of subjects (5/8) had multiple recording sites. This has been added to the Methods (page 21).

Plots should generally be aligned to impromptu turn cue onset. Without this time alignment, it’s unclear whether the observed changes are related to the cue, the prior movement, or the subsequent movement. Even with this alignment there are some questions about prior events spilling into the 'impromptu' turn condition, depending on the latency of the second command and the nature of the kinematic response of the units in a given experiment.

We have added a new supplementary figure (Figure 3 —figure supplement 1) in the revision to show firing rates aligned to the turn cue (panels A-B) and the onset of turn movement (panels C-D) for the two groups of units. When aligned to the turn cue, there is no difference in firing rates of Turn Units on Impromptu and Planned Turn trials (page 11). However, when aligned to the turn movement itself there is a decrease in the firing rate of Turn Units on Impromptu Turn trials compared to Planned Turn trials. This difference precedes the turn movement. Therefore, this finding cannot be related to the turn cue but is instead related to the turn movement.

In addition, in calculating the PSTHs displayed relative to the onset of the turn movement, all spikes before movement onset (including those around the instruction and go cues) and after feedback are removed prior to calculating the firing rate (discussed on page 21). Furthermore, the turn onset occurs 541 ± 263 ms and 731 ± 217 ms after the movement onset on Planned Turn and Impromptu Turn trials, respectively. It is unlikely that the movement onset would cause firing rate responses in STN cells >500 ms later. Indeed, responses to the movement onset peak <200 ms after this event. The comparison of firing rates aligned to the turn cue and the onset of the turn movement, the time course of firing rate differences, and our methodology in calculating the PSTHs boosts confidence that our results are not contaminated by the effects of prior events.

Because our goal in this study was to determine neural responses in STN to planned and unplanned movements, our key comparisons are with respect to turn onset not the turn cue. Past human studies in STN frequently used go/no-go or stop-signal tasks with button presses or joystick movements as the operant response. Reporting cue-aligned responses is reasonable in those cases. However, our task has the advantage of recording a full, extended movement, which could also be multi-segmented (Planned and Impromptu Turn trials). We can thus parse when exactly movement onset and turn onset occur. Indeed, this is where we also found significant modulation of firing rates. We now include and discuss the cue-aligned results for Movement Units (page 7) and Turn Units (page 11) in the revised manuscript

Figure 1D: When is the go or the cue signal. Unclear how the kinematic data are arranged in time.

Figure 1F: What is the purpose of clustering trajectories? Is this used in the subsequent analyses?

Figure 1E-1F: Single subject or across all subjects?

Figure 2B/D: How do these trajectories differ? Following feedback, the traces look identical. Were stats performed across all principal components?

Figure 2B/D shows an individual unit classified as a Turn Unit using hierarchical clustering. A key characteristic of this type of unit (the mean of which is displayed in Figure 3C) is a relatively flat response profile except around the turn where Impromptu Turn trials have the lowest firing rate and Planned Turn trials the highest. Other Turn Units have similar response profiles (Figure 2 —figure supplement 1). Averaged across all Turn Units, the firing rate is significantly smaller on Impromptu Turn trials than Planned Turn trials between 240 ms before to 80 ms after the turn. Indeed, for this example unit, the mean firing rates on Impromptu and Planned turn trials over this interval are significantly different (p=0.019, Wilcoxon signed rank test). The similarity of responses following feedback strengthens our classification. Statistics in Figure 2E-H were performed jointly over all principal components, as described on pages 24-25.

Figure 3C: Turn versus impromptu turn changes are < 0.2 Z-scores – these are not convincing changes.

While the difference in normalized firing rates between Planned and Impromptu Turn trials may seem small to the reviewer, it is a substantial range of the overall modulation of firing rates, which varies from -0.16 to +0.18 z-scores (please note that the z-scoring was done based on the dynamics throughout the whole trial, not aligned to a particular event). Further, the difference is highly significant (p=0.0031, linear mixed effects model that controls for unit and subject effects), and quite sufficient for decoding. In fact, a decoder trained to separate trial types based on population responses achieves an accuracy of 74%. The decoder analysis and clarifications above are now included in the revision (page 11).

Use of terms ipsilateral and contralateral is confusing. I think they mean medial/lateral or flexor/extensor as these tasks appear to be occurring in the same hand.

Ipsilateral and contralateral trials and responses correspond to medial and lateral, respectively, since patients always used the hand contralateral from which recordings are made. We have now adopted the medial/lateral terminology suggested by the reviewer.

Figure 7E-F – description says that turn trials show a decreased spike synchronization, but this is because β disappears. Also it appears that this happens regardless of trial type (not just during turns). They might consider also including high frequency LFP activity for comparison with β.

That’s an excellent question. We do not disagree that part of the effect of decreased spike synchronization to β oscillations is because β disappears. However, while β synchronization disappears during the movement phases of the task, firing rates increase. If STN spiking were driven only by β oscillations we would have expected firing rates to decrease at this point. Because the data go against this prediction, STN neurons appear to switch from entrainment to β to alternative firing patterns during the movement phases of the trial including around the turn onset. We also agree that the β desynchronization happens regardless of trial type, and this in itself is noteworthy as, based on existing theories, we would have expected to see an increase in β before impromptu turns. We have clarified these points on page 15 of the revised manuscript.

We had limited our analysis to the lower frequency bands of the spectrum because of significant contamination of the LFP spectrum at higher frequencies by spiking itself. We have now implemented the spike-removal algorithm of Boroujeni et al. (2020) to attempt to analyze this spectrum. The results are shown in Figure 7. We do not detect any difference with the cluster mass test.

Page 24 – Bootstrap is defined as sampling WITH replacement, please justify/clarify.

We had meant to state “permutation procedure”, as we did sample without replacement. We apologize for the error and have now corrected it in the revision (page 27).

The authors should explicitly state what disease the diagnoses of the patients (assuming Parkinson's disease in all, since these are STN implants).

All patients had medically refractory Parkinson’s disease. Their ages and UPDRS scores are listed Supplementary File 1. We regret the omission and have now stated patients’ underlying disease status in various sections of the manuscript, including the abstract.

Supplemental Video with 3D PCA visualization, why are the plotted dimensions PC1, PC3, and PC4? Where is PC2?

In Video 1, our goal was to show the separation of the trajectories in a clearer fashion than Figure 2E-H. For the 3D plot, it was necessary to omit one dimension. Since PC2 distinguishes the different trial types the least, we chose to omit it from this plot. It was still included in all statistical analyses of principal component data. We now clarify this point in the caption for Video 1 in the revised manuscript (page 43) and will be happy to generate additional videos that show PC2.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Essential revisions:

1) Please comment more on what the 'turning' units are actually encoding, and if they could actually be encoding the acceleration and velocity in the vertical axis? Further, do the results still support the statement that 'one sub-population encoding movement kinematics and direction and another encoding unexpected action switches'. Specifically, Figure 4 indicates that both 'Movement Unit' and 'Turn Unit' firing rate can reflect kinematic parameters with different time lags. The authors argued that (Line 276-278): '… the strength of this relationship (for Turn Units) was considerably weaker than for Movement Units (mean R2 = 0.45 compared to 0.61 for Movement Units). Furthermore, movement kinematics did not accurately predict the PSTHs of these units (Figure 4 —figure supplement 2B, D, and F).' However, it is quite subjective to say the R2 values are 'considerably weaker'. Comparing the plots of A, C and E against B, D, and F in Figure 4 —figure supplement 2, can the author's confidently say that movement kinematics can predict Movement Units (shown in plots of A, C and E), but cannot predict Turn Units (shown in plots of B, D and F)? Furthermore, even if the model fitting is a bit worse for Turn Units compare to Movement Units, did the authors used the same model that fits best for movement units? Or to be more fair, did the authors do model comparison for best fitting models for the two types of Units separately?

Thank you for bringing up these important points. The comment raises three key questions, which we answer in the same order asked by the reviewers. First, we explored the possibility of Turn Units encoding vertical acceleration and velocity, using models containing separate predictors for horizontal and vertical movement, as shown in a new figure, Figure 4 – supplement 4. These models did not indicate preferential encoding of vertical movements and performed poorly compared to the models that used the total acceleration and velocity. For example, linear models with separate terms for vertical and horizontal acceleration achieved an out-of-sample R2 of 0.089. Equivalent models for velocity were even worse. Therefore, we do not see evidence that Turn Units preferentially encode vertical acceleration or velocity. These analyses are now explained on page 7 of the revised manuscript.

Regarding the reviewers’ second question about the superior encoding of movement kinematics by the Movement Units, the firing rates of these units were significantly better explained by the movement kinematics compared to the Turn Units. The 95% confidence interval for the out-of-sample R2 of the kinematic model with total speed and acceleration was [0.45 0.51] for the Movement Units and [0.085 0.18] for the Turn Units. The 95% confidence interval for in-sample R2 was [0.57 0.61] for the Movement Units and [0.44 0.49] for the Turn Units. There is no overlap between these confidence intervals, supporting our conclusion. However, please note that our key point in the manuscript is not that Turn Units are void of information about movement kinematics. Rather, we show that:

1. Movement Units represent movement kinematics better while they lack information about unexpected changes of movement plans (unplanned turns).

2. In contrast, Turn Units have distinct firing rates for unplanned and planned turns but weaker representation of movement kinematics.

This relative double dissociation between the Movement and Turn Units is the basis for the names we have chosen for these two groups of neurons. Finally, it is worth clarifying that we arrived at these two groups of neurons not by searching for the neurons that do or don’t represent movement kinematics and unplanned turns. The groups emerged from an unsupervised clustering of the neurons based on the overall dynamics of their PSTHs. We prefer this natural clustering to an engineered search, even though by implementing such a search to identify neurons with specific response properties we could have improved the double dissociation mentioned above. The emergence of the two natural clusters of neurons with different response dynamics and distinct functional properties bolsters the hypothesis that the STN population is inhomogeneous and has the capacity to contribute to complex mental functions. We now clarify the relative double dissociation on pages 4, 11 and 12, and report confidence intervals of the R2 of the kinematic model in Figures 4 – supplements 3-4.

Regarding the reviewers’ third question, we selected the kinematic model with the best out-of-sample R2 for both Movement Units and Turn Units. Figure 4 – supplement 4 shows the model fits for Turn Units. The out-of-sample R2 for Turn Units in the model with total speed and acceleration is 0.13, and the best fitting model, which has total acceleration as its only predictor, has an out-of-sample R2 of 0.23. The same two models produce out-of-sample R2 of 0.48 and 0.43, respectively, for the Movement Units. We now report 95% confidence intervals on the in-sample and out-of-sample R2 of all competing models in Figure 4 – supplements 3-4. The differential ability of the two groups of neurons for encoding kinematics is clear. Further, for any of the top kinematic models, the residual firing rates of Turn Units showed a significant difference for planned and unplanned turns, 240 ms before to 80 ms after turn onset (for the total acceleration model, p=0.0031, HolmBonferroni corrected linear mixed effect model; discussed on page 7). Thus, our conclusions are robust to the exact choice of kinematic models. We now clarify this robustness on page 21.

2) The authors mentioned that 'Turn Units decreased their firing rates before turn onset on Impromptu Turn trials but increased them before turn onset on Planned Turn trials.' (Line 289-291). However, we don't see this statement being supported by data/figure. Can the authors provide evidence/data for this? Or did we miss anything?

Thank you for another great question. In retrospect, we see that our phrasing has been confusing as it could have implied that the firing rates on planned and impromptu turn trials are increased and decreased, respectively, compared to some baseline period. We simply meant that firing rates on impromptu turn trials are lower than on planned turn trials 240 ms before to 80 ms after turn onset (p=0.0031 Holm-Bonferroni corrected linear mixed effect model, test type). The firing rate for the equivalent times on simple reach trials is in between the other two trial types, but it is not a baseline for the statistical tests. We now clarify this point in the paper and have rephrased this sentence (page 7) and a similar statement in our Discussion (page 9-10) to avoid confusion.

3) Reviewer 3 from the previous round of review raised that 'Conceptually, even if 'impromptu' units are an actual subpopulation of neurons, a decrease in STN

activity just before the execution of the unexpected turn seems consistent, not inconsistent, with the previous hypotheses on STN as a "brake" or "global stop" nucleus. We think the authors should clarify their thoughts on this more in their writing.' We feel this was not sufficiently addressed.

Existing theories about the role of STN are based on studies in which subjects were instructed to withhold uninitiated action plans. A common observation in those studies has been that STN responses increase when subjects halt their planned (but uninitiated) movements, leading to the hypothesis that increased activity is associated with stopping (explained in the first paragraph of the Introduction and page 16 in our Discussion). In general, the reviewer is correct that it is possible for any change in STN activity, including a decrease in firing, to be responsible for stopping. However, to our best knowledge, this has not been observed experimentally in past studies. Our finding of decreased STN firing before changes in ongoing movements does not preclude STN acting to stop or change movements, but they do point to the incompleteness of existing models. We now clarify this point to our discussion (page 10).

4) For the statistical analysis Equation. 1 – Equation. 9, it is not clear whether the authors have used self-programmed scripts for these tests or existing toolbox/functions in Matlab to do this. As far as I understand, Equation. 1 – Equation. 9 are similar to what is implemented in the linear mixed effects models available in Matlab (https://www.mathworks.com/help/stats/linear-mixed-effects-models). This can be applied to each time point, and then using cluster mass test to identify the time window with significant effect. Or alternatively, can the authors specify what is difference between the method used in this study and the linear mixed effects model? The authors mentioned that 'For testing factor ;, the levels of ; are permuted within each unit and without permutation of the other factor. The remainder of the calculation continues as for the one factor case, resulting in separate p-values for each main effect, as in a multi-way ANOVA.' Is it equivalent to repeating the multilevel modeling permuting the levels of m? If the method is based on the framework of multi-level modeling, it might not be necessary to present all the Equations.

Our statistical analyses heavily rely on the math of mixed effects modeling, but we could not have readily used pre-existing packages in Matlab because they do not offer a systematic method to determine the time windows for our analyses. To address this shortcoming, we developed the necessary analyses and their associated code by extending the standard cluster mass test. As a reminder, we sought to determine if different conditions (e.g., trial types and movement laterality) were associated with changes in firing rates, but we did not know a priori at which points in time to expect such differences. Furthermore, simple testing of firing rate differences in different windows chosen by examining the data would have resulted in an uncontrolled multiple comparisons problem. Our newly-developed analyses explicitly deal with these multiple comparisons while preserving statistical power. The analysis framework shares many components with the mixed effect models, as astutely observed by the reviewers. Our scripts are available in a public github repository

(https://github.com/dlondon12/InterruptedReachClustMass) and we are working on a methods manuscript that explains the analysis approach in depth.

5) The accurate rate of different conditions was around 0.7. Did the authors only focus on the 'accurate trials' which satisfy the 'accuracy criteria' or considered all trials? Is there any difference between the neural activities for 'accurate trials' vs 'inaccurate trials'?

We did focus all of our analyses on trials which satisfied the accuracy criteria. Those that did not can be split into fixation breaks (which should obviously be excluded as these represent lapses in task engagement) and trials in which the subject performed a movement that did not meet the accuracy criteria. Since our goal was to determine the neural responses to a cue to change movement, we had to ensure that the movement actually changed and in a way that allowed planned and impromptu turns trials to be compared. The median number of inaccurate trials per unit that were not fixation breaks was 19 (aggregated across the three trial types), unfortunately too few to draw conclusions about the neural responses accompanying the inaccurate movements. As expected, adding these trials to our analyses (skipping the exclusion step) does not change our conclusions (e.g., firing rate is significantly smaller on Impromptu Turn trials compared to Planned Turn trials 240 ms before to 80 ms after turn onset). We have now added these points to the Methods section (page 15) and changed the caption for Figure 1B to clarify that it shows the fraction of trials that met the accuracy criteria.

6) What is 'level s' in Equation 9 and Line 748? The ; trial represented trial types, and N represent the different units. What is level s?

Level s refers to the movement direction (target side) of each trial (medial or lateral). This is an orthogonal factor to the trial type (reach, planned turn and impromptu turn trials). The two factor cluster mass test considers the separate effects of each factor. This is now described on page 17 just prior to Equation 4.

7) Can the authors compare the baseline spiking rate of the 'movement related' and 'turning related' units identified by clustering analysis applied on the PSTHs.

The median baseline spiking rate of Movement Units is 35.5 Hz (IQR: 15.5 – 88.0 Hz) compared to 36.6 Hz (IQR: 22.2 – 56.6 Hz) for Turn Units. These rates are not significantly different (p=0.82, Wilcoxon rank sum test). This is now reported on page 4.

8) Figure 2. It is mentioned in the figure legend that 'Shadings in C-H are SE.' The same for Figure 4, 'Shading indicates SE in A-F'. But I didn't see the shadings in these plots? The same for Figure 3 plot B-C, it will be informative to show the SE in the plots.

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

Article and author information

Author details

  1. Dennis London

    1. Center for Neural Science, New York University, New York, United States
    2. Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Visualization, Writing - original draft, Writing - review and editing
    For correspondence
    dennis.london@nyulangone.org
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8134-2683
  2. Arash Fazl

    Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    Contribution
    Conceptualization, Data curation, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  3. Kalman Katlowitz

    1. Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    2. Neuroscience Institute, NYU Langone Health, New York, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5568-6343
  4. Marisol Soula

    1. Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    2. Neuroscience Institute, NYU Langone Health, New York, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  5. Michael H Pourfar

    Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    Contribution
    Investigation, Resources
    Competing interests
    No competing interests declared
  6. Alon Y Mogilner

    Department of Neurosurgery, Center for Neuromodulation, NYU Langone Health, New York, United States
    Contribution
    Investigation, Methodology, Project administration, Resources, Supervision
    Contributed equally with
    Roozbeh Kiani
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1493-0463
  7. Roozbeh Kiani

    1. Center for Neural Science, New York University, New York, United States
    2. Neuroscience Institute, NYU Langone Health, New York, United States
    3. Department of Psychology, New York University, New York, United States
    Contribution
    Conceptualization, Methodology, Project administration, Resources, Software, Supervision, Writing - original draft, Writing - review and editing
    Contributed equally with
    Alon Y Mogilner
    For correspondence
    roozbeh@nyu.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0614-6791

Funding

Simons Collaboration on the Global Brain (542997 543009)

  • Roozbeh Kiani

McKnight Scholar Award

  • Roozbeh Kiani

Pew Scholarship in the Biomedical Sciences

  • Roozbeh Kiani

National Institutes of Mental Health R01 (MH109180-01)

  • Roozbeh Kiani

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

Acknowledgements

The authors would like to thank Gouki Okazawa, and Saleh Esteki for helpful discussions. This work was supported by the Simons Collaboration on the Global Brain (grant 542997), McKnight Scholar Award, Pew Scholarship in the Biomedical Sciences, and National Institute of Mental Health (R01 MH109180).

Ethics

Subjects signed informed consent including consent to publish. The study protocol was approved by the NYU School of Medicine Office of Science and Research Institutional Review Board. Study ID: S16-01855.

Senior Editor

  1. Richard B Ivry, University of California, Berkeley, United States

Reviewing Editor

  1. Nicole C Swann, University of Oregon, United States

Reviewers

  1. Andrew Sharott, University of Oxford, United Kingdom
  2. Harrison Walker

Publication history

  1. Preprint posted: November 12, 2020 (view preprint)
  2. Received: November 13, 2020
  3. Accepted: September 14, 2021
  4. Accepted Manuscript published: September 14, 2021 (version 1)
  5. Version of Record published: October 8, 2021 (version 2)

Copyright

© 2021, London 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

  • 891
    Page views
  • 121
    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)

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

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

Further reading

    1. Neuroscience
    Qing-Tao Meng et al.
    Research Article Updated

    Histamine-dependent and -independent itch is conveyed by parallel peripheral neural pathways that express gastrin-releasing peptide (GRP) and neuromedin B (NMB), respectively, to the spinal cord of mice. B-type natriuretic peptide (BNP) has been proposed to transmit both types of itch via its receptor NPRA encoded by Npr1. However, BNP also binds to its cognate receptor, NPRC encoded by Npr3 with equal potency. Moreover, natriuretic peptides (NP) signal through the Gi-couped inhibitory cGMP pathway that is supposed to inhibit neuronal activity, raising the question of how BNP may transmit itch information. Here, we report that Npr3 expression in laminae I-II of the dorsal horn partially overlaps with NMB receptor (NMBR) that transmits histaminergic itch via Gq-couped PLCβ-Ca2+ signaling pathway. Functional studies indicate that NPRC is required for itch evoked by histamine but not chloroquine (CQ), a nonhistaminergic pruritogen. Importantly, BNP significantly facilitates scratching behaviors mediated by NMB, but not GRP. Consistently, BNP evoked Ca2+ responses in NMBR/NPRC HEK 293 cells and NMBR/NPRC dorsal horn neurons. These results reveal a previously unknown mechanism by which BNP facilitates NMB-encoded itch through a novel NPRC-NMBR cross-signaling in mice. Our studies uncover distinct modes of action for neuropeptides in transmission and modulation of itch in mice.

    1. Cell Biology
    2. Neuroscience
    Sergio Velasco-Aviles et al.
    Research Article

    The class IIa histone deacetylases (HDACs) have pivotal roles in the development of different tissues. Of this family, Schwann cells express Hdac4, 5 and 7 but not Hdac9. Here we show that a transcription factor regulated genetic compensatory mechanism within this family of proteins, blocks negative regulators of myelination ensuring peripheral nerve developmental myelination and remyelination after injury. Thus, when Hdac4 and 5 are knocked-out from Schwann cells in mice, a JUN-dependent mechanism induces the compensatory overexpression of Hdac7 permitting, although with a delay, the formation of the myelin sheath. When Hdac4,5 and 7 are simultaneously removed, the Myocyte-specific enhancer-factor d (MEF2D) binds to the promoter and induces the de novo expression of Hdac9, and although several melanocytic lineage genes are misexpressed and Remak bundle structure is disrupted, myelination proceeds after a long delay. Thus, our data unveil a finely tuned compensatory mechanism within the class IIa Hdac family, coordinated by distinct transcription factors, that guarantees the ability of Schwann cells to myelinate during development and remyelinate after nerve injury.