## Abstract

A pursuit-tracking task (PTT) and multielectrode recordings were used to investigate the spatiotemporal encoding of hand position and velocity in primate primary motor cortex (MI). Continuous tracking of a randomly moving visual stimulus provided a broad sample of velocity and position space, reduced statistical dependencies between kinematic variables, and minimized the nonstationarities that are found in typical “step-tracking” tasks. These statistical features permitted the application of signal-processing and information-theoretic tools for the analysis of neural encoding. The multielectrode method allowed for the comparison of tuning functions among simultaneously recorded cells. During tracking, MI neurons showed heterogeneity of position and velocity coding, with markedly different temporal dynamics for each. Velocity-tuned neurons were approximately sinusoidally tuned for direction, with linear speed scaling; other cells showed sinusoidal tuning for position, with linear scaling by distance. Velocity encoding led behavior by about 100 ms for most cells, whereas position tuning was more broadly distributed, with leads and lags suggestive of both feedforward and feedback coding. Individual cells encoded velocity and position weakly, with comparable amounts of information about each. Linear regression methods confirmed that random, 2-D hand trajectories can be reconstructed from the firing of small ensembles of randomly selected neurons (3-19 cells) within the MI arm area. These findings demonstrate that MI carries information about evolving hand trajectory during visually guided pursuit tracking, including information about arm position both during and after its specification. However, the reconstruction methods used here capture only the low-frequency components of movement during the PTT. Hand motion signals appear to be represented as a distributed code in which diverse information about position and velocity is available within small regions of MI.

## INTRODUCTION

Neural activity in primary motor cortex (MI) is correlated with aspects of arm motion such as hand position (Georgopoulos et al. 1984; Kettner et al. 1988), speed (Ashe and Georgopoulos 1994; Moran and Schwartz 1999a), direction of motion (Ashe and Georgopoulos 1994; Fu et al. 1995; Georgopoulos et al. 1982), and force (Sergio and Kalaska 1998; Taira et al. 1996). Most MI neurons appear to combine information about multiple movement features (Ashe and Georgopoulos 1994; Moran and Schwartz 1999a) that may be specified separately in time (Fu et al. 1995). The temporal aspects of the encoding process are important both for understanding the neuronal processing of dynamic signals (Buracas et al. 1998; Mainen and Sejnowski 1995; Rieke et al. 1997) and for the problem of decoding information from populations of neurons (Humphrey et al. 1970; Warland et al. 1997), yet previous work has not differentiated temporal patterns imposed by task demands from the underlying temporal dynamics of encoding.

Investigation of the spatiotemporal encoding of motor variables presents several challenges. Tasks used to study movement have most often involved point-to-point movements to a limited number of well-rehearsed targets. Step-tracking tasks, as typically implemented, allow only limited control over kinematic variables because hand motion is a function of the subject's strategy rather than of the experimental design. For example, in a typical point-to-point movement task, any hand velocity can be used to reach a target as long as target acquisition falls within a maximum allotted time. In addition, typical step-tracking tasks limit the size of the parameter space sampled for each variable; studies of target location encoding are typically limited to a small subset of possible locations (8, in the widely used “center-out” task; Ashe and Georgopoulos 1994; Georgopoulos et al. 1982; Kalaska et al. 1989; Moran and Schwartz 1999a). Furthermore, because hand position and velocity are strongly interdependent in these tasks, it is difficult to determine their relative contributions to MI firing. For example, in the standard radial task, any given peripheral position is associated with just one single direction of motion, and with a highly stereotyped set of velocity profiles.

Another problem— especially significant for studies of temporal dynamics—in tasks typically used to study motor coding is that neural and behavioral variables (such as firing rate and hand speed) are statistically nonstationary. Distributions of these measures vary systematically as a function of trial time, so that, for example, peak firing occurs within a narrow interval after a cue to move. Nonstationarities in the underlying data distributions greatly complicate the analysis of temporal encoding processes because lag-dependent interactions (those related to coding delays) are confounded with trial-time-dependent modulations in activity.

In earlier studies in motor cortex, behavioral variables were treated as static, scalar quantities such as average hand direction or speed, and the concomitant time-varying neural activity was summarized as a single number—the mean firing rate. The data were averaged over many trials and/or fit to highly parametric tuning models (e.g., cosine functions), thereby collapsing what may be more information-rich tuning functions (Sanger 1996). These multiple averages eliminate most of the dynamic, trial-specific information needed to characterize spatiotemporal encoding properties. In contrast, more recent studies have explicitly examined the temporal aspects of kinematic coding in MI using center-out-type tasks (Ashe and Georgopoulos 1994; Fu et al. 1995; Moran and Schwartz 1999a; Sergio and Kalaska 1998) or curved drawing tasks (Moran and Schwartz 1999b; Schwartz and Moran 1999), and treating the kinematics and neural activity as time-varying data. These studies avoid the issues of collapsing data across time but still suffer from the inherent statistical constraints described above, that is, interactions between variables of interest, and the confounding of *time*-dependent with *lag*-dependent properties, where lag is the delay between spiking and its manifestation as behavior. Temporal dynamics, as they have been studied in the context of these tasks, could be an indication of the temporal evolution of task demands and are not necessarily an indication of the underlying dynamics of encoding.

Finally, the serial recording techniques employed in previous work preclude the direct comparison of spatiotemporal encoding properties between neurons because units are recorded under behavioral and state conditions that vary from trial to trial (and therefore from cell to cell). Furthermore, serial recordings of neural data necessitate assumptions of statistical independence between neurons (because the dependencies cannot be observed without simultaneous recording), and these assumptions have been shown to be inaccurate in general (Maynard et al. 1999; Oram et al. 2001).

The present study characterized spatiotemporal encoding of hand motion using a random, continuous pursuit-tracking task (PTT) designed to facilitate evaluation of the spatial and temporal characteristics of MI neurons, while minimizing dependencies and nonstationarities. Using continuous tracking of a randomly moving stimulus, position and velocity encoding is characterized within a systems analysis framework. In this context, hand trajectory is viewed as a random “stimulus” to the system and neural activity is the “response.” Each stimulus is drawn from an experimenter-determined distribution that broadly and continuously covers velocity and position space, and is stationary with respect to trial time. This design effectively controls hand motion at all times and reduces statistical dependencies among variables across the experiment. These statistical properties of the PTT permit the rigorous application of information-theoretic and signal-processing methods to the analysis of position and velocity coding. The relationship between kinematics and firing rate can be characterized in a nonparametric (model-free) manner, without assumptions about the underlying tuning properties of the sampled neurons. The multielectrode recording approach taken here allows quantitative comparisons of encoding between cells, because multiple neurons are recorded under completely identical conditions. Finally, the systems analysis approach further permits a direct quantification of hand trajectory information using signal reconstruction methods that can demonstrate planned motions from population activity. In this paper we describe the spatiotemporal tuning functions of MI neurons for velocity and position during pursuit tracking and we compare the information coded within single cells and across the population. We also demonstrate that MI neurons contain sufficient position and velocity information to reconstruct novel hand trajectories based on information available from the firing of a small sample of MI neurons.

Part of this work appeared in abstract form (Society for Neuroscience Meeting 1999; abstract 665.9; Society for Neuroscience Meeting 2001; abstract 940.1).

## METHODS

### Behavioral task

Three monkeys (one *Macaca fascicularis* and 2 *M. mulatta*) were operantly conditioned to track a smoothly and randomly moving visual target. The monkey viewed a computer monitor and gripped a two-link, low-friction manipulandum that constrained hand movement to a horizontal plane. Manipulandum position was sampled on a 30 × 30-cm digitizing tablet (Wacom Technology, Vancouver, WA) at 167 Hz, with an accuracy of 0.25 mm, and recorded to disk. Hand position was continuously reported on the monitor by a black, 0.2° visual angle circle (0.5 cm radius on the tablet) (Fig. 2*A*).

