## Abstract

Motor cortex plays a substantial role in driving movement, yet the details underlying this control remain unresolved. We analyzed the extent to which movement-related information could be extracted from single-trial motor cortical activity recorded while monkeys performed center-out reaching. Using information theoretic techniques, we found that single units carry relatively little speed-related information compared with direction-related information. This result is not mitigated at the population level: simultaneously recorded population activity predicted speed with significantly lower accuracy relative to direction predictions. Furthermore, a unit-dropping analysis revealed that speed accuracy would likely remain lower than direction accuracy, even given larger populations. These results suggest that the instantaneous details of single-trial movement speed are difficult to extract using commonly assumed coding schemes. This apparent paucity of speed information takes particular importance in the context of brain-machine interfaces (BMIs), which rely on extracting kinematic information from motor cortex. Previous studies have highlighted subjects' difficulties in holding a BMI cursor stable at targets. These studies, along with our finding of relatively little speed information in motor cortex, inspired a speed-dampening Kalman filter (SDKF) that automatically slows the cursor upon detecting changes in decoded movement direction. Effectively, SDKF enhances speed control by using prevalent directional signals, rather than requiring speed to be directly decoded from neural activity. SDKF improved success rates by a factor of 1.7 relative to a standard Kalman filter in a closed-loop BMI task requiring stable stops at targets. BMI systems enabling stable stops will be more effective and user-friendly when translated into clinical applications.

- motor control
- neural coding
- brain-machine interface

previous studies have investigated the extent to which motor cortex encodes kinematic variables, including movement direction (Ashe and Georgopoulos 1994; Georgopoulos et al. 1982; Schwartz et al. 1988) and movement speed (Churchland et al. 2006; Moran and Schwartz 1999b; Schwartz 1992, 1994). The ability to accurately read out direction and speed from motor cortex takes particular importance in the context of brain-machine interfaces (BMIs), which translate neural activity into control signals for driving prosthetic devices, such as robotic limbs (Carmena et al. 2003; Chapin et al. 1999; Velliste et al. 2008; Wessberg et al. 2000), or computer cursors (Gilja et al. 2012; Mulliken et al. 2008; Serruya et al. 2002; Suminski et al. 2010; Taylor et al. 2002). Despite impressive advances in BMI technologies in recent years, BMI control of robotic limbs and cursors is still inferior to able-bodied control of natural limbs and physical pointing devices, especially with respect to the stability of stopping, as pointed out in previous studies (Carmena et al. 2003; Ganguly and Carmena 2009; Gilja et al. 2012; Hochberg et al. 2006; Kim et al. 2008). To better understand the origin of this poor control of BMI movement speed, we looked for signatures of a robust representation of instantaneous movement speed in single-trial reaching movements.

We analyzed spike trains recorded simultaneously across primary and premotor cortices of rhesus monkeys during a three-dimensional (3D) center-out reaching task. Using standard information theoretic and population decoding techniques, we found substantially less speed-related information than direction-related information in neural activity at the levels of both single units and simultaneously recorded populations. We also performed a unit-dropping analysis, which suggests that our ability to decode movement speed might not improve substantially with access to larger numbers of neurons. None of our analyses revealed a substantial representation of the moment-by-moment details of movement speed on single-trial bases.

The finding that speed information is difficult to extract from motor cortical population activity informed a novel approach to implementing movement speed when driving BMI devices. This decoder, termed the speed-dampening Kalman filter (SDKF), incorporates the assumption that movement speed and angular velocity should be inversely related. Rather than relying on neural activity to provide the complete details of movement speed, which may be difficult to extract in the real-time setting of BMI, SDKF enhances control of movement speed by using directional signals, which are more easily extracted from neural activity. Since movement direction can be reliably inferred from motor cortical population responses, angular velocity (the temporal derivative of direction) can be extracted reliably as well. SDKF uses angular velocity of the decoded cursor trajectory to modulate cursor speed, thus reducing the system's reliance on cortical activity to directly provide the moment-by-moment details of movement speed.

## METHODS

### Neural Recordings and Behavioral Tasks

Monkeys performed two tasks: an arm reaching task and a BMI cursor control task. All animal procedures were performed with the approval of the Institutional Animal Care and Use Committee of the University of Pittsburgh.

#### Neural recordings I: arm reaching.

Two male rhesus macaques (*Macaca mulatta*) were implanted with 96-channel microelectrode arrays (Blackrock Systems, Salt Lake City, UT) in motor cortex contralateral to the reaching arm used in the behavioral task. Neuronal activities were manually sorted (Plexon, Dallas, TX) from single- and multineuron units, and spike times were recorded throughout the behavioral tasks. *Monkey F* arm reaching data have been previously described in Fraser and Schwartz (2012). Briefly, *monkey F* had two arrays: one array targeted proximal arm area of primary motor cortex, and a second array targeted ventral premotor cortex. Across both arrays, 119 units were identified and tracked across 4 experimental sessions using techniques described in Fraser and Schwartz (2012). *Monkey T* had a single array targeting arm proximal area of primary motor cortex. Across 5 experimental sessions, 67.8 ± 11.4 units were identified.

#### Behavioral task I: arm reaching.

Both monkeys were trained to perform 3D arm reaching movements. Arm movements were tracked at 60 Hz using an infrared marker (Northern Digital, Waterloo, ON, Canada) taped to the wrist of the reaching hand. Tracked positions were displayed to the subjects as a spherical virtual cursor (radius: 8 mm) on a stereoscopic display (Dimension Technologies, Rochester, NY). Movements were either from a workspace-centered virtual target to 1 of 26 virtual peripheral targets (center-out) or from a peripheral target to the central target (out-center). All targets were displayed as spheres (radius: 8 mm), and peripheral targets were distributed roughly evenly about the surface of a virtual sphere (radius: 66 mm, *monkey F*; 75 mm, *monkey T*). A trial was initiated by the subject acquiring visible overlap of the cursor with the start target for 400–600 ms. Next, a virtual target was presented, and the subject was required to acquire that target with the virtual cursor within 800 ms of presentation and hold with visible overlap for another 400–600 ms. Trials were deemed successful upon completion of this sequence and were followed by a water reward of 60 μl (*monkey F*) or 150–190 μl (*monkey T*). Failed trials were not rewarded.

We analyzed 1,040 successful trials from *monkey F* and 1,316 successful trials from *monkey T*. All analyses were performed on data recorded during the period between completion of the start hold and beginning of the target hold.

#### Neural recordings and behavioral task II: BMI control.

*Monkey F* also performed an eight-target 2D center-out BMI cursor task, whereby recorded neural activity was translated in real-time into movements of a BMI cursor. The cursor (radius: 7 mm) and targets (radius: 7 mm) were displayed to the subject on a frontoparallel display. Target directions were chosen pseudorandomly from one of eight directions spaced uniformly about the perimeter of a workspace-centered circle (radius: 85 mm), and unless noted otherwise, target hold times were randomly drawn from a uniform distribution (range: 0–600 ms). The subject initiated a new trial by modulating neural activity to drive the cursor to visibly overlap a workspace-centered target for 150 ms. After this initial hold, a peripheral target appeared, instructing the subject to acquire the target with the BMI cursor. A trial was deemed successful if the subject acquired and maintained target acquisition for the trial-specific hold period. A trial was deemed failed if the cursor left the target within the hold period following target acquisition, or if the target was not acquired within 3 s after target onset. The subject was naive to each trial's target hold requirement until trial completion. Successful trials were rewarded with 150–180 μl of water. To initiate the next trial, the subject needed to return the cursor to again visibly overlap a workspace-centered target for 150 ms. The cursor was automatically returned to the workspace center only following trials that the subject failed by exceeding the 3-s time limit on target acquisition.

