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

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 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 Figure 1: Subjects completed cued single-and multi-step 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) Accuracy (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. 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.

Diverse Responses of STN Units
We recorded single and multiunit spiking activity from the STN while human subjects performed single-or multi-step 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 ( Fig. 1B-C), with their hand trajectories following the instructed paths ( Fig. 1D-F). However, the latencies of turns in Impromptu Turn trials were systematically longer than movement onset latencies ( Fig. 1B) 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.
We recorded from thirty-nine units in 9 STNs of 8 subjects (1 subject completed separate recording sessions from STNs of both hemispheres). The recorded neurons showed a variety of response dynamics ( Fig. 2A-D), 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 (Fig 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 (Supplementary 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.
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 (Fig. 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 (Fig. 3): units of one group modulated their firing rates with ongoing movement kinematics (Movement Units; 41%), and units of the second group changed their firing rates prior to turn onset, distinguishing different trial types (Turn Units; 46%). Units in the same group were not perfectly homogenous and demonstrated diversity that could suggest sub-clusters within each group. To preserve our four principal components reveals distinct response patterns associated with different task events and trial types. Shadings in C-H are SE. Asterisks indicate Holm-Bonferroni corrected p<0.001 for pairwise permutation tests. 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 Methods).
statistical power, we focus on the two top clusters as a basis for deeper exploration of STN responses in the following sections.

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. Furthermore, firing rates were significantly lower on Planned Turn trials beginning 140 ms after movement onset and lasting until 130 ms before turn onset (Fig. 3B, p=0.0069 and p=0.0064 for Planned Turn vs Reach and Planned Turn vs Impromptu Turn, respectively, Holm-Bonferroni corrected Wilcoxon signed rank tests). 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  (G-H) The firing rate of each unit was modeled as a linear function of time-lagged speed and acceleration (Eq. 10). 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. Whereas Movement Unit responses are well defined by movement kinematics, Turn Unit responses are minimally influenced by them. 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 S1. During the initial fixation period, hands were stationary, but after the go signal, horizontal acceleration and horizontal speed increased on all trial types ( Fig. S1A-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 (Fig, S1C-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 speed) decreased (Fig.  S1E). 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 ( Fig. 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 (Fig. 4E).
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 (Eq. 10). 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 Fig. 4G and Fig. S2A, C, E; R 2 =0.61, averaged 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 ( Fig. 4G and S2G), 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.
To determine if the firing rate differences across trial types were fully explained by movement kinematics, we calculated the residual firing rate around the means that the models above predicted for the movement kinematics of each trial type in a session. Systematic trends in the residuals of different trial types would indicate factors beyond movement kinematics that shape the neural responses. However, such trends were absent in the data. The mean residual firing rates remained indistinguishable from 0 in different task epochs (no significant clusters identified by cluster mass test comparing actual and predicted firing rates), supporting sufficiency of kinematic models for explaining the response dynamics of these units.
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 R 2 for the held-out data 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 (R 2 =0.48 for held-out data). 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 (Fig. S3). 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 Ipsilateral 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 ipsilateral movement direction but similar kinematic-dependent response modulations around the mean for different movement directions. In our task, the instructed movement direction was either ipsilateral or contralateral 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 (Ipsilateral vs. Contralateral). 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 ipsilateral to the recorded STN, as shown in Figure 5A-D, H (p=0.021 Wilcoxon signed rank test).
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 ipsilateral compared to contralateral movements. However, the speed and acceleration were comparable on ipsilateral and contralateral trials (Fig. S4). Further, we used the kinematic model of each unit to predict its firing rate for ipsi and contralateral movements. Residual firing rates around these predictions offer a test for the hypothesis that kinematic differences account for the mean firing rate differences on ipsilateral and contralateral 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 ( Fig. 5E-H) with the mean residual firing rate across this interval greater for ipsilateral trials than contralateral trials (p=0.034, Wilcoxon signed rank test). Therefore, Movement Units showed a sizeable difference in firing rates around movement onset of ipsilateral and contralateral 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 ycomponents of speed and acceleration perform poorly predicting the firing rate? Those models could accommodate differential modulation of firing rates with the kinematics of ipsilateral and contralateral movements. For example, if firing rates changed by different amounts with changes of speed on ipsilateral compared to contralateral 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 ipsilateral and contralateral trials, but the mean firing rate was overall lower on contralateral trials. These differences between ipsilateral and contralateral 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 (Fig. 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.011, Holm-Bonferroni corrected Wilcoxon signed-rank test).
Can this difference be explained by movement kinematics? Like Movement Units, Turn Unit firing rates were also correlated with the lagged kinematics (Fig. 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 ( Fig. 4H and S2H). However, the strength of this relationship was considerably weaker than for Movement Units (mean R 2 =0.45, compared to 0.61 for Movement Units). Furthermore, movement kinematics did not accurately predict the PSTHs of these units (Fig. S2B, D, and F). 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 less than zero on Impromptu Turn trials (p=0.029, Holm-Bonferronicorrected Wilcoxon-signed rank test). Furthermore, the mean residual firing rate over this interval was smaller on Impromptu Turn trials than on Reach or Planned Turn trials (p=0.031 and p=0.022, Holm-Bonferroni corrected Wilcoxon signed rank tests). Therefore, independent of changes in movement kinematics, Turn Units decreased their firing rates before turn onset on Impromptu Turn trials but increased them before turn onset on Planned Turn trials. Unlike Movement Units, Turn Units did not display any difference when comparing ipsilateral and contralateral trials (Fig. S5). 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 (Figs. 5 and S5). Figure 6A-I shows distinct population activity for ipsilateral and contralateral 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 (Fig. 6K, p<0.001, permutation test; also see Supplementary Video 2 for state space trajectories). Thus, ipsilateral and contralateral trials were associated with distinct population responses well before movement onset and for tens of milliseconds following the movement onset.
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 ( Fig. 6J; p=0.024, Holm-Bonferroni corrected permutation test). A similar difference was also present between Planned Turn and Impromptu Turn trials (Fig. 6J, p=0.002, Holm-Bonferroni corrected permutation test 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. 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 ipsilateral and contralateral 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.
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 S1).
Both Movement and Turn Units spiked in synchrony with beta oscillations as quantified by the Point-Phase Consistency (PPC), shown in Figure 7. PPC averaged across all trial types 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 (no significant clusters identified by cluster mass test). (I-J) Average beta band PPC across all trial types in the analysis windows for each task event. Error bars are standard error (* indicates PPC significantly different from 0, p<0.05, Holm-Bonferroni corrected Wilcoxon signed rank test). and times was significantly different from zero for Movement and Turn Units (mean PPC, 0.011±0.006, p=0.0023 for Movement Units, and 0.014±0.005, p=0.0033 for Turn Units, Wilcoxon signed rank tests). This corresponded to a firing rate peak-to-trough modulation of 54% and 63%, respectively, for Movement and Turn Units.
However, contrary to past theoretical predictions, PPC was not significantly different between trial types (no significant clusters detected by cluster mass test, mean PPC across all times not different between trial types, p=0.65 and p=0.68, for Movement and Turn Units, respectively, Friedman's test). Further, the action potentials of Movement Units were uncoupled from the beta oscillation around the time of the turn onset (mean PPC=0.0011±0.0021, Wilcoxon signed-rank test p=0.88). Notably, the uncoupling occurred even earlier for Turn Units, around the time of movement onset (mean PPC=0.0085±0.0043, p=0.064, Wilcoxon signed-rank test; compare Fig. 7I and J). Compatible with these observations, the induced beta power did not differ across trials types (no significant clusters detected by cluster mass test), decreased below baseline shortly before movement onset, and returned to baseline after feedback ( Figure S6). 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. These observations 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 are only weakly related to movement kinematics; they primarily respond to an unpredictable need to switch movement trajectory.
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 population activity decreased beginning 240 ms before the switch in movement trajectory on Impromptu Turn trials, in contrast to the increased population activity in the same time window around the turn onset on Planned Turn trials. Further, 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 been commonly associated with increased activity in the STN population (Bastin et al., 2014;Isoda and Hikosaka, 2008;Pasquereau and Turner, 2017;Schmidt et al., 2013). 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;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 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 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 modulaition 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 (Figs. 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 reservoire 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 ipsilateral and contralateral movements (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 ipsi and contralateral 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 (Fig. 6). Within the STN population, the mean residual activity of Movement Units after removal of kinematic modulations distinguished ipsi and contralateral movement directions until 270 ms after Movement Onset (Fig. 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 (Figs. 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. Interestingly, these sub-populations are partially spatially segregated within STN, with Movement Units located more laterally, consistent with a more prominent role for dorsolateral STN in movement circuitry compared to ventromedial STN (Hamani et al., 2017). 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 plan 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). Further, stimulation of glutamatergic cortical inputs to STN (i.e., the hyperdirect pathway) worsens bradykinesia (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 Pasquerau and Turner note, this richness in neural code appears incompatible with a narrowly defined role for STN that includes only movement inhibition (Pasquereau and Turner, 2017). 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.

Lead Contact
Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Dennis London (dennis.london@nyulangone.org).

Materials Availability
This study did not generate any new unique reagents or other materials.

Data and Code Availability
The dataset and custom code generated in this article are available from the Corresponding Author upon request.

Subjects
Subjects were patients undergoing implantation of deep brain stimulation electrodes in the subthalamic nucleus. 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 intraoperative session. All experimental procedures were approved by the Institutional Review Board at NYU Langone Medical Center.

Surgery and Electrophysiological Recordings
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. Preoperative localization of the subthalamic nucleus was guided by standard anatomical landmarks on MRI, including the craniocaudal midpoint of the red nucleus and the indirect anterior commissure -posterior commissure (ACPC) based coordinates, roughly 11 mm lateral to midline, 5 mm inferior to the ACPC 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 5 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 microelectode 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. Entrance into the STN was determined by an increase in the high frequency background electrophysiological activity (Novak et al., 2007) and the presence of rhythmically firing cells with high firing rates (Hutchison et al., 1998) as assessed by a trained neurophysiologist (author M.P.).

Behavioral Task
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 towards 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). 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 seconds.
A session was terminated when the subject showed signs of fatigue or after approximately 10 minutes. 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
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 principle components of spike waveforms. Single units, multi-units, and unsorted background spiking (hash) were considered as separate units. 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 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. On Reach trials, subjects did not make a turn, and therefore, there was no turn onset event. To compare responses around turn onset on turn trial types with the 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. Interevent 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 from the end of the preceding trial were also censored. For smoothing the PSTHs, we used a causal kernel (alpha function, α ! τe -!" for > 0 and 0 for τ ≤ 0, withα=20 s -1 ). PSTH error bars were calculated as the standard error over trials. Units were excluded from analyses if there were fewer than 4 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). For analysis of PSTHs on contralateral and ipsilateral trials, we only examined units with at least 4 correct ipsilateral and contralateral trials of each type (mean±s.d., 13.5±3.3 trials per trial type per subject).
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
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 the Wilcoxon signed-rank test. 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 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 signficacnce. 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 , consider a firing rate response profile, r !,! t , where m represents the trial type. Let's denote the mean response across units with ! , across trial types with ! , and across units and trial types with ( ). The total sum of squares, !"! , the sum of squares across the trial types, ! , and across the the units, ! , are: Considering each unit as a repeated measure, the residual sum of squares is Eq. 2 and the statistic for the difference across trial types is This statistic was thresholded at the 90 th percentile of the F-distribution with parameters − 1 and − 1 − 1 . 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 !""#" within each cluster in this paper, but we also obtained similar results with other relevant statistics including the sum of ! or of !!!, !!! !!! . The null hypothesis is that the 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 -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 , with 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, !,!,! ( ). For factor , the sum of squares is For factors and , we must also consider the interaction sum of squares: where !,! is the mean across trials for the factors and . Analogous expressions can be constructed for each individual factor and each pair of factors. -statistics are calculated for each main effect: 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 -statistic calculation. For each tested effect, there is a separate permutation distribution depending on the "exchangeable units" (Anderson, 2001). 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.
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 units which represent repeated measures. Principal components analysis of these units results in dimensionality reduction to 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 !,!,! , we consider -dimensional data of the form !,! = [ !,!,! , … , !,!,! ], where !,! ( ) is the projection of !,!,! onto the first principal components, and !,!,! ( ) is the projection onto the -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 ! , where each time point is the centroid across the levels of s, and , 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: where . is the Euclidean norm.
We threshold these sums of squares above the 90 th percentile as was done with the Fstatistics. To determine the value of the 90 th 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, 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 ( Fig. 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).

Mixed Effect Models
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 ( ! , ! ) on the screen. The cursor position was filtered using a Gaussian kernel (10 ms standard deviation) and used to calculate instantaneous speed and acceleration time : where is the time interval between consecutive measurements. We then calculated the magnitude of the instantaneous speed and acceleration vectors on the screen: where . is the Euclidean norm. To combine data across sessions, we normalized ! ( ) and ! ( ) by z-scoring across all trials in each session. Figure 4 shows averages of normalized speed and acceleration for each trial type, , 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: where ! is the mean firing rate of the unit at time of trial type , ! and ! are the mean movement speed and acceleration, respectively, and ∆ is the time lag. The parameters were fit as separate fixed effects for each unit. We systematically varied ∆ , searching for the model with the highest R 2 as the best fitting model. Across units and across ∆ that were modeled, the firing rate time points that were fitted were kept constant. The maximum positive and negative ∆ determined which data points were fitted and which were held out for model testing. Consider an ! on the range [ !"#$" , !"# ] and ∆ modeled on the range [ !"# , !"# ]. There is a subset of points in [ !"#$" , !"# ], such that there is data for ! + ∆ and ! + ∆ (i.e. + ∆ is also on the interval [ !"#$" , !"# ]). Therefore, ! was only fit on the range [ !"#$" − !"# , !"# − !"# ]. The optimal value of ∆ 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 [ !"#$" − , !"# − ] (note this range contains the best fit ∆ ) were used as a testing set for the models. Eq. 10 was chosen as a good descriptive model of firing rates after extensive search of possible linear models of speed and acceleration.
Overall, we compared 14 alternative models that contained a combination of speed, acceleration, and their interaction (Fig. S3). 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

LFP and Joint LFP-spiking analysis
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 seconds before to 2.5 seconds 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 first 250 milliseconds after the fixation onset to calculate a baseline for the average evoked and induced power for each trial. We then used a bootstrap 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. Repeaing 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 ! 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, 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 seconds, 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 spiketriggered 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 event-aligned analysis intervals were limited to at most 500 ms, we limited our frequency-based analysis to above 4 Hz.

Movement Trajectory Analysis
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 { ! } is: which can be rewritten as the mean of the mean of independently-sampled pairs of data points, and : 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: 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. Figure S2: Firing rates were regressed on movement acceleration and speed. (A-F) The actual firing rates (light solid lines) and the model fits (dark dotted lines) are shown for Reach (A, B), Planned Turn (C, D), and Impromptu Turn (E, F) trials. Regressions were performed on a fraction of data; out of sample firing rate predictions are shown with dark xmarked lines. The fits were performed using time-lagged acceleration with the lags optimized to achieve the best fit. The adjusted R 2 of the fits with different time lags for different units are shown for Movement Units (G) and Turn Units (H). The highest adjusted R 2 for each unit determined the best-fit time lags in Fig. 4G, H. Figure S3: Total speed and acceleration provide a better explanation of the Movement Unit responses compared to more complex models with distinct parameters for different movement directions. (A-P) Titles of each panel indicate the parameters of each model. and indicate speed and acceleration, respectively. x and y subscripts indicate movement along the horizontal and vertical axes, respectively. |. | indicates Euclidean norm, which equals the absolute value for onedimensional variables (e.g., ) and the vector magnitude for two-dimensional variables (e.g., total speed). The models were trained for time points shared across all lags and tested held-out times not included in training. Training and testing R 2 are displayed for each model ( for in-sample, training data; for out-of-sample, test data). Negative indicate residual variance that exceeds the variance of underlying data, indicating extremely poor model performance for test trials, likely due to overfitting the training data. (Q) The best-fitting time lags for each unit were similar across models. Columns correspond to models in panels A-P. Each row corresponds to a specific Movement Unit with units ordered by time lag of our chosen model (the best-fitting model containing total speed and acceleration, shown in panel A).

Figure S4
Speed and acceleration are similar on ipsilateral (lighter color) compared to contralateral trials. Speeds are plotted in the left two columns and accelerations in the right two columns. Shading is SE. Figure S5: There were no significant differences between firing rates on ipsilateral and contralateral trials on Turn Units (no significant clusters identified by cluster mass test). Figure panels are data from Turn Units arranged analogously to

Figure S6
: Induced power spectrograms (A-C) and induced power in the beta range (D) aligned to task events are shown for different trial types. Induced power decreases from baseline beginning before movement onset and lasting until after feedback. There were no significant differences between trial types (no significant clusters identified by cluster mass test).