At the beginning of each trial, a red, 0.6° (1.5 cm tablet radius) tracking target appeared in a random position, drawn from a 2-D, zero-covariance Gaussian (up to the cutoff imposed by the edge of the screen) distribution with mean located at the workspace center. The monkey was required to align the feedback and target cursor within 1.5 s (4 s for monkey Ra); if the target was not acquired, the trial was aborted and the target reappeared at a new, independently, identically distributed (i.i.d.) position to begin the next trial. A 700-ms hold period followed target acquisition, after which the target began to move in a smooth, but random fashion. If the monkey continuously maintained the feedback cursor within the target for 8-10 s, a juice reward was delivered. Each target trajectory stimulus was a randomly generated i.i.d. signal that was presented only once: the target position (and thus to first-order, hand position) during the tracking period was generated by, in essence, running Gaussian white noise through a band-pass filter, with the horizontal and vertical components generated independently. More specifically, a spectrum was constructed, consisting of 2^{17} integer frequency components, such that the power was 1/f within the band-pass and 0 otherwise. Each frequency component was assigned a different, random phase. This spectrum was then inverse Fourier transformed producing the position signal in the time domain. This signal was then scaled appropriately for the workspace and resampled at 8 Hz. The power spectrum of the resulting signal for one experiment is shown in Fig. 3*C*. Note that this is not identical to the original spectrum because of the finite length of the signal. Spectra for other experiments were qualitatively similar, because, by construction, they were identical with the exception of the bandwidth, which was left as a free parameter and varied between experiments (see Table 1). Importantly, the 1/f characteristic of the band-pass filter for the position signal means that the velocity signal is approximately white within the band-pass, and thus has minimal autocorrelation width for that given band-pass.

For comparison, 2 of the monkeys were also trained to perform a standard “center-out” task (Fig. 1; see Georgopoulos et al. 1982; Maynard et al. 1999 for details). In these experiments, radial and tracking trials were randomly interleaved. The monkeys used in this study had been trained on the center-out task before introduction to the continuous tracking task. All 3 animals were able to perform the tracking task within the first 2 days of training, with varying degrees of proficiency; performance (as measured by the length of time for which the monkey could consistently track a target of given mean speed; see Table 1) continued to improve throughout the training period. Data analyzed here were collected 8-11 mo after introduction to the tracking task.

### Recordings

Details of the basic recording hardware and protocols are available elsewhere (Donoghue et al. 1998; Maynard et al. 1999). After task training, a Bionic Technologies LLC (BTL, Salt Lake City, UT) 100-electrode silicon array was implanted in the arm representation of MI. The array was placed on the precentral gyrus medial to a line extending from the genu of the arcuate sulcus posteriorly to the central sulcus, and parallel to the sagittal fissure, a region previously localized as the MI arm representation (Georgopoulos et al. 1982). The neurons showed modulations around movement commonly observed in radial task experiments, thus confirming, physiologically, placement of the array in the MI arm area. The BTL arrays consisted of 100 platinizedtip silicon probes (about 200-1,000 kΩ at 1 kHz; Nordhausen et al. 1996), arranged in a square grid (400 μm on center). The electrodes were 1 mm in length, corresponding in MI to recordings near the layer III/V boundary. In the 3 monkeys there were 74 (Monkey Ra), 47 (Er), and 24 (Co) possible active recording electrodes, a number limited by the connectors used. All procedures were in accordance with Brown University Institutional Animal Care and Use Committee-approved protocols and the Guide for the Care and Use of Laboratory Animals (National Institutes of Health publication no. 85-23, revised 1985).

Signals were amplified and sampled at 30 kHz/channel using a commercial recording system (Bionic Technologies, Salt Lake City, UT). All waveforms that crossed a manually set threshold were digitized and stored (from 0.33 ms before to 1.17 ms after threshold crossing); spike sorting to isolate single units was performed off-line. Single units with signal-to-noise (SNR) ratios >2.5 were stored as spike times referenced to the stimulus signal for further analysis. Analysis of spiking was confined to data recorded from 1 s after tracking began to 1 s before the end of trial, to eliminate nonstationarities associated with trial beginning and end.

### Analysis

spatiotemporal tuning. We summarized the spatiotemporal tuning of the recorded cells as follows. We computed functions *N*(, τ) and *N*(, τ) to describe the firing rate as a function of position () and velocity (), respectively, at a series of time leads and lags (τ). These functions are defined as the conditional mean firing rate of a cell at time *t*, given that a particular kinematic value ( or ) occurred at time *t* + τ. That is and where *E*(· | ·) denotes conditional expectation, *R* is the spike rate, and τ defines the delay between the spike count bin and the kinematic bin [i.e., *N*(, -100 ms) gives the expected firing rate 100 ms after the particular hand position was observed].

To compute these tuning functions, data were taken at all times {*t*_{i}} when the hand was moving with a particular velocity (or was located at a particular position) (ρ ± dρ, θ ± dθ) cm/s, for some (ρ, θ) in polar coordinates. The bin widths 2dρ and 2dθ were chosen to be just large enough to ensure adequately sampled data in all bins; we typically took >50 samples per bin. For example, we set bin widths in one experiment to 0.4 radians × 0.7 cm/s (velocity), and 0.4 radians × 0.5 cm (position). We then calculated the mean firing rates at {*t*_{i} - τ} for the lags τ shown. We represent this lag variable by the symbol τ throughout, reserving *t* for the time since the beginning of the behavioral trial.

We used polar instead of rectangular coordinates for the discretization for 3 reasons. First, polar coordinates respect the radial symmetry of the (properly scaled) observed Gaussian joint distributions of hand position and velocity (Fig. 4): all bins at a given radius ρ are roughly equiprobable, whereas the corresponding statement is false for any fixed value of horizontal or vertical position or velocity. Second, the size of the bins in polar coordinates (approximately ρ dρ dθ) grows with ρ, partially correcting for the falloff of the probability distribution of these behavioral variables at the extremes of their ranges. Finally, in polar coordinates firing rates are represented as a function of direction, a convention that facilitates comparisons with prior studies. The origin for these curves was taken to be the mean of the distributions of the behavioral variable; for the velocity tuning functions, the origin was at (0, 0) cm/s, whereas for position the origin was at the center of the tablet. We fit planes and other parametric families to the tuning curves by a standard least-mean-squares optimization procedure (Nelder-Mead simplex search). In addition, we used a Monte Carlo procedure to obtain conservative significance levels for the presence of good fits, under the null hypothesis that the spike trains were homogeneous Poisson processes (i.e., that the apparent fluctuations in firing rate observed in Fig. 2 were random, had a trivial probabilistic structure, and were independent of the behavior of the hand). We simulated spike trains (homogeneous Poisson processes with rates matched to the observed individual neural firing rates), estimated *N*(, τ) and *N*(, τ) using real kinematic data for each instantiation of these simulated spike trains, and computed the mean-square deviation for each resulting fit. The significant fit level was taken as the point at which the cumulative empirical probability distribution of the random goodness-of-fit value reached 0.99.

information-theoretic analysis. Mutual information is a nonparametric measure of dependency that is capable of detecting dependencies that correlational measures ignore. The mutual information between the random signals *N* and *S* is defined as (Cover and Thomas 1991) 1 where *p*(·) and *p*(· | ·) denote marginal and conditional probabilities, respectively, and ∫_{X} is the integral over some space *X*. Information is difficult to compute in general because full knowledge of the joint distribution *p*(*N*; *S*) (where *N* and *S* are functions of time) is needed. This presents a possibly infinite-dimensional learning problem; in the present experiment one would be required to know the probability of a given spike train given any time-varying position signal. Consequently, we do not attempt to estimate the information rate between spike trains (denoted *N*, for neuron) and the behavioral signal (*S*); rather, we address the simpler problem of computing the information between the observed neuronal firing rate and the behavioral signal (hand velocity or position, here) at discrete (*single*) time lags τ; that is 2 *N*(0) here denotes the activity of the given neuron in the current time bin, and *S*(τ) denotes the state of the behavioral signal (e.g., the position of the hand) at time lag τ after the present time; computing *Eq. 2* requires only an integral in 2-D space (one dimension each for horizontal and vertical), instead of the high-dimensional integral required to compute the full information (*Eq. 1*) between spike trains and the time-varying position signal.

To simplify *Eq. 2* even further, we modeled the conditional distributions of the behavioral signal given an observed spike count per bin, *p*[*S*(τ) | *N*(0) = *i*], *i* ∈ 0, 1, 2,..., as Gaussian, with mean μ_{τ,i} and covariance matrix σ_{τ,i}. This simplification makes the computation of *Eq. 2* tractable, given the size of the available data set. Thus for *Eq. 2*, we calculate 3 numerically, where *G*(μ, σ) is the (2-D) Gaussian density with mean μ_{τ,i} and covariance σ_{τ,i}. The Gaussian model was motivated by empirical observations and gave a sufficient fit to the data for many observed cells and spike count bins, according to a 2-D Kolmogorov-Smirnov test (bivariate Kolmogorov-Smirnov-type test; Press et al. 1992; *P* < 0.05). In the cases in which the Gaussian fit was inadequate, we applied a nonparametric binning approach (computing the integral in *Eq. 2* as a finite sum) instead; the Gaussian and binned-information estimates were highly correlated (correlation coefficient > 0.95) across all cells and all time bins, indicating that the Gaussian method provides an adequate information estimator for this set of data.