Two-dimensional cursor velocity was decoded from recorded neural population activity using either a novel speed-dampening Kalman filter (SDKF) or a standard velocity-only Kalman filter (VKF). Each experimental session consisted of alternating blocks of trials under SDKF control and blocks of trials under VKF control. The decoder applied during the first block was selected randomly at the beginning of each session. Across 6 experimental sessions, neural responses from 86.1 ± 12.2 single and multiunits were sorted, and spike counts were recorded in 33-ms nonoverlapping bins. In total, the subject performed 1,216 successful trials with each decoder.

During four additional experimental sessions, target hold requirements were relaxed to 50 ms. In these sessions, cursor movements were decoded by VKF only, using 95.3 ± 9.7 units. All other experimental methods match those described above with the exception that in two of these sessions, the cursor automatically snapped back to the workspace center following trial success. The subject performed 2,352 trials under this 50-ms hold condition.

### Data Discretization

Movement speed and direction are continuous-valued quantities expressed using different units and numbers of degrees of freedom. To enable an unbiased comparison of the relationships between neural activity and these two kinematic quantities, we discretized movement speed and direction such that their statistical properties were matched. Arm movements were segmented into nonoverlapping 30-ms intervals, and average movement speed and movement direction were computed across each 30-ms interval. We labeled each 30-ms interval with 1 of 26 candidate speed labels and 1 of 26 candidate direction labels. For a given movement speed, the applied speed label corresponded to the nearest of 26 candidate speed centroids. This set of speed centroids was chosen for each experimental session such that each label was applied to approximately the same number of data points (Fig. 1*B*). For a given movement direction, the applied direction label corresponded to the direction centroid, of 26 candidate centroids, whose angle with the actual movement direction was smallest. Direction centroids were optimized such that each of the 26 direction labels was applied to approximately the same number of data points (Fig. 1*C*). The details of this optimization procedure are provided in the Appendix. This discretization procedure results in a uniform prior distribution of movement speeds and a matched uniform prior distribution of movement directions. By matching these distributions, chance prediction accuracy is thus also matched between the two kinematic quantities.

### Information Analysis

To characterize the speed-related and direction-related information carried by spike trains from individual neural units, we computed the mutual information between neural activity and movement kinematics. Mutual information is a statistic describing the extent of dependence between two random variables. In contrast to correlation coefficients determined from linear regression, mutual information does not require specifying a form for the relationship between two variables, allowing it to capture arbitrary nonlinear relationships if they are indeed present in the data.

Over a series of time lags, τ ∈ [−300 ms, 300 ms], we assembled tuples of kinematic and neural measurements, {*s*_{t}, *d*_{t}, **y**_{t−τ}}, where *s*_{t} ∈ {*s*_{1}, . . . , *s*_{26}} is the discretized speed for the 30-ms interval beginning at time *t*, *d*_{t} ∈ {*d*_{1}, . . . , *d*_{26}} is the discretized direction for the 30-ms interval beginning at time *t*, and **y**_{t−τ} ∈ ℤ^{q} is a vector containing the spike counts across the *q* units during the 30-ms interval beginning at time *t* − τ. For unit *j* and time lag τ, the mutual information between spike counts and discretized kinematics is (1)
where **Y**_{j} is the set of unique, nonoverlapping 30-ms spike counts observed for unit *j* (e.g., **Y**_{j} = {0, 1, 2} for a neuron that spiked at most twice during any 30-ms bin), and **X**_{τ} is the set of kinematics labels for time intervals that lag the spike counts in **Y**_{j} by τ ms. The labels in **X**_{τ} correspond to discretized movement speeds, {*s*_{1}, . . . , *s*_{26}}, when movement speed is the kinematics variable of interest, and similarly, to discretized movement directions, {*d*_{1}, . . . , *d*_{26}}, when movement direction is the kinematics variable of interest. The terms *p*(*x*), *p*(*y*), and *p*(*x*, *y*) are the normalized frequencies of kinematics label *x*, spike count *y*, and joint pair (*x*, *y*), respectively. In the event that the pair (*x*, *y*) does not appear in a data set [i.e., *p*(*x*, *y*) = 0], we evaluate the summand in *Eq. 1* to be 0.

For each unit we determined the lag at which mutual information was maximized between spike counts and lagged direction, τ_{direction}, and lagged speed, τ_{speed}, where positive lags correspond to causal relationships between neural activity and movement kinematics. Henceforth, we refer to the mutual information at these optimal lags as maximal direction information (MDI) and maximal speed information (MSI).

#### Significance testing for information analysis.

Information measures are known to be biased such that one can measure positive values of information when in fact the variables are independent (Treves and Panzeri 1995). To determine whether a unit's MDI and MSI values were greater than expected by chance, we performed the following permutation test. To determine a null information value, we shuffled the correspondence between spike counts and discretized kinematics and then computed mutual information. This shuffle preserves the marginal distributions of spike counts and kinematics but destroys any relationship between the quantities. To obtain a distribution of null information values, we repeated this procedure 10,000 times using speed as the kinematic variable and, similarly, using direction as the kinematic variable. We determined a *P* value for each kinematic variable to be the fraction of null information values that were larger than the single mutual information value determined from the nonshuffled data.

To determine whether a unit's MDI and MSI values were significantly different from each other, we performed the following bootstrap procedure (Efron and Tibshirani 1993). For a data set consisting of *n* timesteps, we generated a resampled data set by randomly drawing *n* timesteps with replacement from the original data set and then computed mutual information values across all lags for speed and direction. We repeated this computation 20,000 times and determined *P* values by computing the fraction of resampled computations resulting in MDI greater than MSI or MSI greater than MDI.

### Regression Analysis

To establish a link between this information analysis motor cortical tuning, we also performed a linear regression analysis (Ashe and Georgopoulos 1994; Georgopoulos et al. 1982; Lebedev et al. 2005; Perel et al. 2013; Schwartz 1992). We fit the following direction-only, speed-only, and velocity tuning models:

Direction-only tuning:
(2)
Speed-only tuning:
(3)
Velocity tuning:
(4)
where *y*_{t} is the spike count during timestep *t*, **v**_{t} = [*v*_{t,1} *v*_{t,2} *v*_{t,3}] is a 3D reach velocity, and ‖**v**_{t}‖ is the corresponding reach speed. The {*b*} are coefficients fit to data. Each model was fit separately across a range of time lags, τ ∈ [−300 ms, 300 ms].

### Neural Decoding for Arm Reaching

To characterize the kinematic information carried by simultaneously recorded population activity, we performed a population-decoding analysis. We trained Poisson naive Bayes (PNB) classifiers (Shenoy et al. 2003) to independently predict discretized movement speed and movement direction from a 30-ms population spike count vector. In the current analysis, PNB assumes that *1*) each neuron fires at a characteristic rate determined by the current kinematics (either movement direction or movement speed), *2*) given these kinematics, each neuron fires independently, and *3*) observed spike counts are Poisson noise-corrupted instantiations of the characteristic rates. PNB, while explicitly specifying the structure of the relationship between neural activity and kinematics, can capture nonlinear tuning effects and Poisson-like signal-dependent noise.