A Monte Carlo procedure identical to the one described in the previous section was used to estimate significance levels for the observed information values. This procedure produced information values <10^{-4} bits. A different procedure, in which we shuffled the neural data with respect to the behavioral data, so that neural data from one trial was associated, in a random manner, with the behavioral data from a different trial, led to similar results. The significance bound was therefore defined as *I*[*N*(0); *S*(τ)] > 10^{-4} (see Fig. 12).

signal reconstruction. The ability to reconstruct aspects of hand motion from multiple, simultaneously recorded spike trains was used as a test of availability of position or velocity information in the recorded population. We used a multiple linear regression approach (Neter et al. 1985): our estimate *R* (for reconstruction) of the position at the current time *t* is given by a linear combination 4 where *i* indexes time; *j* is the cell number; *N*(*i, j*) denotes the activity of cell *j* at time *i*; *a*_{i,j} represents the corresponding “weight”; *C*, the number of cells; *T*_{pre} and *T*_{post}, the time before and after the current time *t* used to estimate the current position, respectively; and d*t*, the width of the time bins used. The filter coefficients *a*_{i,j} were computed as in Warland et al. (1997). Two filters were generated, one each for the horizontal and the vertical positions.

The analytical solution to the optimal linear estimation problem in the time domain involves the inversion of a correlation matrix (**N ^{T}N**) that can be fairly large [matrix size =

*D*

^{2}, where

*D*= 1 +

*C*(

*T*

_{pre}+

*T*

_{post})/d

*t*]; we used standard singular value decomposition (Press et al. 1992) techniques to check the numerical stability of this matrix inversion. The data showed no evidence of overfitting such as a decrease in performance as

*D*became large. None of the results shown was smoothed, nor were any relevant parameters subjectively selected (e.g., to select the “best” neurons for analysis). Cross-validation methods were used to estimate the expected error of our reconstructions: we fit the regression model to a “training” set consisting of all but 10 trials of the data set, then computed the mean-square error of the regression on this “test” set, the 10 held-out trials. This process was iterated multiple times as successive, disjoint blocks of 10 trials were used to test the regression; we report the regression coefficient computed by this procedure, where this coefficient is defined as usual as

*r*

^{2}= 1 - {

*E*[(

*R*-

*S*)

^{2}]/

*E*(

*S*

^{2})}, where

*R*is the reconstructed hand position and

*S*is the true hand position.

A frequency domain regression analysis (Haag and Borst 1998; Rieke et al. 1997) was used to estimate a lower bound on the frequency content of the information contained in the MI population (Fig. 15). Neural and position signals were Fourier transformed, and the neural Fourier coefficients at a given frequency ω, *N̂*(ω), were regressed onto the coefficients of position, *Ŝ*(ω), to obtain the reconstruction of *S* at ω, *R̂*(ω). Goodness of reconstruction was plotted as the SNRs obtained at each frequency 5 where *E*(·) denotes the sample mean (with the number of samples here equal to the number of trials), and ^{*} denotes a complex conjugate. The bound on information rate was calculated, as usual, from Shannon's formula (Cover and Thomas 1991; Rieke et al. 1997).

Finally, the reconstruction error was examined as a function of *1*) the total length of time (*T*_{pre}) spike trains were observed and *2*) the number of neurons included in the analysis. We examined the dependency of the estimation error on *T*_{pre} by recalculating *r*^{2} for several different values of *T*_{pre} (Fig. 16*A*). The analysis of *r*^{2} versus the number of cells (Fig. 16*B*) is slightly more complicated, given that the regression error is a function of not only how many cells one chooses to observe, but also which subset of cells is chosen. Therefore neurons from a simultaneously recorded data set were randomly selected and the range of *r*^{2} obtained for each such randomly selected subset was plotted. For reasons of computational efficiency, we did not use the cross-validation method to compute *r*^{2}, but rather used the equation , which gives the expected *r*^{2} given that the true covariance matrix of *S* is σ_{ss} and the cross-correlation between *N* and *S* is σ_{ns}; *N* here is a vector-valued signal, with each element corresponding to the firing rate of a single cell, and *E*(·) denotes expectation. In practice, σ_{ss} and σ_{ns} must be estimated from data, and because of sampling error, the *r*^{2} computed by cross-validation tends to be of lower magnitude than the *E*(*r*^{2}) calculated here; therefore we normalize the curves in Fig. 16 by the maximal observed *E*(*r*^{2}).

neural stationarity. We tested neural activity for trends in both the firing rate over the course of each experiment and the firing rate across trial time. The firing rate as a function of time (intratrial or across the experiment) was fit by a line and the slope was tested to see whether it was significantly different from zero. This was done through a bootstrap procedure. Tests were done separately for each cell. See the appendix for details.

Cells exhibiting significant trends in rate over experimental time were further tested for significant changes in their spatiotemporal tuning functions over experimental time. Those cells with significant rate changes and significant tuning changes were discarded. Cells exhibiting significant intratrial rate changes were not excluded (see results). Of an original 120 cells, we excluded 7 because of nonstationarities, leaving the 113 we use in all subsequent analyses.

## RESULTS

Eleven data sets from 3 monkeys were analyzed. These data consisted of 6-17 min of tracking behavior recorded simultaneously with neural data from 3-19 single units (median = 11 min and 11 cells; see Table 1). In total 113 (of an original 120; 7 were not analyzed because of nonstationarities; see *Neural stationarity* below) neurons were analyzed. We first describe behavior and neural activity during the PTT and compare them with data from the radial task. Next, we report spatiotemporaland temporal-tuning functions for individual MI neurons during the PTT, and finally, we discuss results of a linear reconstruction technique for extracting behavioral signals from these neurons.

### Pursuit-tracking task

The pursuit-tracking task (PTT) and typical point-to-point movement tasks vary considerably in the extent of parametric space explored, the dependencies among variables, and the stationarity of kinematic and neural signals. Figure 1 illustrates kinematic and neural activity data obtained from one monkey performing the center-out task, to provide explicit comparison with the PTT. The center-out task, by design, results in movements from a constant location to one of a fixed set (here, 8) of discrete locations. Although there is no specific trajectory requirement, the need to end at a specific location within task-time constraints generally results in roughly straight, stereotyped hand trajectories. Figure 1*A* shows hand paths for trials to each of the 8 directions. This task design results in strong dependencies between horizontal and vertical position (Fig. 1*A*) and horizontal position and velocity (Fig. 1*B*). Note, also, that many (x, y) pairs, even near the center of the workspace, are never sampled. Figure 1, *C* and *D* illustrate the nonstationarity of kinematic and neural variables in the center-out task: mean hand speed shows a sharp transient increase with movement onset, irrespective of target location (Fig. 1*C*), and mean firing rates show similar large *t*-dependent modulations (recall that *t* denotes time relative to the start of the trial).

By contrast, the PTT covers the kinematic space more fully and achieves considerably improved independence of kinematic variables and stationarity of kinematic and neural activity (Figs. 2, 3, 4). Figure 2*A* provides an example of PTT performance for a single trial. Tracking was smooth, with continuous modulation of hand speed and direction. Mean hand speed, which followed that of the visual target set in the experimental design, ranged from 2.5 to 4.7 cm/s across this set of experiments (Table 1). Tracking movements were largely determined by the visual stimulus, as demonstrated by the close temporal relationship of the hand and visual cue (Fig. 2*A**, inset*). The peak of this cross-covariance was consistently located within 50 ms of zero with a peak correlation coefficient that exceeded 0.97 in each data set, consistent with the conclusion that the animals tracked the stimulus. The short visuomotor “reaction time” indicates that the animal is at times actively predicting the smoothly evolving stimulus trajectory. The relatively high tracking accuracy over time can also be appreciated in the individual plots of x and y position versus time across a trial (Fig. 2, *B* and *C*). The overall smoothness of hand movement during tracking is evident in the autocovariogram (Fig. 3*A*), and in the power spectrum of hand position (Fig. 3*B*); most of the power in the hand position signal was below 1 Hz (Fig. 3*B*; the autocovariogram and power spectra in Fig. 3 were computed from data from a single experiment, but these functions were qualitatively similar in each other data set). For comparison the power spectrum of the horizontal position of the stimulus signal is shown (Fig. 3*C*); again, most of the power is below 1 Hz.

Figure 4 presents the statistical properties of the PTT for comparison with those of the center-out task (cf. Fig. 1). The joint distributions of 2-D hand position and 2-D velocity in the PTT were well approximated by Gaussian distributions with zero covariance (modified Kolmogorov-Smirnov test; *P* < 0.05), as expected given the task design. No significant correlation was observed between any of the pairs of velocity and position variables (Pearson test; *P* < 0.05). Thus the PTT samples the kinematic space more densely than does the center-out task. In addition, kinematic variables such as hand speed and position are effectively stationary across the task. Mean hand speed does not vary as a function of trial time (*P* < 0.05; compare Figs. 4*C* and 1*C*) and average firing rate does not depend on the time relative to the start of tracking for the cells shown (test on correlation with linear trend over the first or last 2.5 s of the trial; *P* < 0.05; compare Figs. 1*D* and 4*D*). Figure 4*D* is shown for illustrative purposes because, for some cells in our database, the average firing rate was not constant over trial time (e.g., some cells displayed anticipatory “ramp-up” activity near the end of successful trials). Any intratrial rate nonstationarities during the PTT cannot be explained as a function of the variables of interest (i.e., the kinematics) because these variables are stationary. The comparison between Figs. 1*D* and 4*D* is meant to show that the center-out task *induces* rate nonstationarities, whereas the PTT does not.

### Neural activity during tracking

Figure 2*D* shows a representative example of the spiking patterns of 21 cells recorded simultaneously during a single pursuit-tracking trial. Qualitatively, randomly selected MI neurons typically showed varying modulation patterns in the PTT; these same neurons showed marked mean rate modulations in step-tracking tasks (compare Figs. 1*D* and 4*D*). Mean firing rates during the PTT ranged over 1.5 log units (about 2-40 Hz; Fig. 5) and were not correlated with overall mean hand speed (Spearman rank-order correlation coefficient; *P* < 0.05). The relationship between the spike count mean and variance (per 50-ms bin) is largely linear with unity slope, except at the highest mean firing rates, where the Fano factor (the ratio of the variance to the mean) falls slightly below the unity level.

### Neural stationarity

Our results depend on the stationarity of the underlying data. By construction, the stimulus (i.e., the motion of the tracking target) is stationary; thus the animals' hand motions are approximately stationary. This does not, however, guarantee the stationarity of the neural activity associated with these motions. In averaging over the entire experimental time period to derive our tuning measures we are implicitly assuming that tuning is constant on this time scale. Because the subjects are well trained on the task before recording, and the task requirements are held constant across the experiment, there is good reason to think that this is true—no learning is likely to be occurring. However, changes in the animal's overall behavioral state (e.g., motivation) might cause average spike rates to drift up or down over a recording session. To test for this we looked for linear trends in the average spike rate for each cell across experimental time.

Cells with a linear trend whose slope was not significantly different from zero, or with less than a 20% change in rate, were deemed stationary on the experimental time scale and included in the other analyses. Cells with a significant nonzero slope and a change in rate of >20% over the experiment were further tested for trends in their spatiotemporal tuning functions (see following text). Of an original 120 cells we found 44 (37%) with significant (by bootstrap shuffling of time bins, *P* < 0.05) rate trends over the experiment. Of these, 7 (5%) were found to have tuning functions that differed significantly (see methods) over experimental time. These cells were excluded from further analysis, leaving the 113 reported here.

We also tested for stationarity of rate as measured across trial time. For each experiment we aligned trials on the beginning of the tracking phase and averaged the neural activity for each cell across trials to get a mean firing rate for each time bin. We tested for linear trends in the average rate over the course of trial time. We found 27 (23%) of 120 cells with significant (by bootstrap shuffling of time bins, *P* < 0.05) rate trends of >20% over trial time. No cells were excluded based on these intratrial rate trends. Because the kinematics are stationary over trial time these intratrial trends in rate are unlikely to be linked to the tuning that we report. The fact that intratrial trends, when they were present, were different for different, simultaneously recorded cells (e.g., some cells had a positive rate trend, whereas others showed a negative rate trend) also supports the idea that it is not the kinematics that are inducing these changes. It is likely that other, uncontrolled and unobserved variables (e.g., reward expectation) are inducing these rate trends. For these reasons, we argue that these effects may be interesting in their own regard, but do not detrimentally influence the results reported here.

### Spatiotemporal tuning

The spatiotemporal tuning properties of MI neurons were defined from the time (lag)-varying tuning of the cell with respect to velocity or position signals (see methods). Conceptually, using each spike time as a reference point for sampling of the kinematic variable, one can determine the spatial information provided by firing about that variable at any time in the future or the past, relative to that spike time. Spatiotemporal tuning functions for 113 single MI neurons were generated for velocity and position [denoted *N*(, τ) and *N*(, τ), respectively]. These functions summarize a neuron's instantaneous firing rate dependency on hand velocity or position, , at different delays τ, where τ is the time difference between a particular hand motion variable sample and the observed firing rate sample. A lead (τ > 0) is the amount of time the neuron was firing in advance of that kinematic measurement, whereas a negative τ represents a lag.