The probabilistic model for PNB is given by (5) (6) (7)
where *Eq. 5* defines the prior probability of kinematics *x*_{t}, *Eq. 6* is the probability of the observed spike count *y*_{t,j} for unit *j* given the current kinematics, and *Eq. 7* is the probability of the observed population spike count vector, **y**_{t} = [*y*_{t,1}, . . . , *y*_{t,q}]′, across q simultaneously recorded units given the current kinematics. The parameters of the PNB model are the *p*_{k} for *k* ∈ {1, . . . , 26}, representing the prior probability of kinematics label *k*, and the firing rate parameters {γ_{j,k}} for each neuron *j* given kinematics label *k*. These model parameters were determined via maximum likelihood over the training data. By design of the kinematics discretization, *p*_{k} ≈ 1/26 for both speed and direction, where the correspondence is approximate rather than exact due to the fact that training data were chosen randomly from each data set (see description of cross-validation below). To predict movement kinematics given an observed spike count vector, we compute (8)
Spike counts and discretized kinematics used in this analysis were identical to those used in the information analysis.

As described, PNB enables a prediction of the current kinematics given the current neural activity. To predict kinematics based on a history of neural activity, we used an augmented PNB model that incorporates the entire causal history of nonoverlapping 30-ms spike counts beginning 300 ms before the corresponding movement kinematics. The probabilistic model for this history-based PNB includes the prior from *Eq. 5* and replaces *Eqs. 6* and *7* with (9) (10)
where *Eq. 9* gives the probability of having observed *y*_{t−τ,j} spikes from unit *j* at time *t* − τ, given the current kinematics, and *Eq. 10* is the joint probability of having observed the history of spike counts over the past 300 ms across the *q* recorded units, given the current kinematics. The parameters {γ_{j,k,τ}} are now indexed for each neuron *j*, kinematics label *k*, and time lag τ ∈ *T*, where *T* is the set of lags from 0 to 300 ms in 30-ms intervals. To predict movement kinematics given the spike count history, we compute (11)

To assess how well the speed- and direction-based PNB models would generalize to unseen data, we performed twofold cross-validation. Data were randomly partitioned into two subsets. First, we trained PNB models on the first subset and evaluated predictions using the second subset. Next, we reversed this process, training PNB models on the second subset and evaluating predictions on the first subset. In this fashion, each model was evaluated using data not seen during model fitting.

### Simulated Neural Populations

To provide intuition for the information and prediction analyses, we simulated neural population activity under four parametric encoding models. For each model, we fit parameters to actual recorded neural activity and arm kinematics from a representative data set (F081909). We then generated spike counts from each encoding model, again using the actual recorded kinematics from the example data set. This procedure provides simulated data sets that exactly match the example real data set with respect to the number of units, number of reaches, duration of reaches, and statistics of kinematics.

We simulated positive firing rates according the following encoding models.

Direction-only tuning:
(12)
Speed-only tuning:(13)
Velocity tuning:
(14)
Independent speed and direction:
(15)
where *r*_{t} is the firing rate at timestep *t*, **v**_{t} = [*v*_{t,1} *v*_{t,2} *v*_{t,3}] is a 3D reach velocity, and ‖**v**_{t}‖ is the corresponding reach speed. The {*b*} are coefficients fit to data. For a given unit, the time lags, τ, in the direction-only and speed-only models were chosen to be the lags at which that unit achieved its MDI and MSI, respectively, from the information analysis over the real data. In the velocity and independent speed and direction models, the time lags were chosen to be the lag associated with the larger quantity between that unit's MDI and MSI. After simulating these firing rates, we then generated noisy spike counts, *y*_{t}, according to (16)
where Δ*t* = 30 ms matched the binning used in the information and decoding analyses.

### Neural Decoding for BMI Control

Two-dimensional cursor velocity was decoded from binned spike counts using either a VKF or a novel SDKF. For both decoders, 2D cursor positions were computed by integrating the corresponding decoded velocity.

#### Velocity-only Kalman filter.

For BMI control, we implemented a Kalman filter (Kalman 1960) to predict intended movement velocity given a sequence of recorded neural activity. Kalman filter predictions combine knowledge from a trajectory model describing the relationship between velocities from one timestep to the next, and from an encoding model describing the relationship between spike counts and intended velocity. When trajectory and encoding models are linear-Gaussian, the Kalman filter velocity predictions are optimal with respect to mean square error of predicted velocities.

The trajectory model underlying the Kalman filter takes the form (17)
where **v**_{t} ∈ ℝ^{2} is the velocity intention at timestep *t*, **A** ∈ ℝ^{2×2} maps beliefs about the velocity at timestep *t* − 1 into beliefs for timestep *t*, and **Q** ∈ ℝ^{2×2} is a covariance matrix describing the uncertainty corresponding to this mapping. *N* denotes a Gaussian (normal) distribution. The encoding model is defined as (18)
where **y**_{t} ∈ ℝ^{q} is the vector of spike counts simultaneously recorded across *q* units at timestep *t*, **C** ∈ ℝ^{q×2} maps intended velocity to expected spike counts, **d** ∈ ℝ^{q} accounts for baseline firing rates, and **R** ∈ ℝ^{q×q} is the observation noise covariance. We fixed **A** = **I** and estimated **Q**, **C**, **d**, and **R** (constrained to be diagonal) via linear regression over data collected from a calibration session (described below). In every session, these estimated parameters resulted in a stable VKF decoder. The details of the VKF algorithm and its stability are provided in the Appendix.

Kalman filters have been applied effectively toward decoding movement kinematics from neural activity in the context of both offline reconstruction of natural arm movements (Wu et al. 2006) and online control of a BMI (Gilja et al. 2012; Hochberg et al. 2012; Kim et al. 2008; Koyama et al. 2010; Li et al. 2011; Orsborn et al. 2012; Wu et al. 2004). When implementing a Kalman filter, one must select the state variables to be modeled by the trajectory and encoding models. Previous studies have shown that Kalman filters with a velocity-only state representation provide superior online BMI control compared with a position-only state representation (Kim et al. 2008). Thus, as a baseline for comparison, we implemented a velocity-only Kalman filter (VKF).

Training data for building the VKF decoder were collected during a closed-loop calibration session prior to each experiment (Chase et al. 2012; Velliste et al. 2008). Calibration sessions consisted of center-out trials with decreasing levels of assistance, whereby cursor velocities orthogonal to the center-to-target direction were automatically attenuated. In an initial block of eight trials, model parameters were chosen randomly and complete error reduction was applied, resulting in straight-to-target cursor trajectories. VKF parameters were fit to the recorded neural activity and velocity intentions, which were assumed to be in the center-to-target direction with constant speed. For this initial fitting step, each trial contributed roughly 30 timesteps of both intended velocity and spike count vectors (recorded in 33-ms nonoverlapping bins). This initial quantity of data appears to be sufficient for determining an initial set of VKF parameters. In the second block of eight trials, these VKF parameters were implemented and error attenuation was decreased slightly. We repeated this cycle for typically five blocks (40 trials), fitting new VKF parameters after each block using all previous trials. All error attenuation was eliminated by the last calibration block such that the subject operated the BMI under complete neural control.

#### A speed-dampening Kalman filter for closed-loop BMI control.

The SDKF extends the VKF by enforcing a tradeoff between movement speed and magnitude of angular velocity. When cursor trajectories exhibit large absolute angular velocities, SDKF constrains decoded speeds to be closer to 0 in a graded fashion depending on the magnitude of the angular velocity. SDKF implements this tradeoff through an adaptive trajectory model (19)
where λ_{t} ∈ [0, 1] is a time-varying speed-dampening factor that is given values near 1 when the cursor trajectory has been straight and shrinks toward 0 as angular velocity increases. The details of the SDKF algorithm are provided in the Appendix. In a given experimental session, the remaining decoding parameters for SDKF, **A**, **Q**, **C**, **d**, and **R**, were identical to those used for VKF.

### Simulated Closed-Loop Control of Movement

To establish a link between the information analyses and the closed-loop BMI control experiments, we simulated closed-loop BMI control. Two-dimensional control of a cursor was driven by a simulated population of neurons with log-linear tuning curves parameterized by the independent speed and direction model of *Eq. 15*. Parameters of these tuning curves were fit to neural and kinematics data recorded in the arm reaching task (data set F081909) between target onset and target acquisition. Speed and direction data were not discretized, and movement directions were truncated from 3D to 2D to match the 2D BMI task. Simulated BMI movements were decoded using VKF, which was trained on these same arm reaching kinematics and corresponding simulated spike counts. The simulated task was matched to the real BMI behavioral paradigm with respect to target positions, target hold requirements, and conditions for task success.

At each simulated timestep, desired kinematics were chosen based on target position and the most recent cursor position. Desired movement direction was straight from the most recent simulated cursor position to the target position. Desired movement speed depended on the distance between the target and the most recent simulated cursor position. Desired speed was 0 if the cursor and target visibly overlapped by at least one-half of the cursor radius. Otherwise, desired speed was drawn from a normal distribution whose mean and standard deviation were matched to real arm movement data for similar cursor-to-target distances. Simulated spike counts were drawn from Poisson distributions with rates determined by these desired kinematics and the log-linear tuning curves. Finally, a cursor update was decoded from the simulated spike counts using VKF.

## RESULTS

Here we present our findings from the information, regression, and decoding analyses over the arm reaching data. The results consistently suggest that movement speed is substantially more difficult than movement direction to extract from the moment-by-moment details of neural firing in motor cortex. Next, we link these findings to BMI control, which we show to suffer from inadequate speed control, especially with respect to stopping, relative to direction control. Finally, we demonstrate improved BMI stopping and speed control under a novel BMI decoder, the speed-dampening Kalman filter (SDKF).

### Single-Unit Activity Carries More Information About Direction Than About Speed

Spike trains from single units contained speed- and direction-related information in a variety of forms. Figure 2 shows speed and direction information as a function of lag for a number of representative units. Direction information curves were unimodal for nearly every recorded unit, whereas it was not uncommon for speed information curves to be bimodal. These bimodal speed information curves are likely a reflection of task-induced autocorrelation in movement speed (e.g., through bell-shaped speed profiles). Most units' MDI and MSI values were significant relative to null information levels, although some units were exceptions, as detailed in Table 1. The number of units with significant MDI was larger than the number of units with significant MSI for all nine experiments.

Figure 3 shows the lags at which each unit achieved its MDI and MSI. Optimal lags for direction information were most often causal, meaning information-carrying spikes tended to lead movement direction. Optimal lags for speed-related information were casual and acausal in roughly equal frequencies. Many units thus had substantial discrepancies between their optimal direction and speed lags, and, on average, optimal direction lags were more positive (i.e., causal) than optimal speed lags within individual units.

Perhaps most striking, however, were the differences between MDI and MSI values within individual units. Of 119 units from a representative data set (F081809), 45 (38%) had MDI values that were significantly greater than their MSI values, whereas only 12 (10%) showed the opposite relation (Fig. 4*A*). Consistent with this breakdown, we found significantly more direction-related information than speed-related information on average across all recorded units with significant differences in MDI and MSI from this data set (Fig. 4*B*). This breakdown of unit types was consistent across data sets from both subjects, with 3.01 ± 0.63 times as many direction-dominated cells as speed-dominated cells (see Table 2). For all *monkey F* data sets, average MDI was significantly greater than average MSI (*P* < 0.001, 1-tailed *t*-test). For *monkey T* data sets, average MDI was always greater than average MSI, but because of lower unit counts, these differences were statistically significant for only three of five data sets (*P* < 0.05, 1-tailed *t*-test).

To aid in interpreting this uneven breakdown of direction- vs. speed-encoding units, we simulated spike counts from several relevant encoding models. We fit each model to the neural activity and movement kinematics (nondiscretized) from the representative data set (F081909) and then simulated spike counts using the same real kinematics. As expected, when we simulated from the direction-only encoding model, no units were identified with MSI significantly greater than MDI (Fig. 5*A*), and similarly, when we simulated from the speed-only encoding model, no units were identified with MDI significantly greater than MSI (Fig. 5*B*). The information pattern from the velocity-encoding population (Fig. 5*C*) resembled that from the direction-only population, but since speed is a fundamental component of velocity, MSI values were slightly larger in the velocity-encoding population. Even so, none of these simulated units had an MSI value that was significantly larger than its corresponding MDI value.

The information signature from the simulated speed- and direction-encoding population (Fig. 5*D*) best resembled that of the real data (Fig. 4*A*), with a large number of units with significantly greater MDI than MSI in addition to a small number of units showing the opposite trend. This similar breakdown of direction- vs. speed-dominated units should be expected since the data were generated from a model fit to the real data. The direction-only, speed-only, and velocity encoding models result in theoretically prescribed distributions of MDI and MSI values. This distribution for the independent speed and direction model, however, can favor either direction or speed unit types, depending on the data. Also, note that the distributions of MDI and MSI values from these simulations are biased toward slightly smaller values than those from the real data in Fig. 4*A*. These differences speak to the fact that the real neural activity contains movement information not captured by the parametric tuning models used in these simulations [as observed, for example, by Churchland and Shenoy (2007)], yet this information is captured by the mutual information computations employed over the real neural data in this analysis.

Tuning indices (TI) from the linear regression analysis, defined as from fits to the tuning models of *Eqs. 2–4*, are shown in Fig. 6 as a function of lag between kinematics and neural activity. Direction TI curves closely matched the direction information curves of Fig. 2, and similarly, speed TI curves closely matched the speed information curves. Velocity TI curves typically had maxima that exceeded both the corresponding speed and direction TI maxima (although a few exceptions can be found). These velocity TI curves were more closely matched to the TI and information curves of direction than of speed but did not appear to be a simple function of one or the other. We note that TI values for direction and velocity are not directly comparable to those for speed because of differences in numbers of parameters between models and because TI values were computed over the same data used to fit the models. Rather, this analysis was motivated *1*) to help carry over intuition from previous studies framed from a regression perspective and *2*) to demonstrate that the appearance of velocity tuning does not necessarily predict the quantity of speed- or direction-related information that may be extracted from a population.

### Population Activity Enables Better Predictions of Direction Than of Speed

Results from the information analysis suggest that, at least in single-unit activity, the encoding of movement speed is substantially weaker than that of movement direction. To determine whether this finding holds true when considering the joint population activity, we applied a series of PNB classifiers toward predicting kinematics from population responses. Classifiers were trained to predict discretized kinematics based on *1*) a single 30-ms spike count aligned in time with movement kinematics (instantaneous) or *2*) the entire causal history of nonoverlapping 30-ms spike counts beginning 300 ms before the movement kinematics (history). Direction predictions were significantly more accurate than speed predictions under both the instantaneous and history conditions and across all data sets. For the representative data set detailed in previous sections (F081909), instantaneous direction accuracy was 26.0%, whereas speed accuracy was only 9.8%, as shown in Fig. 7, *A* and *G*. Incorporating spike count history into predictions for this data set increased direction accuracy to 38.7%, whereas speed accuracy only increased to 13.1%. Although these prediction accuracies may seem low on an absolute scale, they are actually relatively high given that predictors had to choose from 26 possible labels for both speed and direction, and as such, chance prediction accuracy was only 3.8%. These trends were consistent across all data sets, with direction accuracy 2.34 ± 0.25 times higher than speed accuracy for instantaneous predictions and 2.80 ± 0.23 times higher for predictions based on spike count history. Prediction accuracies for all datasets are tabulated in Table 3.