Figure 6 illustrates the spatial features of velocity [*N*(, τ)] and position [*N*(, τ)] tuning, at a single value of τ, for 2 different neurons. Tuning functions are plotted first in rectangular coordinates (Fig. 6, *A1, B1*) and then transformed into polar coordinates (Fig. 6, *A2, B2*; see methods). Polar coordinates are adopted for the remaining figures to simplify comparisons between position and velocity tuning. The origin for these tuning surfaces is taken as (0, 0) for velocity, and the center of the tablet workspace for position (in each case, the origin was the mean and mode of the observed kinematic distribution (see Fig. 4).

In polar coordinates the velocity tuning function plots firing rate against speed (ρ) and direction (θ); θ = 0 corresponds to movement to the right. The cell shown in Fig. 6*A* is approximately sinusoidally (i.e., cosine-) tuned for direction [i.e., the function *N*_{v}(ρ, θ, τ) can be fit by a cosine for any speed ρ]. The phase of this cosine is constant as a function of ρ, so that the direction tuning curve is approximately cosine as well (here *R* is some sufficiently large constant). Finally, the amplitude of this tuning curve scales approximately linearly with speed; the cell is in a sense more strongly tuned for direction at higher tangential velocities. A first-order model of this tuning function can be given by 6 where *a*_{0}, *a*_{1} > 0 are the baseline firing rate and constant “gain” parameters, respectively, and θ_{PD} is the cell's “preferred direction.” Because *Eq. 6* defines a plane in velocity space, we will refer to this model as the “planar model,” with *a*_{1} termed the “planar slope” parameter and θ_{PD} the “major axis.” This model has been shown to apply to MI firing during reaching (center-out) movements as well (Moran and Schwartz 1999a). For our data, the planar model for velocity gave a significant fit for 99% of the neurons in our sample (see methods). The data for Fig. 6*A* were recorded during an experiment in which pursuit-tracking and center-out trials were interleaved; by plotting the center-out target location tuning curve (Fig. 6*A3*) next to the PTT velocity tuning function (Fig. 6, *A1, A2*), we see that, for this neuron—although not necessarily for all neurons—the 2 concepts of tuning effectively coincide.

Neurons in MI were also tuned for hand position (Fig. 6*B*) during the PTT. For the position tuning functions in polar coordinates, the firing rate is plotted against distance from the origin (ρ) and direction (θ), where θ = 0 corresponds to rightward locations. Sinusoidal tuning in θ, similar to that observed in Fig. 6*A* for velocity, is evident. The firing rate increases linearly with ρ but maintains constant phase; that is, tuning functions for position are significantly fit by planes as well (98% of neurons). A planar model significantly fit MI tuning functions for both velocity and position for 90% of the cells in our database. In comparison, Kettner et al. (1988) found that 64% of neurons they recorded in the motor cortex arm area showed a linear relationship between firing rate and hand position, although, in their case, the hand was held static at each position. To examine whether tuning peaked at a particular value (e.g., akin to tuning of hippocampal place cells), we tested the fit of 2-D Gaussian functions for these tuning curves. The Gaussians provided a better fit to the position tuning functions for only 5 (4%) of the cells, and a better fit to velocity tuning for only 2 (2%) of the cells, despite the fact that the Gaussian function had 4 extra free parameters. Moreover, in each of these 7 cases, the width parameter in the Gaussian function was quite large, indicating the shallowness of the observed “peaks.” Thus the simple planar model in *Eq. 6* appears to be a reasonable first-order description of the 2-D tuning of MI cells for both position and velocity. The distribution of *R*^{2} values for fits to *Eq. 6* are shown in Fig. 9, *D* and *E*. In the following, the fit parameters of the planar model are used to summarize the tuning properties of the observed MI population.

Spatial tuning functions shown in Fig. 6 are representative of a single delay (τ), which fails to show the temporal dynamics of this tuning. Consequently, tuning was examined over multiple lags and leads τ. Figures 7 and 8 each show an example of spatiotemporal tuning functions for velocity *N*(, τ) and position *N*(, τ) for a single cell. These figures illustrate the heterogeneity of the temporal dynamics of MI tuning for these variables. Figure 7 depicts the most common MI tuning type. First, the cell is strongly velocity-tuned, especially at nonnegative delays (τ ≥ 0). Second, velocity tuning peaks at approximately τ = 100 ms, a lead consistent with the hypothesis that these cells signal upcoming observed hand velocity. Tuning begins to emerge several hundred milliseconds before this time and fades several hundred milliseconds afterward. Throughout this time the overall tuning structure remains essentially phase (θ) invariant. The temporal structure of this velocity tuning function *N*(, τ) is, for many cells, largely explained by a modification of *Eq. 6*, expressed as 7 where *a*_{1}(τ) is a smooth function of τ, with a maximum at 100 ms, such that *a*_{1}(τ) ≈ 0 for τ > 1 s. *Equation 7* is a useful heuristic for understanding how tuning evolves for most cells, in that it implies a fixed orientation (PD) over all τ. In no case do we see a smooth shift in PD over τ. That is, over τ, the gain (i.e., *a*_{1}) may go from positive to zero to negative—thus effectively abruptly flipping the PD by 180°— but the θ_{PD} term does not vary as a function of τ.

Position tuning showed a spatiotemporal structure that appeared to be directly related to velocity tuning for some neurons, but unrelated for others. The position tuning *N*(, τ) of the neuron in Fig. 7 can be explained in terms of the inherent dependencies between velocity and position (when considered as time-varying signals, not as static variables; cf. Fig. 4). To see why, assume that this cell's firing rate depends only on hand velocity. Nevertheless, hand velocity and position are necessarily correlated for most nonzero lags (although for PTT data this correlation is fairly weak for all lags, and zero for zero lag, as shown in Fig. 4). Whenever the hand is moving to the right at time *t* = 0, the mean position at time *t* = -ϵ will be to the left of the mean position at time *t* = +ϵ, for all sufficiently small positive times ϵ. Thus if we have a neuron signaling rightward *velocity* of the hand at τ ≈ 100 ms, as does the cell shown in Fig. 7, we should expect this neuron to signal the leftward *position* of the hand at negative time lags (τ = -1 s) and the rightward position at more positive lags (τ = +1 s), as observed here. Thus in this case, the position “tuning” of this cell can be explained parsimoniously in terms of its velocity tuning.

In contrast, Fig. 8 shows an example of a neuron whose position tuning cannot be readily explained from velocity tuning, suggesting that it specifically encodes position separately from velocity. In this example, position tuning is more pronounced and more temporally invariant than velocity; peak position tuning remains stable at θ ≅ π/4, whereas the velocity tuning peak changes from θ ≅ π/4 to θ ≅ -2π/3 between τ = -1 and τ = 0.88 s. Note that this change in phase is not a continuous shift, with peaks at intermediate angles, but a bimodal function in which, at intermediate values of τ, the tuning diminishes and then reappears. As described above, and consistent with *Eq. 7*, phase shifts of a more continuous (i.e., rotational) nature were not observed in this population. Having recorded this cell during an experiment in which pursuit-tracking and center-out trials were interleaved, we can observe that the center-out target location tuning (Fig. 8, *inset*) matches closely that predicted by integrating the spatiotemporal tuning function for *position*, but not velocity, over τ.

Figure 9 summarizes the spatial aspects of these velocity- and position-tuning functions. The distribution of the optimal planar angle (*a*_{1} in *Eq. 7*) and major axis (θ_{PD}) is shown for both position (Fig. 9*A*) and velocity (Fig. 9*B*). The distributions of θ_{PD} were indistinguishable from uniform on [0, 2π] for both variables (Kolmogorov-Smirnov test); that is, even within the small patches of MI sampled by the electrode array, a broad representation of hand position and velocity is present. The position and velocity major axes are weakly statistically dependent: when the differences modulo π between the major axes (Fig. 9*C*) are plotted, the position and velocity major axes for a neuron tend to be close [Kolmogorov-Smirnov deviation from uniformity (i.e., independent velocity and position θ_{PD}), *P* < 0.0001], as shown by the peak at 0. Position and velocity appear, for about half our recorded population, to be encoded essentially independently (Δ*PD* > π/8). For the other half (corresponding to the peak at zero in Fig. 9*C*) position and velocity tuning mirror each other, as in Fig. 7.

### Temporal dynamics of encoding

An information-theoretic analysis was used to provide a direct measure of position and velocity information available from the recorded neurons and to describe more quantitatively the temporal evolution of this encoding. The results in Figs. 6, 7, 8 demonstrate that by observing the position or velocity of the hand it is possible to derive information about the activity of a given MI neuron. The converse, by Bayes's rule, is also true: information about position or velocity can be decoded from MI firing rates. Figure 10 shows the conditional probability distributions, with corresponding Gaussian fits, of the horizontal hand velocity at *t* + τ, τ = 100 ms, given that this cell fired zero (Fig. 10*A*), one (*B*), 2 (*C*), or 3 (*D*) spikes within a 50-ms window around time *t*. The marked overlap in the set of curves demonstrates that the firing rate of MI neurons typically conveys highly ambiguous information with the small numbers of spikes observed in a narrow time window.

These conditional probability distributions can be used to quantify the temporal evolution of tuning in individual neurons. For this analysis the mutual information between the cell's firing rate and the kinematics of the hand is computed as a function of τ, *I*[*N*(0); *S*(τ)]. Here *N*(0) represents the cell's activity in a given short time interval (here, 5 ms; the interval is taken to be short to avoid redundancy effects induced by the fact that the hand position and velocity change relatively slowly) and *S*(τ) denotes the value of position or velocity some time τ before or after the current time, *t* = 0. This information statistic is an objective measure of how well these neurons are tuned for these behavioral variables; the more tuned a given cell is at a given value of τ, the more highly separated are the probability distributions corresponding to those shown in Fig. 10, and the higher the value of *I*[*N*(0); *S*(τ)]. Because this quantity is calculated directly from the underlying probability distributions it does not depend on any underlying assumptions about the linearity of the relationship between the neural firing rate and the behavioral variable, as do standard correlational statistics. The resulting curves, as functions of τ, discard all spatial tuning properties (e.g., preferred direction) and therefore show only temporal (τ-dependent) tuning features.

Figure 11 shows examples of information curves for hand velocity (Fig. 11, *A1*-*C1*) and position (Fig. 11, *A2*-*C2*), for 3 experiments. Individual curves within a panel (*A, B*, or *C*) and between panels (e.g., *A1* vs. *A2*, etc.) can be directly compared because the neurons shown were recorded simultaneously (and therefore the information curves were constructed using identical kinematic data). These temporal tuning curves were heterogeneous, especially in the position domain; some are unimodal, others multimodal, some peak at τ > 0 and others at τ < 0, all within the same set of simultaneously recorded data. The widths and shapes of the curves vary widely (note that the position curves change more slowly than do the velocity curves, partially because of the autocorrelation structure, as discussed above) and there does not appear to be any simple rule relating the curves for velocity and position. The width of the velocity information curves is uncorrelated with those of the corresponding position curve (Spearman's rank-order correlation coefficient, *P* < 0.05; test performed only on the 77 cells with significant velocity and position information content). This analysis also showed differences in the time at which peak information was available about position and velocity (Fig. 11, *D* and *E*). Temporal tuning peaks are always markedly more clustered for velocity than for position, with velocity curves consistently peaking near τ = 100 ms (i.e., firing leads behavior by 100 ms), and position peaks more temporally dispersed, suggesting that cells carry feedback as well as advance position information.

Figure 12 summarizes the information content of the observed MI cells and confirms quantitatively the considerable heterogeneity of these neurons. Information content ranges over 2 orders of magnitude. A weak but statistically significant correlation was present between velocity and position information content (Spearman's rank-order correlation coeffi-cient = 0.69; *P* < 0.05) (Fig. 12*A*). Moreover, information values did not cluster by mean hand speed during the experiment (Spearman's rank-order correlation coefficient; *P* < 0.05), indicating that position and velocity tuning is an intrinsic property of these cells, independent of the details of the particular sample of movements tested. On average, cells carry small amounts of information for position and velocity and only 10% more information for velocity than for position (median peak velocity information = 0.0011 bits per 5-ms bin; median peak position information = 0.0010 bits; 64% of cells carried more velocity than position information). Neurons with low average activity conveyed about as much information as those with high average firing rates, indicating that there is no dependency of information content on mean firing rate (Fig. 12*B*; Spearman's rank-order correlation coefficient = -0.06, n.s.). By extension, the information content of a given cell did not depend on firing rate variance (recall Fig. 5). It should be noted that, even though we find no dependency between mean hand speed and either *1*) position information, *2*) velocity information, or *3*) mean firing rate (Fig. 12), the range of mean speeds sampled is relatively small. It may be that dependencies would be observed if a broader range of mean speeds were tested.

Figure 13 graphically illustrates position information (Fig. 13*A*) or velocity information (Fig. 13*B*) versus planar gain (see Fig. 9), as derived from the fit to *Eq. 6*. In general, gain and information are correlated (correlation coefficient 0.65 for velocity, 0.73 for position). Gain increases as information increases, in keeping with the standard notion that a cell is more strongly tuned if its firing rate is more modulated by the variable of interest. Note, however, that there are cells that are *1*) well fit by the planar model (*R*^{2} > 0.5; “+” symbols in figure) and *2*) provide a relatively large amount of information, but *3*) have a relatively small gain. This is consistent with the idea that cell tuning is a function of not only the depth of modulation but also the variability in firing rate. This means that cells can convey large amounts of information about a variable even if they do not exhibit large, obvious rate modulations.