To summarize the full distribution of predictions, we also computed mutual information between predicted and actual discretized kinematics. As shown in Fig. 7, *B* and *H*, direction predictions carried more information than did speed predictions. As a performance metric, information complements prediction accuracy in that information provides a summary of the structure of predictions, including both the predictions that matched the actual kinematics as well as those that did not. Specifically, if two sets of predictions have the same fraction of correct predictions, information will be higher for the set whose incorrect predictions are less uniformly distributed across labels. In this regard, information can appropriately account for near misses, for example, if a classifier frequently predicts an incorrect label corresponding to a speed that is only slightly higher than the speed corresponding to the true label.

Confusion matrices are shown in Fig. 7, *C* and *D*, and Fig. 7, *I* and *J*. Whereas adjacent speed labels correspond to adjacent speed ranges, no such natural ordering exists for three-dimensional directions. To compensate, we provide column-reordered confusion matrices in Fig. 7, *E* and *F*, and Fig. 7, *K* and *L*, whereby the rows in each column have been sorted by angle (direction) or absolute difference (speed) between kinematics corresponding to actual and predicted labels. The confusion matrices show that incorrect prediction labels typically clustered around the correct label for both speed and direction. These distributions were tighter for direction than for speed, resulting in greater information values for direction than for speed.

We chose to include causal neural activity as the input to classifiers to mimic the real-time prediction problem that a BMI is required to solve. However, the information analysis revealed maximal speed information at acausal lags for many units (Fig. 3). When repeating the decoding analysis using both the causal and acausal histories of neural activity, we found prediction accuracies were largely unchanged compared with the corresponding accuracies using only casual neural activity.

To ensure that our discretization procedure is not responsible for these discrepancies between direction and speed prediction accuracies, we analyzed speed prediction accuracies as a function of speed bin widths used for discretization (data not shown). For predictions based on instantaneous spike counts, there was never a significant effect of bin width on prediction accuracy. For predictions based on spike count history, bin width had a small but significant effect on two of nine analyzed data sets; however, these two experiments had low unit counts relative to the other experiments.

To determine the effect of population size on prediction accuracy, we performed a unit-dropping analysis. As expected, predictions become more accurate with increased population size for both movement speed and movement direction (Fig. 8*A*). We note that the confidence intervals in the latter portion of the neuron-dropping curves (i.e., for numbers of units approaching the actual recorded population size) will be biased to be smaller than they actually are due to the similarity across draws from the actual population. However, even when accounting for this, extrapolation of these accuracy curves beyond the numbers of units we recorded suggests that, had we recorded a larger sample of neurons, speed prediction accuracy would likely remain substantially lower than that of direction predictions (data not shown).

We also computed PNB prediction accuracy as a function of the number of contributing units for simulated population recordings (Fig. 8, *B–E*). Consistent with the single-unit information analyses, the independent speed- and direction-encoding model resulted in a population with PNB prediction accuracies best matched to those of the real data. The key corresponding features are *1*) the ratio of direction-to-speed prediction accuracy across population size and *2*) the nearly saturated speed prediction accuracy when all units are incorporated into predictions. However, this simulated population gives systematically lower prediction accuracies for both speed and direction relative to those given by the real recorded data. This discrepancy again speaks to the fact that the real recorded neural data contain movement information not captured by the parametric tuning models used for these simulations and that PNB classifiers are capable of extracting this information from neural activity.

As described, the information and prediction analyses have treated speed and direction separately, characterizing each kinematic variable's relationship to neural activity independent of the other variable. We also performed the prediction analysis using a joint discretization scheme whereby classifiers were trained to jointly predict speed and direction from each of 26^{2} possible pairs of discretized speeds and directions. Because these classifiers required learning many more parameters, some data sets were not large enough to support the analysis. For the data sets that were large enough, marginal prediction accuracies for direction were again significantly greater than those for speed, although overfitting of the increased numbers of parameters produced absolute accuracies that were slightly lower than those reported in our main results (data not shown). To mitigate overfitting, we restricted the analysis to in-plane trials and decreased the number of speed and direction bins to eight each. This joint decoding analysis was well defined for all data sets and produced similar results to an analogous analysis where speed and direction were each decoded independently (data not shown).

### Difficulties Extracting Speed May Explain Deficiencies in BMI Control

Previous BMI studies have noted subjects' difficulties in controlling BMI cursor speeds, especially with respect to stopping and holding a cursor at a desired target location (Carmena et al. 2003; Ganguly and Carmena 2009; Gilja et al. 2012; Hochberg et al. 2006; Kim et al. 2008). To align with these studies, we implemented a BMI cursor control task using a velocity-only Kalman filter (VKF), a state-of-the-art neural decoder for BMI applications. Example cursor trajectories under VKF are shown in Fig. 9, as well as in Supplemental Video 1. (Supplemental material for this article is available online at the *Journal of Physiology* website.) During experimental sessions when all trials had minimal target hold requirements (50 ms), cursor trajectories were swift and straight into the targets (Fig. 9*A*). To quantify the subject's ability to crisply stop at targets, we introduced substantial hold requirements such that trial success required the cursor to overlap the target for a randomized hold time (0–600 ms), and a trial was failed if the cursor exited this acceptance region during the hold time. During these experimental sessions with substantial target hold requirements, the subject demonstrated poor control of movement speed under VKF, as evidenced by frequent trial failures due to overshooting through the target (Supplemental Video 1). Successful trials often involved meandering trajectories (Fig. 9*B*) such that cursor speed was relatively low upon initial acquisition of the target. When target hold requirements were minimal, this meandering behavior was not observed, and cursor speeds were substantially higher upon target acquisition (Fig. 10), suggesting that the subject adopted a VKF-specific strategy whereby stable stops were replaced by slow movements over the target. Consistent with this strategy, the subject's performance under VKF decreased substantially as target hold time requirements increased (see Fig. 12*A*).

We have shown that *1*) in single-trial arm reaches, speed information is relatively deficient in motor cortical activity compared with the abundant levels of direction information, and *2*) motor cortical activity alone cannot support precise control of BMI cursor stability under VKF. Several studies have highlighted the differences between closed-loop BMI control and offline analyses of arm control (Chase et al. 2009; Lebedev et al. 2005). To establish a link between the apparent deficiency of speed information offline and deficient speed control online, we simulated closed-loop BMI control as decoded by VKF. The underlying simulated neural population had log-linear tuning to the independent speed and direction model (*Eq. 15*), which was fit to real recorded neural activity from the arm reaching experiments. Thus the simulated population encoded substantially more information about desired direction than about desired speed. As shown in Fig. 11, simulated control was consistent with the real VKF BMI behavioral data in that success rates decreased rapidly as target hold requirements increased. This result suggests that the amount of speed information found during arm movements is consistent with the subject's deficient ability to hold at targets, and furthermore, that the amount of speed information required for stable stopping under VKF BMI control exceeds the amounts that we found in recorded motor cortical activity.

### SDKF Restores Stopping Ability During Closed-Loop BMI Control

From the information and prediction analyses, we found our ability to extract speed-related information from motor cortical activity to be relatively poor, despite using methods that require only mild assumptions about how movement speed might be encoded in neural activity. The resulting implication for BMI is that, even if it were possible to perfectly extract the limited speed information available in recorded neural activity, this may never enable reliable closed-loop control of BMI cursor speed because the encoding for movement speed is simply not that strong. We designed the SDKF to overcome this limitation in the decodability of movement speed. SDKF leverages the subject's reliable control of movement direction to improve control of movement speed by implementing a tradeoff between speed and angular velocity in the decoded velocity signal.

Example cursor trajectories under SDKF are shown in Fig. 9*C*, as well as in Supplemental Video 2. SDKF significantly enhanced the subject's ability to stop and hold for the duration of target hold requirements, as shown in Fig. 12*A*. For target hold times between 300 and 600 ms, success rates were 1.7 times higher under SDKF control than under VKF control. Cursor movements under SDKF were typically straight toward the target, rather than meandering (Fig. 9), indicating that the subject could instruct a crisp stop upon acquiring the target. We applied a constant speed gain to decoded SDKF velocities such that movement times were matched between VKF and SDKF trials (Fig. 12*B*). In this setting, SDKF achieved improved cursor stability at targets with movement times that were not significantly different from VKF movement times for hold requirements longer than 100 ms (for hold requirements under 100 ms, movements did not need to slow down substantially at targets to achieve task success).

SDKF improves BMI performance by leveraging natural features of goal-directed movements, as well as by potentially encouraging strategies specific to the feedback equations defining SDKF. Goal-directed movements, in both natural reaching and BMI settings, tend to begin with a high-speed, straight ballistic phase and tend to end with low-speed corrective movements relying heavily on sensory feedback. SDKF detects these corrective movements in the form of increased absolute angular velocity and correspondingly slows or stops the BMI cursor. This feature of SDKF is akin to a “hockey stop,” whereby a fast-moving hockey player makes a quick rotation of the skates to bring about a crisp stop on the ice. For SDKF, the result is a BMI cursor that automatically slows down near the target in response to corrective movements, in contrast to the overshooting behavior typically produced by standard BMI decoders. Figure 13 shows decoded angular velocity for SDKF and VKF as a function of distance to target. For both decoders, absolute angular velocities are low during the ballistic phase when the cursor is far from the target, and as the cursor approaches the target, angular velocities increase. Interestingly, SDKF trajectories showed larger angular velocities near the target compared with VKF, suggesting that the subject may have adopted a strategy of exaggerating turns near the target because doing so would be decoded by SDKF as a crisp “hockey stop.”

SDKF was designed to improve online BMI control of cursor speed. However, SDKF's performance benefits for online control, especially those that may be attributable to SDKF-specific control strategies, need not result in improved performance when reconstructing arm movements offline. We applied both SDKF and VKF toward offline reconstruction of arm velocity and found root-mean-squared reconstruction error to be nearly twice as large for SDKF relative to that for VKF (Fig. 14). This result highlights the fact that decoding algorithms with superior performance in online cursor control do not necessarily achieve superior performance in offline reconstruction, and as such, optimizing neural decoders offline cannot always be expected to yield the best decoders for online BMI applications (Chase et al. 2009).

## DISCUSSION

We asked whether the moment-by-moment details of movement speed could be extracted from motor cortical activity, and although we did find significant speed-related activity, our ability to extract movement speed was substantially worse than our ability to extract movement direction. In single-unit information analyses, we found roughly threefold higher frequencies of direction-dominated units compared with speed-dominated units. In population decoding analyses, we were able to predict movement direction with more than double the accuracy of corresponding speed predictions. These results are problematic for BMI systems, which depend on the ability to extract kinematic variables, including movement speed, from population activity on a moment-by-moment basis. To address this problem, we designed a BMI decoding algorithm, SDKF, which increased the ability to stop and hold the BMI cursor at instructed targets by 70.8%.

### Information and Prediction Analyses: Sensitivity to Modeling Choices

The information and prediction analyses required specific data processing to ensure a fair comparison of the extractability of speed vs. direction from neural activity. Because speed and direction are continuous-valued quantities expressed in different units and with differing numbers of degrees of freedom, we discretized speed and direction such that their discretized distributions had matched marginal statistics. We chose a discretization that resulted in roughly equal numbers of data points assigned to each of 26 speed labels and each of 26 direction labels.

With 26 discretization labels, chance prediction accuracy is 1/26 = 3.8% for both speed and direction, i.e., the accuracy of the best predictor that does not have access to the underlying neural activity. We chose to use 26 labels because there were 26 targets in the reaching task. Movement speeds lie in a continuum with no natural set of boundaries, and thus the number of speed labels must be arbitrary. When the prediction analysis was repeated with different numbers of discretization labels, {10, 15, 20, 40}, the results were consistent with those that we report when using 26 discretization labels. The numbers of units with significantly greater MDI than MSI remained consistent across all data sets, as did the numbers of units showing the opposite trend. Similarly, direction prediction accuracy was significantly higher than speed prediction accuracy for all data sets.

A logical alternative binning scheme is to discretize movement speeds using bins of constant width, which would result in substantially different numbers of data points across speed labels. With this alternative discretization scheme, we found that a single low-speed label can account for up to 20% of data points. Here, a chance predictor that always predicts the most frequent speed label will have an accuracy of 20% without considering the neural activity, thus complicating the ability to compare speed predictability with direction predictability. Even so, we found that speed prediction accuracy under this alternative discretization scheme was only a few percent better than chance.

Finally, for simplicity, the information and prediction analyses ignore temporal autocorrelation in the kinematics data. For the center-out reaching task, movement direction tends to be very similar across timesteps in a single trial. Movement speed tends to have more temporal variability, because speeds tend to have bell-shaped profiles throughout each trial. Incorporating such structure ought to improve direction prediction accuracy by more than it would improve speed prediction accuracy, thus only conservatively biasing our findings.

### Movement Representations in Motor Cortex

Previous studies have identified speed-related information in the activity of single motor cortical neurons (Churchland et al. 2006; Ifft et al. 2011; Moran and Schwartz 1999b; Schwartz 1992, 1994), as well as in motor cortical signals recorded from intracortical local field potentials (Heldman et al. 2006), electrocorticography (Anderson et al. 2012), magnetoencephalography (Jerbi et al. 2007), positron emission tomography (Turner et al. 2003), and functional magnetic resonance imaging (Rao et al. 1996). We found speed-related information as well, but when quantified relative to direction-related information, the extracted speed signals appear surprisingly weak. We can think of three potential interpretations of these results. The first is that instantaneous speed is robustly represented in motor cortex, but our analysis techniques were incompatible with the details of the neural encoding. The second is that the motor cortical representation of instantaneous speed, although weak relative to direction, is strong enough to enable robust control of movement speed. The third is that instantaneous speed is not robustly represented in motor cortex, and other factors that we did not consider combine to implement movement speed. We discuss each of these possibilities in turn below.

If motor cortex does encode the fine-timescale details of movement speed, why did our analyses not reveal a robust speed signal? One possibility is that the subset of units carrying reliable speed information might change depending on other kinematic parameters such that the population of neurons actually controlling movement speed dynamically changes. For example, it has been demonstrated that there is an interaction between direction and speed such that speed modulations are only apparent in the firing rates of a neuron during movements in the preferred direction of the neuron (Moran and Schwartz 1999b; Schwartz 1992). Decoding techniques that explicitly account for such dependencies might enable a more robust extraction of movement speed. Another possibility is that speed may be encoded more broadly across motor cortical populations that are substantially larger than those recorded in this study. Our unit-dropping analysis in Fig. 8*A* shows a shallow slope in the speed accuracy curve as more neurons were added to the decoder, suggesting that even if we had recorded from greater numbers of neurons, speed predictions might still be substantially less accurate than direction predictions. However, array recordings are typically biased toward monitoring populations of neurons located on cortical gyri, and it may be that speed information can be more readily extracted from neurons in the less accessible banks of cortical sulci. A third possibility is that speed is carried through a different neural code than we assumed. We applied methods to identify kinematic information in spike counts of motor cortical neurons. Although suggestive, our results leave open the possibility that movement speed is encoded in patterns of spike timing rather than spike counts. Spike timing information has been found in other systems, such as the rat whisker system (Panzeri et al. 2001), the mouse visual system (Jacobs et al. 2009), and the primate auditory system (Chase and Young 2008). However, reports of spike timing codes in the motor system have been limited (although, see Hatsopoulos et al. 1998). Finally, we note that the center-out task does not impose explicit requirements on movement speed throughout a reach, and it might be possible to design a reaching task that better modulates population activity with respect to movement speed.