### Signal reconstruction

The preceding analyses demonstrate that individual MI neurons carry information about hand position and velocity. To determine what information is present in the MI population, we attempted to reconstruct, or decode, hand position from the activity of the population, using simultaneously recorded MI neurons. Hand position reconstruction at any given time *t* was estimated using a weighted linear combination of the neural activity from all observed cells, some time *T*_{pre} before and *T*_{post} after time *t* (Neter et al. 1985; Paninski et al. 1999; see methods). This linear correlation approach returned a moderately good reconstruction of the hand trajectory with no a priori assumptions (e.g., cosine tuning) on the tuning process other than linearity. Figure 14 shows 3 reconstructed signals for x and y position over time, as well as an example of reconstructed hand path (x vs. y), with the corresponding true signals for comparison. The quality of reconstruction is summarized by the usual correlation statistic *r*^{2} in Table 1. The performance of the linear estimator ranged from marginal up to about 50% of variance captured. The data in Table 1 also show that the observation of neural data after the kinematic event occurred (i.e., *T*_{post} > 0) robustly improves the reconstructions (Wilcoxon paired-sample rank test, *P* < 0.05), as expected given the results in Fig. 11*E*; this suggests that MI encodes something akin to a feedback copy of the ongoing hand motion in addition to the feedforward “drive” signal embodied in corticomotoneuronal cells.

The linear regression technique can be used to quantitatively evaluate which aspects of the tracking are contained in MI activity. For example, Fig. 15 indicates that reconstructions capture information only about the lower-frequency components of the hand trajectory. We quantified this observation by reconstructing the position signal directly in the frequency domain. Across the 11 experiments, the SNR of the reconstruction consistently dropped to the unity level by about 0.5 Hz (Fig. 15), indicating that the linear technique fails to extract information about hand position above this frequency from neural activity recorded during PTT. By contrast, the monkey tracks the stimulus at much higher frequencies, as shown in the average coherence plot between the hand and visual target position (dashed line, Fig. 15). This rapid falloff of SNR sharply limits the overall information rates for hand position in the observed MI populations; using Shannon's formula, we obtain rates of only about 2 bits/s, even in the most informative populations observed. This translates to a maximal rate of about 0.1 bit cell^{-1} s^{-1}.

The PTT makes it possible to examine the effect of the duration of spike observation (filter duration) and number of neurons on reconstruction quality. Reconstructions improved as spiking over longer times was considered (Fig. 16). To compare between experiments, the raw *r*^{2} values were normalized by the peak *r*^{2} observed during the given experiment, so that these curves range from 0 to 1. Figure 16*A* gives a sense of the typical trade-off between how quickly the reconstruction can be computed and reconstruction accuracy: the more time bins examined, the better the reconstructions, but at the cost of a greater delay in the reconstruction output. The slope of this *T*_{pre} versus *r*^{2} graph is quite sharply peaked near zero, indicating a kind of “diminishing returns” in *T*_{pre}: when *T*_{pre} is small, we have to observe relatively fewer neural data to achieve a given increase in *r*^{2} than when *T*_{pre} is already large. Increasing the number of neurons considered also improves reconstruction (Spearman rank-order correlation coefficient between number of cells observed and *r*^{2}; *P* < 0.05). However, the degree of improvement depends on which population of cells is observed (Fig. 16*B*). In particular, it is difficult to extrapolate from the curves shown here, to make any quantitative statements about the asymptotic behavior of the estimator as the number of cells observed becomes large (cf. Wessberg et al. 2000).

## DISCUSSION

The pursuit-tracking task, coupled with the multielectrode recording technique, enabled us to characterize 3 novel features of the relationship between motor cortical activity and hand movement. First, we were able to describe the temporal dynamics of position and velocity tuning as a function of lag and compare the resulting “spatiotemporal” tuning functions directly for simultaneously recorded cells. We found that these tuning functions wax and wane over time (lag); velocity information typically leads behavior and peaks within a narrow temporal window, whereas position tuning curves are much more heterogeneous. We did not find evidence for temporal segregation of tuning, either within individual cells or across the population tuning, because position and velocity could overlap in time. Furthermore, we confirmed, under novel (dynamically varying) behavioral conditions, that individual neurons encode both position and velocity information, and that both types of tuning can be locally summarized as cosine functions that are “gain modulated” by speed (Moran and Schwartz 1999a) or distance (Kettner et al. 1988). Second, information-theoretic analyses suggested a distributed representation of this kinematic information in MI; neurons carry approximately equal amounts of information for position and velocity, with absolute information values small for any single cell. Finally, linear regression and frequency domain techniques indicated that ensembles of primary motor cortical neurons best encode the low-frequency components of the hand position signal, but that a limited number of MI neurons is sufficient to predict random, smooth, 2-D hand trajectories with a moderate degree of precision.

### Pursuit-tracking task

In the PTT each hand path can be considered as a novel, time-varying “stimulus” for the motor system, with the neural activity representing the observed response. The theoretical strength of this analogy is debatable, but its empirical utility should be clear. For example, once we view movement in terms of a collection of time-varying signals (whether these signals are hand position and velocity, as analyzed here, or muscle tensions, joint angles, or any other behavioral signal), many of the points we emphasized above are immediate. First, it becomes clear why trial-averaging (averaging neural data over trials during which the temporal details of the relevant behavioral signals differ) might obscure essential details of the encoding process. Figure 4, *C* and *D* makes this point dramatically; here, trial-averaging destroys *all* information about the relationship between neuronal activity and behavior. Second, it is clear why control over the animal's movements is essential (for the same reasons that control over stimulus parameters is essential to a sensory physiologist); the PTT provides control over movements and attentive state because it demands continual visual monitoring of the stimulus to correctly guide the hand. Similarly, we see why it is important to study the response of the system to inputs from as large a portion of the relevant parameter space as possible, and why we need to be able to vary the multiple parameters of interest independently. For instance, the fact that directional tuning emerged from an analysis of the very large ensemble of random movements used here demonstrates that this property is not an epiphenomenon attributed to overtraining on a limited movement repertoire, or the statistical idiosyncrasies of radial-type tasks (cf. Fig. 1).

Most important, the stationarity of the PTT enabled us to treat movements as samples from a stochastic process. Each sample could be treated in a uniform (i.e., identically distributed) manner. There was no need to attempt to create a period of stationarity by dividing trials into behaviorally distinct epochs as done in step-tracking paradigms. This at once increases the effective size of our data set and allows the use of powerful statistical tools for systems analysis that depend on stationarity, such as frequency domain methods (Fig. 15) and all analyses of τ-dependent properties performed here (Figs. 7, 8, and 11). In contrast, during the radial task trial-time linked rate modulations occur on time scales of the order of hundreds of milliseconds (Fig. 1). These nonstationarities would contaminate the τ-dependent properties of the spatiotemporal tuning functions, which vary on a time scale of seconds in the case of position. These features make the PTT a potentially useful framework (albeit, of course, not the only such framework) to study other aspects of movement encoding.

### MI tuning functions

The tuning functions *N*(, τ) and *N*(, τ) examined in this study (Figs. 6, 7, 8) are analogous to the “spatiotemporal receptive fields” analyzed in various visual areas (DeAngelis et al. 1999), or “spectrotemporal” auditory fields (Kowalski et al. 1996), with one exception: the long correlation times of natural movement, compared with the signals used as stimuli in these sensory studies, cause our tuning functions to change more slowly in τ than do the functions derived in the sensory domain. The term “tuning function” is meant to be more neutral than “receptive field”; the results have not established what is actually *directly* encoded by the neuron, only what can be recovered from firing. Systems analysis approaches in sensory systems have revealed a similar diversity of tuning functions when neuronal firing is considered across the temporal and spatial domain (DeAngelis et al. 1999).