Could the amount of speed information we found be enough to support precise control of arm speed? It is difficult to know how much information is necessary to enable the degree of speed control that our subjects exhibited during arm reaching. Furthermore, control over the moment-by-moment details of speed was not an explicit requirement in the arm reaching task. Rather, our subjects simply had to move to the target within a specified amount of time and maintain stability in the target for a specified hold period. Further experiments will be required to determine how much information is necessary to enable to precise control of speed in arm movements. In simulation, however, we found that populations of neurons carrying the amounts of speed information that we measured in real neurons demonstrated deficiencies in cursor stability at targets similar to those demonstrated in real BMI control by a monkey (Fig. 11). Although the simulation analysis was framed in the context of BMI, we believe it also has implications for natural movement control. Since the simulated neural encoding was fit to real neural activity underlying arm movements, these results suggest that arm movements rely on more than a readout of movement kinematics from motor cortical activity.

A third interpretation is that motor cortex is not the sole arbiter of movement speed, but may coordinate with other brain areas that contribute toward movement speed control (Tan et al. 2009). The role of M1 in driving movements has been extensively debated (for reviews, see Schwartz 2007; Scott 2003). Whereas much evidence has been presented in support of M1 encoding instantaneous movement details (Georgopoulos et al. 1982; Morrow and Miller 2003; Scott and Kalaska 1997), several studies have suggested a more dynamics-based encoding (Aflalo and Graziano 2007; Churchland et al. 2012). For example, it may be that M1 specifies a desired peak speed for a particular movement and that the motor periphery is responsible for generating the fine-timescale dynamics of movement speed. More generally, the speed signal in cortical activity may not be isomorphic with arm speed. The signal from motor cortex must be understood as only one factor combined with additional neural processing in the many other neural structures with speed-dependent activity, transformed by musculoskeletal action to produce arm movement. A better understanding of motor cortical operations and their contribution to arm movement will make it possible to develop more accurate extraction algorithms for decoding the details of this behavior.

An important distinction in the present study is that we sought to quantify the information in single-trial, simultaneously recorded population neural activity about kinematics during a single 30-ms timestep, because this is the relevant timescale for online BMI decoding. Moran and Schwartz (1999b) identified a robust speed representation in trial-averaged data generated from sequentially recorded units whose responses were combined as a population. We repeated that analysis with the data from the current study and also found a robust speed representation (data not shown). This correspondence suggests that although speed is encoded across sequentially recorded populations, it is difficult to extract from populations of simultaneously recorded units in the real-time setting of BMI, possibly because of correlated noise in single-trials that can be suppressed when averaging across trials.

### Implications for BMI Control

Similar to natural reaching movements, BMI cursor movements require precise speed control. Typical approaches to decoding BMI movements assume a relatively simple encoding of speed (e.g., linear through a velocity tuning model as in VKF). We found that an independent speed and direction model matched the neural data better than direction-only, speed-only, and velocity-only models. However, the information and prediction analyses in this study imply that adjusting these modeling assumptions, e.g., by using a nonlinear encoding model, may still result in limited BMI performance with respect to decoded speed because the moment-by-moment details of neural firing do not appear to carry requisite levels of speed information.

To overcome the apparent limitation in available speed information, we designed SDKF using a novel approach toward achieving high-fidelity control of speed in a BMI. SDKF incorporates a well-controlled neural signal, that of movement angular velocity, to improve on the low-fidelity speed signal present in neural activity. When a tradeoff is incorporated between movement speed and angular velocity, speed accuracy is improved without neural activity being required to supply an accurate speed signal. This SDKF design feature was informed in part by natural arm movements. First, as previously mentioned, natural arm movement speeds are influenced by the dynamics of the muscles and spinal cord, which may possibly alleviate the need for M1 to specify the moment-by-moment details of movement speed. In this sense, SDKF has a biomimetic interpretation in that the history dependent trajectory model (*Eq. 19*) imposes speed dynamics that are not directly specified by the neural activity. Second, natural arm movements have been shown to demonstrate a tradeoff between speed and curvature, often referred to as the two-thirds power law (Lacquaniti et al. 1983). Neural correlates of this relationship have been reported in previous studies of motor cortical activity underlying arm movements (Moran and Schwartz 1999a; Schwartz 1994). We used angular velocity (the temporal derivative of direction) as a proxy for curvature (the spatial derivative of direction) to simplify the BMI implementation.

Although the implemented tradeoff between speed and angular velocity enables SDKF to apply speed information not directly specified by the neural activity, the tradeoff alone is not sufficient to supply all of the requisite speed information. For example, changes in movement speed may be desired when changes in direction are not, especially for straight movements typical in a center-out task. For this reason, we incorporate the tradeoff between speed and angular velocity as an additional mechanism to complement the speed control implicit in SDKF's velocity encoding model. Future work will be needed to determine how well the SDKF trajectory model generalizes to tasks requiring curved movements (e.g., pursuit, circle drawing). It remains to be seen whether speed dampening is assistive or restrictive in these tasks.

A trivial means of improving cursor stability at targets is to simply slow down the decoded cursor movement. A slower cursor provides the subject with more time to instruct corrective movements to avoid inadvertent overshoot upon target acquisition. However, this approach increases movement times and decreases the overall throughput of the BMI. We applied a constant speed gain to SDKF such that movement times were matched between SDKF and VKF. With higher success rates for the same movement time, SDKF achieves a substantially higher throughput than does VKF. Rather than choose speed gains to match movement times, we could have matched mean movement speeds (VKF speeds were slightly faster on average than SDKF speeds). In this case, we would expect SDKF movement times to be shorter than those for VKF, but potentially at the expense of success rates for longer hold times. BMI decoders are inherently subject to this speed-accuracy tradeoff (Gowda et al. 2012), and in future experiments it may be worthwhile to specifically probe this tradeoff by evaluating decoders across a range of speed gains.