The recorded MI neurons typically showed spatial and temporal structure in their tuning for both velocity and position. Where comparisons could be made, our results on the spatial properties of MI tuning for position and velocity were generally consistent with previous reports. Both types of tuning showed a directional dependency fit by a cosine (Ashe and Georgopoulos 1994; Georgopoulos et al. 1982, 1984; Maynard et al. 1999; Todorov 2002). Direction tuning is stable across delay (τ) for a majority (60%) of the cells reported here. Direction tuning during the PTT showed a linear dependency on speed (or distance, ρ, for position tuning curves), as observed in center-out-like tasks (Ashe and Georgopoulos 1994; Hamada 1981; Hamada and Kubota 1979; Moran and Schwartz 1999a; Schwartz 1993). Our results show that speed and distance scale the directional tuning curve without affecting its shape; this relationship can be described locally by a simple planar model (see also *Eq. 6*; Georgopoulos et al. 1984; Schwartz 1993). Note that the near-planar form of these tuning functions implies that single MI neurons do not encode a particular location, which is quite different from the place fields of, e.g., hippocampal neurons (Brown et al. 1998). Our results complement previous work on spatial tuning in 2 main respects: first, by showing that planar fields persist in a dynamic behavior setting. Second, they provide greatly enhanced local detail about the tuning structure because of the higher-density sampling properties of the PTT. In addition, our work emphasizes the heterogeneity of the slope and orientation of the position and velocity planes within the relatively small region of cortex covered by the electrodes (Fig. 9), suggesting that these parameter spaces are fully represented within any given small patch in the MI arm area. This is consistent with the view that representations of arm control are very broadly distributed in MI (Sanes and Donoghue 1997).

The results of this study provide significant new information concerning the temporal properties of MI neurons, especially in the context of their spatial tuning (Figs. 7 and 8). Most of the cells showed spatiotemporal tuning for both position and velocity, with a continuum between strong velocity encoders (like the cell shown in Fig. 7) and strong position encoders (Fig. 8). The population of cells showed a broad range of properties. In some cases, position tuning could be understood as a feature of velocity tuning, whereas in other cases position seemed to be an independently coded variable. Although these features have not been demonstrated to be the actual variables encoded by MI neurons, at a minimum our results constrain the types of mechanistic models of this encoding process (Pugh et al. 2000; Todorov 2000).

Our evidence shows that neurons with heterogeneous velocity and position coding features are commingled even within a small volume of cortex, but does not support the hypothesis that position and velocity neurons form separate classes. Rather, encoding of these variables appears to be represented across a continuum in which these features are differentially weighted. Salinas and Abbott (2001) suggest that a mixture of cells with these sorts of encoding properties is well suited to form translation-invariant representations. In MI, this could mean that neural ensembles could represent particular kinds of movements irrespective of their particular location in space or, conversely, particular locations, regardless of motion. Such a network might also account for motor equivalence where the same action is produced, with structural similarity, from multiple effectors. Distributed, multiple representation with gain fields is also thought to be useful to provide signals that can be readily decoded by their target structures (Salinas and Abbott 2001).

Our data differ significantly from previous studies of the temporal properties of motor cortical cells. For example, Johnson et al. (1999b), using a hybrid pursuit/center-out paradigm, reported that speed tuning was specified before direction, with little overlap in speed and direction coding in time, suggesting that the 2 signals are not combined. Distinct temporal ordering was also found in a study using a center-out task with multiple radii (i.e., 8 directions with 6 distances each) (Fu et al. 1995). Cell discharge was first correlated with direction, then target position, and finally with distance.

We believe that the discrepancies between these findings and ours arise because the studies quoted address a fundamentally different question than that examined here. Specifically, these previous studies examined the way in which tuning tracks the evolution of task requirements. That is, the differences in temporal ordering of encoding can be attributed to the fact that nonstationary tasks, such as the standard center-out task, impose a particular temporal order due to the ordering of task requirements [i.e., the variables of interest are highly dependent on trial time (*t*); thus so is the neural activity with which they are correlated]. For example, as Johnson et al. (1999) suggest with regard to their results, preferred direction shifts during the delay period were related to alignment of visual and movement signals—an occurrence temporally linked to a behavioral epoch. In other words, it is likely that the demands of the task *evoke* an early correlation between hand speed and firing rate, followed later by a correlation with direction, given that a judgment of target speed would aid the animal in timing its interception before tracking. Similar arguments apply to the results presented in Fu et al. (1995). In contrast, as described below, the PTT does not impose any temporal ordering in the coding of kinematic parameters. This means that the dynamics we report describe the evolution of tuning as a function of the delay (τ) between spiking and behavior [rather than as a function of trial time (*t*)], a description that, to the best of our knowledge, has not been examined in detail before in the motor system. Note that lag-dependent tuning could also be examined in the context of nonstationary behavioral tasks (that is, tasks for which the distributions of the variables of interest depend explicitly on trial time *t*; recall, for example, the *t*-dependency of the mean hand speed in Fig. 1). However, to compute tuning dynamics in this case, one must, in general, examine a new tuning function not just for each τ, but also for each *t* (necessitating an average over a large number of trials, instead of the average over *t* we took advantage of here). Of course, this does not solve the undersampling problem facing tasks of radial type (recall Fig. 1); because a small set of similar hand paths are repeated in these tasks, any tuning functions computed from such data will implicitly be dependent on these particular trajectory histories, and might therefore function poorly as a general description of the cell's encoding properties.

### Temporal dynamics of position and velocity information

The purely temporal (that is, τ-dependent) properties of MI spatiotemporal tuning functions have not, to our knowledge, been previously studied in detail. We introduced a way to measure these temporal tuning properties, without any assumptions of linearity in the encoding process, by computing the “temporal tuning curves” for velocity *I*[*N*(0); (τ)], and position *I*[*N*(0); (τ)] (Fig. 11). Analysis of these objects quantified the heterogeneity visible in Figs. 7 and 8; the shapes of these temporal tuning curves varied considerably from cell to cell. This diversity was evident even when the curves were constructed using exactly the same behavioral data and thus cannot be explained in terms of kinematic or motivational differences between experiments; the ability to remove these confounds represents an important advantage of simultaneous multielectrode recording. Temporal heterogeneity was not predicted by previous work; nevertheless, this wide range of tuning properties is consistent with previous descriptions of the diversity of the correlation strength between neural activity and various behavioral parameters (Kakei et al. 1999; Porter and Lemon 1993).

Temporal tuning curves were heterogeneous not just in their shape but also in their overall amplitude (Fig. 12); information content for hand velocity and position varied over 2 orders of magnitude. Moreover, these information values did not depend on the mean firing rate of a given neuron (or on the variance of the firing rate; Fig. 5), or the dynamic range of the behavioral signal. In other words, this measure of information content appears to quantify an intrinsic property of MI cells, one that is relatively insensitive to these gross neural and behavioral parameters. However, information content for these variables might depend on other parameters (e.g., the posture of the animal or orientation of the arm) (Scott and Kalaska 1997), which were not systematically varied in the present study.

Finally, the information values computed for instantaneous position or velocity during tracking were perhaps surprisingly small, compared with those previously computed for static target location in the center-out task (Hatsopoulos et al. 1998; adjusted for differences in bin size). Firing rate appears to vary smoothly as a function of position and velocity, and the conditional distributions of the kinematic signal given spike count depended only weakly on the spike count (Fig. 10); in other words, MI cells are broadly, not sharply, tuned for these variables. This is in agreement with previous results, including those of Ashe and Georgopoulos, (1994), where ANOVA techniques performed on radial task data indicate that (static) target location accounted for much more of the variance in firing rate than did time-varying hand position, velocity, or acceleration. However, we found that the information content of MI cells for velocity was only 10% greater, on average, than the information content for position (Fig. 12), whereas Ashe and Georgopoulos (1994) found a much stronger preference for hand velocity than position. These discrepancies may be attributed either to the many statistical differences between the tracking and center-out paradigms or to differences between the (linear) ANOVA procedure employed in Ashe and Georgopoulos (1994) and the slightly more general information-theoretic analysis employed here. In addition, it is worth remembering that single-unit recording studies typically search for highly modulated cells, whereas we derived data from any well-isolated cell that could be retrieved from the electrode array, and thus these differences might result partially from selection bias. A more quantitative analysis of how these properties depend on cortical layer and area would be useful.

In sum, position and velocity are weakly encoded in the observed population of MI cells, when one compares the encoding of single-joint-related motor variables (Humphrey et al. 1970), or of higher-level variables such as target location (Ashe and Georgopoulos 1994; Hatsopoulos et al. 1998). Additionally, information content for velocity and position are weakly correlated (Fig. 12), which indicates that MI cells directly encode variables that are, in turn, indirectly linked to hand velocity and position. Correlation with other arm motion variables would help to determine which parameters are best represented by MI neurons.

### Signal reconstruction

We demonstrated that a linear algorithm, given a small, randomly chosen set of neurons and <20 min of training data, can reconstruct the random trajectory of a monkey's hand through 2-D space (Figs. 14, 15, 16, Table 1). (Neurons were “randomly chosen” in the sense that no preselection of well-tuned neurons was performed; all well-isolated units that happened to be within the recording range of our chronically implanted electrode array during a given experiment were analyzed.) Moreover, relatively small subpopulations of cells can capture significant fractions of the available information (Fig. 16*B*). The ability to reconstruct a trajectory using a simple algorithm from small sets of neurons suggests that it would be relatively straightforward to control devices in complex ways using limited neural samples from MI (Kennedy and Bakay 1998; Moxon et al. 1999; Wessberg et al. 2000) and our decoding approach (Serruya et al. 2002). Such neural prosthetics could be used to restore movement to paralyzed individuals as indicated by recent real-time control studies in monkeys (Serruya et al. 2002; Taylor et al. 2002).

Our approach was closest to that of Humphrey et al. (1970), Paninski et al. (1999), and Wessberg et al. (2000) (see Rieke et al. 1997 and references therein for similar approaches in various sensory domains, and Brown et al. 1998 for applications of various reconstruction methods in the study of hippocampal place cells). Humphrey et al. successfully estimated various single-joint-related parameters with time-domain linear regression techniques, using 3 to 8 simultaneously recorded MI units. More recently, Wessberg et al. (2000) obtained results similar to ours, using one task in which the hand was constrained to move along a 1-D track and another in which the monkey made stereotyped reaching movements. Our work is the first to show that nonstereotyped, random, multidimensional hand motion can be reconstructed with moderate accuracy. Although our estimation efforts, as well as those of Wessberg et al. and Taylor et al. (2002), successfully extracted a significant amount of information about hand position in one, two, and three dimensions the linear method does miss a large fraction of the variance in the hand position signal. There is a great deal of variability across our experiments (Table 1), some of which can be at least partially accounted for by differences in the time we observed the cells and in the total number of cells observed (see our Fig. 16*B* and Wessberg et al.'s Figs. 2, *E* and *F* and 3, *F* and *G*; note that they plotted *r* in these figures, not *r*^{2} as in ours). However, even in our best experiments, when we compare our results to those from other cortical data sets with different associated behavioral or sensory signals, typically consisting of fewer neurons (e.g., Buracas et al. 1998; Humphrey et al. 1970), we were able to account for a perhaps surprisingly low percentage of the available variance (and almost none >0.5 Hz). The rate of information extracted from the population by linear estimation reaches a maximum of approximately 1 bit/s in the experiments analyzed here. See Rieke et al. (1997) and references to compare this finding with sensory coding, where information rates are 2 orders of magnitude larger than for MI.

This low information rate might simply mean that linear estimation is not an optimal solution for reconstruction of hand motion. However, there are a few reasons to believe that nonlinear estimators will not drastically outperform linear ones. First, our information-theoretic results above indicate that these cells contain a limited amount of information about position or velocity; no nonlinear operation can extract nonexistent information (Cover and Thomas 1991). Second, Wessberg et al. were unable to find a neural network that performed significantly better than the simple linear estimator. There are more sophisticated nonlinear methods for signal estimation given neural activity; work on the application of more elegant Bayesian and/or recursive estimators, which require the development of explicit models of the information encoding process, is in progress (Gao et al. 2002). However, in the absence of any robustly superior nonlinear estimator, we must provisionally conclude that information about hand position is only weakly (perhaps only indirectly) encoded in the activity of these MI neurons. It seems plausible that a more indirect reconstruction approach, using a model that incorporates the dependency of neural firing rate on joint angle, might significantly improve the quality of the achievable hand position reconstructions; this experiment has yet to be carried out.

The linear technique fails to capture effectively all of the variance above 0.5 Hz (Fig. 15). This does not rule out the hypothesis that MI cells encode higher-frequency information about other aspects of the kinematic behavior of the arm, including joint angles and muscle activity (see, e.g., Lemon 1988 for clear examples of higher-frequency EMG information available in cortical spike trains). Nor does it rule out the idea that higher-frequency information is contained in the spike trains, but is not being extracted by the linear algorithm, although as above, we can argue that this is unlikely. Our results here seem to disagree with those of Wessberg et al. who show results of a frequency coherence analysis in their Fig. 2, *A* and *B*. Their cells appear to be significantly coherent with the position of the hand at frequencies ≤5 Hz, and none of their cells appears to be significantly correlated with frequency components below 0.5 Hz, where our signal-to-noise ratio is greatest. It is unclear what accounts for these discrepancies; they could be related to the particular tasks, or to analytical or statistical differences. This apparent absence of high-frequency hand position has important consequences both for the design of neural prosthetic systems and for our understanding of general neural coding in MI.

## APPENDIX

### Method for assessing neural stationarity

We tested neural activity for trends in both *1*) firing rate over the course of each experiment and *2*) firing rate across trial time. The firing rate as a function of time (intratrial or across the experiment) was fit by a line and the slope was tested to see whether it was significantly different from zero. This was done through a bootstrap procedure, described below. Tests were done separately for each cell.

For trends across experimental time the total number of spikes in the tracking phase of each trial was used. A line was fit to these spike counts as a function of experimental time to determine a slope. To test whether this slope was significantly different from zero we performed 400 bootstrap reshufflings. For each iteration of the bootstrap the relationship between the experimental time and the firing rate at that time was broken and randomly shuffled. A line was fit for each iteration to determine a slope, and a distribution of the absolute values of these slopes was constructed. If the absolute value of the slope of the rate trend for a cell was >95% of the randomly shuffled slopes it was considered significantly different from zero. Otherwise, it was not considered to exhibit a rate trend.

For intratrial rate trends a similar procedure was used. In this case, trials were aligned on the start of the tracking phase and the neural activity was averaged over the ensemble of trials to get a mean spike count in each time bin. Spike count as a function of trial time was then subjected to the bootstrap procedure as described above.

Cells with significant rate trends across experimental time were further tested to determine whether the trend in rate affected the tuning functions we compute. A spatiotemporal tuning function (as described above) was computed for each of these cells using data from the first half of the experiment, and a second tuning function was computed using data from the second half of the experiment. A plane (*Eq. 6*, above) was fitted to each and the gain (slope) and orientation (major axis) were compared. If both of these differed significantly between the first and second halves of the experiment the cell was considered to be nonstationary and was not included in any further analyses. Determination of significance was carried out by generating surrogate data to construct a distribution of gains and orientations. The procedure carried out for each cell was as follows.

*1*) Divide the experiment into two parts so that one part, D1, contains the first half of the trials, and the other part, D2, contains the rest.

*2*) For trials in D1 compute the spatiotemporal tuning function, and fit *Eq. 6* to obtain the planar fit parameters, *a*_{0}, *a*_{1}, and θ_{PD} (offset, gain, and orientation, respectively). This was repeated for trials in D2.

*3*) Generate surrogate data is as follows.

*a*) Using trials in D1, generate tuned Poisson spike counts as follows: For each kinematic sample (θ, ρ), generate a Poisson spike count with parameter λ, given by λ = *a*_{0} + *a*_{1}ρ cos (θ - θ_{PD}).

*b*) Using the Poisson counts from *Step 3a*, the spatiotemporal tuning function was calculated, and *Eq. 6* was fit to determine the offset (), gain (), and orientation () parameters for the surrogate.

*c*) *Steps 3a* and *3b* were repeated, this time using only trials from D2, to determine , , and .

*d*) The difference in the gains between D1 and D2 [i.e., ] and the percentage change in orientation i.e., , were calculated, where *abs*(·) means absolute value.

*4*) *Step 3* was repeated 100 times for each cell to obtain a distribution of gain differences and orientation changes under the Poisson assumption.

A cell was considered nonstationary, and excluded from all further analyses, if its gain difference, or orientation change was >97% of the surrogate values [i.e., *P* < 0.03 (Bonferroni-corrected, *P* < 0.05)].

## Acknowledgments

We thank S. Geman, D. Mumford, E. Bienenstock, S. Shoham, Y. Gao, and M. Black for many interesting discussions and M. Hawken for critical comments on earlier drafts of the manuscript.

GRANTS

This work was supported by Defense Advanced Research Projects Agency (MDA972-00-1-0026) and National Institute of Neurological Disorders and Stroke Grants N01-NS-9-2322 and R01-NS-2-5074. L. Paninski was also supported by a Howard Hughes Medical Institute Predoctoral Fellowship.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.↵* L. Paninski and M. R. Fellows contributed equally to this work.

- Copyright © 2004 by the American Physiological Society