Previous BMI studies have proposed alternative approaches to solving the “cursor-stopping problem.” One approach is to directly decode a discrete target variable, such as movement end-point (Shanechi et al. 2013; Srinivasan et al. 2006; Yu et al. 2007), which could then be used to either constrain a subsequent trajectory estimate or to generate automatic control signals for acquiring the target. Additional approaches are to decode using a nonlinear neural tuning model that directly incorporates intended movement speed (Li et al. 2009) or to decode a binary stop signal (Hwang and Andersen 2009; Kim et al. 2011; Nuyujukian et al. 2012; Velliste et al. 2010). Taking an alternative approach, Gilja et al. (2012) recently demonstrated improved cursor stopping by applying assumptions based on feedback control. Finally, providing other modalities of sensory feedback (in addition to visual feedback) might help the subject better control BMI movement speed (O'Doherty et al. 2011). SDKF offers a complementary solution that enables the user to continuously guide and stop the cursor, while relying relatively little on the capacity for neural activity to directly specify movement speed.

Gilja et al. (2012) provide a comparison across studies in terms of Fitt's throughput, as computed by (20) (21) In the present study, the distance between workspace center and target center was 85 mm, and cursor and target radii were each 7 mm. Because target acquisition was defined by cursor-target overlap, the effective window size was 14 mm. For SDKF trials with required holds between 0 and 600 ms (expected hold time was 300 ms) and without cursor recentering, mean acquire time was 1.24 s, resulting in a throughput of 2.28 bits/s. The algorithm of Gilja et al. (2012) achieved throughputs of 1.48 and 1.81 bits/s in a task requiring 500-ms holds and without cursor recentering. We provide these numbers to approximately align between studies; however, differences between subjects and differences in trial structure may make exact comparisons impossible. Another important distinction is that in our task, trials were failed if target acquisition was lost at any time before the hold requirement was satisfied, whereas Gilja et al. (2012) continued trials until the hold requirement was satisfied while keeping track of the “dial-in” time between the initial target acquisition and completion of the target hold. For our study, success rate summarizes the subject's ability to stop and hold, whereas “dial-in” time is the analogous metric used by Gilja et al. (2012). Neither of these statistics factor into Fitt's throughput, and as such, Fitt's throughput cannot be used in this form to summarize stopping ability. A more direct comparison between SDKF and the aforementioned approaches might prove insightful in future work, and we believe that a combination approach leveraging the innovations presented across these studies is likely to yield the best results.

The general design principles underlying SDKF demonstrate the potential for performance gains when highly controllable neural modulations, previously used to drive one subset of control dimensions in the BMI task space (e.g., movement direction), are tapped to improve control across other dimensions of the BMI task space (e.g., movement speed). Importantly, these performance gains may be enhanced through subjects' adoption of cognitive control strategies that are effective in an online setting when paired with a decoder designed to be compatible with these strategies (e.g., instruct a sharp turn when a crisp stop is desired).

## GRANTS

This work was supported by National Science Foundation Integrative Graduate Education and Research Trainee Fellowship, National Institutes of Health Collaborative Research in Computational Neuroscience Grant R01 EB005847, and Pennsylvania Department of Health Research Formula Grant SAP#4100057653 under the Commonwealth Universal Research Enhancement program.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

M.D.G. and S.M.C. conception and design of research; M.D.G. performed experiments; M.D.G. analyzed data; M.D.G., B.M.Y., A.B.S., and S.M.C. interpreted results of experiments; M.D.G. prepared figures; M.D.G. drafted manuscript; B.M.Y., A.B.S., and S.M.C. edited and revised manuscript; M.D.G., B.M.Y., A.B.S., and S.M.C. approved final version of manuscript.

## ACKNOWLEDGMENTS

We thank Angus McMorland, George Fraser, and Jeong-Woo Sohn for arm reaching data.

## Appendix

#### Optimization for Direction Discretization

To discretize movement direction, we labeled each timestep in a data set according to the direction centroid, of 26 candidate centroids, whose angle with the measured movement direction was smallest. We designed the following optimization procedure to choose the set of direction centroids such that the discretization would result in approximately the same number of data points for each discretization label.

Given a set of direction centroids, we discretized the actual movement directions and computed the entropy of the resulting discretization (A1)
where **X** is the initialized set of direction labels and *p*(*x*) is the fraction of data points whose minimum angle direction label was *x*. Next, we selected one direction centroid and applied a small random rotation. We recomputed the entropy after discretizing the data using this perturbed direction label and the remaining 25 unperturbed labels. If this entropy was greater than the entropy prior to the random rotation, the rotated direction label was accepted. Otherwise, the rotation was rejected, and the set of direction centroids reverted back to the set prior to this random rotation. This process was repeated over 5 million iterations.

We initialized the procedure with the set of 26 target directions from the arm reaching task. In each iteration, we randomly selected the direction centroid to perturb to be either the centroid that labeled the most number of data points or the centroid that labeled the least number of data points. If the entropy had not increased after 1,000 consecutive iterations (i.e., no random rotations were accepted), the centroid to perturb was selected at random from the full set of 26 centroids. Random rotations were applied by *1*) defining a 3D unit vector in the direction of the unperturbed centroid, *2*) perturbing this unit vector by adding to each coordinate an independent draw from a Gaussian distribution with standard deviation 1 × 10^{−4}, and *3*) projecting the perturbed unit vector back onto a unit sphere.

By construction, this procedure is guaranteed to produce a sequence of nondecreasing entropies and is thus guaranteed to converge to either a local or global maximum. The theoretical maximum entropy is log_{2}(26) ≈ 4.7 when the *p*(*x*) are equal for all direction labels. If the number of data points *n* in a data set does not divide evenly into 26, the theoretical maximum is achieved when the *p*(*x*) differ by at most 1/*n*. In practice, optimized discretizations resulted in entropies that were within 4 × 10^{−5} of the theoretical maximum entropies.

#### VKF Algorithm

The Kalman filter predicts the subject's intended movement velocity given all recorded neural activity up to the current timestep. The Kalman filter prediction is a distribution over intended velocities, which takes the form of a multivariate normal distribution, i.e., *P*(**v**_{t}|**y**_{1}, . . . , **y**_{t}) = *N*(**v̂**_{t}, Σ_{t}). At each timestep *t*, the Kalman filter algorithm estimates the expected velocity, **v̂**_{t}, and a corresponding uncertainty, Σ_{t}, given all neural activity up to the current timestep.

Kalman filter predictions are computed recursively such that the prediction at a given timestep is computed using the prediction from the previous timestep. First, the trajectory model from *Eq. 17* is used to project previous predictions through a one-step update (*Eqs. A4* and *A5*). To determine the relative contributions of this trajectory-only update and the current neural activity, the Kalman gain is computed (*Eq. A6*) by integrating the uncertainties due to the trajectory and encoding models. This gain term is then used to incorporate the current neural activity into the current prediction (*Eq. A7*). Finally, the uncertainty of this prediction is computed based on the uncertainty from the one-step update, but reduced to reflect the information gained from the current neural activity (*Eq. A8*). VKF provides a stable decoding system (residual velocities will decay to 0 if neural inputs remain constant at baseline) when the maximal eigenvalue of **A**_{t} − **K**_{t}**CA**_{t}, from *Eq. A7*, is less than 1. The mathematical description of the complete VFK algorithm is as follows:

initialize: (A2)
for *t* ∈ {1, 2, . . . }
(A3)
(A4) (A5)
(A6)
(A7)
(A8)

#### SDKF Algorithm

To incorporate a tradeoff between speed and angular velocity, SDKF dampens the decoded speed when the recently decoded cursor trajectory exhibits a large absolute angular velocity. Since angular velocity is ill-defined at near-zero speed, speed dampening is reduced when cursor speeds are low, enabling the cursor to accelerate from stops. These design features are implemented through an extension to VKF, whereby SDKF incorporates the adaptive trajectory model described in *Eq. 19*. At each timestep, SDKF computes the direction of the most recently decoded velocity (*Eq. A9*) and the change in direction since the previous timestep (*Eq. A10*), wrapped to remain between −180 and 180 deg. Next, the mean angular velocity is defined as the average change in direction over the most recent three timesteps, which corresponds to 100 ms (*Eq. A11*). Angular velocity-based speed dampening (*Eq. A12*) and speed-based speed dampening (*Eq. A13*) are combined using *Eq. A14*, where λ_{t} ∈ [0, 1] is a time-varying speed-dampening factor. SDKF exactly reproduces VKF decoding when λ_{t} = 1. For 0 ≤ λ_{t} < 1, the one-step update in *Eqs. A4* and *A5* effectively shrinks the velocity prior toward 0, dampening the decoded speed relative to the corresponding VKF decode. The mathematical description of the SDKF algorithm is the same as that for the VKF algorithm but replacing *Eq. A3* with the following sequence:
(A9) (A10)
(A11)
(A12) (A13) (A14) (A15)

We manually selected α = ⅓ and β = 8 to achieve the desired speed dampening during preliminary experiments and fixed the parameters during all analyzed experiments. Speed dampening is shown as a function of speed and angular velocity in Fig. A1. As defined, SDKF's speed dampening can decrease decoded speeds but cannot increase them. To match movement times between VKF and SDKF, we multiplied SDKF-decoded velocities by a constant speed gain factor of 3.

- Copyright © 2014 the American Physiological Society