| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Center for Neural Science, New York University, New York, New York 10003; 2Department of Neuroscience, Brown University, Providence, Rhode Island 02912; and 3Committee on Computational Neuroscience, University of Chicago, Chicago, Illinois 60637
Submitted 22 July 2003; accepted in final form 8 September 2003
 |
ABSTRACT |
|---|
|
|
|
INTRODUCTION |
|---|
|
; 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 |
|---|
|
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 x 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. 2A).
|
|
|
; 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.
|
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
 | (1) |
 | (2) |
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 {ti} 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 x 0.7 cm/s (velocity), and 0.4 radians x 0.5 cm (position). We then calculated the mean firing rates at {ti - } 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.
|
)
 | (1) |
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) |
) 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) |
) 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).
|
): our estimate R (for reconstruction) of the position at the current time t is given by a linear combination
 | (4) |
). 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 (NTN) that can be fairly large [matrix size = D2, where D = 1 + C(Tpre + Tpost)/dt]; 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 r2 = 1 - {E[(R - S)2]/E(S2)}, 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 
, 
(), were regressed onto the coefficients of position, 
(), to obtain the reconstruction of S at , 
(). Goodness of reconstruction was plotted as the SNRs obtained at each frequency
 | (5) |
; Rieke et al. 1997).
|
, which gives the expected r2 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 r2 computed by cross-validation tends to be of lower magnitude than the E(r2) calculated here; therefore we normalize the curves in Fig. 16 by the maximal observed E(r2).
|
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 |
|---|
|
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 1A shows hand paths for trials to each of the 8 directions. This task design results in strong dependencies between horizontal and vertical position (Fig. 1A) and horizontal position and velocity (Fig. 1B). 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. 1C), 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 2A 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. 2A, 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. 3A), and in the power spectrum of hand position (Fig. 3B); most of the power in the hand position signal was below 1 Hz (Fig. 3B; 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. 3C); 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. 4C and 1C) 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. 1D and 4D). Figure 4D 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. 1D and 4D is meant to show that the center-out task induces rate nonstationarities, whereas the PTT does not.
Neural activity during tracking
Figure 2D 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. 1D and 4D). 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.
|
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).
|
) and direction (); = 0 corresponds to movement to the right. The cell shown in Fig. 6A is approximately sinusoidally (i.e., cosine-) tuned for direction [i.e., the function Nv(, , ) 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
 |
 | (6) |
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 a1 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. 6A 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. 6A3) 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. 6B) 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. 6A 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 R2 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.
|
), 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) |
) is a smooth function of , with a maximum at 100 ms, such that a1() 
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., a1) 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 .
|
|
, ) 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 (a1 in Eq. 7) and major axis (PD) is shown for both position (Fig. 9A) and velocity (Fig. 9B). 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. 9C) 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. 9C) 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. 10A), 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.
|
, 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 13 graphically illustrates position information (Fig. 13A) or velocity information (Fig. 13B) 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 (R2 > 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.
|
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 Tpre before and Tpost 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 r2 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., Tpost > 0) robustly improves the reconstructions (Wilcoxon paired-sample rank test, P < 0.05), as expected given the results in Fig. 11E; 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 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 r2 values were normalized by the peak r2 observed during the given experiment, so that these curves range from 0 to 1. Figure 16A 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 Tpre versus r2 graph is quite sharply peaked near zero, indicating a kind of "diminishing returns" in Tpre: when Tpre is small, we have to observe relatively fewer neural data to achieve a given increase in r2 than when Tpre is already large. Increasing the number of neurons considered also improves reconstruction (Spearman rank-order correlation coefficient between number of cells observed and r2; P < 0.05). However, the degree of improvement depends on which population of cells is observed (Fig. 16B). 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 |
|---|
|
) 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. 16B). 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. 16B and Wessberg et al.'s Figs. 2, E and F and 3, F and G; note that they plotted r in these figures, not r2 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 |
|---|
|
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, a0, a1, 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 = a0 + a1 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 |
|---|
|
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 |
|---|
* L. Paninski and M. R. Fellows contributed equally to this work. 
Address for reprint requests and other correspondence: M. R. Fellows, Department of Neuroscience, Brown University, Box 1953, Providence, RI 02912 (E-mail: mrf{at}brown.edu).
|
|
REFERENCES |
|---|
|
Brown EN, Frank LM, Tang D, Quirk MC, and Wilson MA. A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. J Neurosci 18: 7411-7425, 1998.
Buracas GT, Zador AM, DeWeese MR, and Albright TD. Efficient discrimination of temporal patterns by motion-sensitive neurons in primate visual cortex. Neuron 5: 959-969, 1998.
Cover T and Thomas J. Elements of Information Theory. New York: Wiley, 1991.
DeAngelis GC, Ghose GM, Ohzawa I, and Freeman RD. Functional micro-organization of primary visual cortex: receptive field analysis of nearby neurons. J Neurosci 19: 4046-4064, 1999.
Donoghue JP, Sanes JN, Hatsopoulos NG, and Gaal G. Neural discharge and local field potential oscillations in primate motor cortex during voluntary movements. J Neurophysiol 79: 159-173, 1998.
Fu QG, Flament D, Coltz JD, and Ebner TJ. Temporal encoding of movement kinematics in the discharge of primate primary motor and premotor neurons. J Neurophysiol 73: 836-854, 1995.
Georgopoulos AP, Caminiti R, and Kalaska J. Static spatial effects in motor cortex and area 5: quantitative relations in a two-dimensional space. Exp Brain Res 54: 446-454, 1984.[Web of Science][Medline]
Georgopoulos AP, Kalaska JF, Caminiti R, and Massey JT. On the relations between the direction of two-dimensional arm movements and cell discharge in primary motor cortex. J Neurosci 2: 1527-1537, 1982.[Abstract]
Georgopoulos AP, Lurito JT, Petrides M, Schwartz AB, and Massey JT. Mental rotation of the neuronal population vector. Science 243: 234-236, 1989.
Georgopoulos AP, Schwartz AB, and Kettner RE. Neuronal population coding of movement direction. Science 233: 1416-1419, 1986.
Haag J and Borst A. Active membrane properties and signal encoding in graded potential neurons. J Neurosci 18: 7972-7986, 1998.
Hatsopoulos NH, Ojakangas CL, Paninski L, and Donoghue JP. Information about movement direction obtained from synchronous activity of motor cortical neurons. Proc Natl Acad Sci USA 95: 15706-15711, 1998.
Humphrey DR, Schmidt EM, and Thompson WD. Predicting measures of motor performance from multiple cortical spike trains. Science 170: 758-762, 1970.
Johnson MTV, Coltz JD, Hagen MC, and Ebner TJ. Visuomotor processing as reflected in the directional discharge of premotor and primary motor cortex neurons. J Neurophysiol 81: 875-878, 1999.
Kakei S, Hoffman DS, and Strick PL. Muscle and movement representations in the primary motor cortex. Science 285: 2136-2139, 1999.
Kalaska JF, Cohen DAD, Hyde ML, and Prud'homme MA. Comparison of movement direction-related versus load direction-related activity in primate motor cortex using a two-dimensional reaching task. J Neurosci 9: 2080-2102, 1989.[Abstract]
Kennedy PR and Bakay RA. Restoration of neural output from a paralyzed patient by a direct brain connection. Neuroreport 9: 1707-1711, 1998.[Web of Science][Medline]
Kettner RE, Schwartz AB, and Georgopoulos AP. Primate motor cortex and free arm movements to visual targets in three-dimensional space. III. Positional gradients and population coding of movement direction from various movement origins. J Neurosci 8: 2938-2947, 1988.[Abstract]
Kolmogorov AN. Information Theory and the Theory of Algorithms, edited by Shiriaev AN. Boston, MA: Kluwer Academic, 1993.
Kowalski N, Depireux DA, and Shamma SA. Analysis of dynamic spectra in ferret primary auditory cortex. J Neurophysiol 76: 3524-3534, 1996.
Lemon R. The output map of the primate motor cortex. Trends Neurosci 11: 501-506, 1988.[CrossRef][Web of Science][Medline]
Mainen ZF and Sejnowski TJ. Reliability of spike timing in neocortical neurons. Science 268: 1503-1506, 1995.
Maynard EM, Hatsopoulos NG, Ojakangas CL, Acuna BD, Sanes JN, Normann RA, and Donoghue JP. Neuronal interactions improve cortical population coding of movement direction. J Neurosci 19: 8083-8093, 1999.
Moran DW and Schwartz AB. Motor cortical representation of speed and direction during reaching. J Neurophysiol 82: 2676-2692, 1999a.
Moran DW and Schwartz AB. Motor cortical activity during drawing movements: population representation during spiral tracing. J Neurophysiol 82: 2693-2704, 1999b.
Moxon KA, Markowitz RS, Nicolelis MAL, and Chapin JK. Realtime control of a robot arm using simultaneously recorded neurons in the motor cortex. Nat Neurosci 2: 664-667, 1999.[CrossRef][Web of Science][Medline]
Mussa-Ivaldi FA. Do neurons in the motor cortex encode movement direction? An alternative hypothesis. Neurosci Lett 91: 106-111, 1988.[CrossRef][Web of Science][Medline]
Neter J, Wasserman W, and Kutner MH. Applied Linear Statistical Models. Homewood, IL: Irwin, 1985.
Nordhausen CT, Maynard EM, and Normann RA. Single unit recording capabilities of a 100 microelectrode array. Brain Res 726: 129-140, 1996.[CrossRef][Web of Science][Medline]
Paninski L, Fellows M, Hatsopoulos N, and Donoghue J. Coding dynamic variables in populations of motor cortex. Soc Neurosci Abstr 6: 665.9, 1999.
Porter R and Lemon R. Corticospinal Function and Voluntary Movement. Oxford, UK: Clarendon Press, 1993.
Press WH, Teukolsky SA, Vetterling WT, and Flannery BP. Numerical Recipes in C. Cambridge, UK: Cambridge Univ. Press, 1992.
Pugh M, Ringach D, Shapley R, and Shelley M. Computational modeling of orientation tuning dynamics in monkey primary visual cortex. J Comput Neurosci 8: 143-159, 2000.[CrossRef][Web of Science][Medline]
Rieke F, Warland D, de Ruyter van Steveninck R, and Bialek W. Spikes: Exploring the Neural Code. Cambridge, MA: MIT Press, 1997.
Ringach DL, Hawken MJ, and Shapley R. Dynamics of orientation tuning in macaque primary visual cortex. Nature 387: 281-284, 1997.[CrossRef][Medline]
Salinas E and Abbott LF. Coordinate transformations in the visual system: how to generate gain fields and what to compute with them. Prog Brain Res 130: 175-190, 2001.[Medline]
Sanes JN and Donoghue JP. Static and dynamic organization of motor cortex. Adv Neurol 73: 277-296, 1997.[Web of Science][Medline]
Sanger TD. Probability density estimation for the interpretation of neural population codes. J Neurophysiol 76: 2790-2793, 1996.
Schwartz A. Motor cortical activity during drawing movements: population representation during sinusoid tracing. J Neurophysiol 70: 28-36, 1993.
Schwartz A and Moran DW. Motor cortical activity during drawing movements: population representation during lemniscate tracing. J. Neurophysiol 82: 2705-2718, 1999.
Schwartz AB, Taylor DM, and Tillery SI. Extraction algorithms for cortical control of arm prosthetics. Curr Opin Neurobiol 11: 701-708, 2001.[CrossRef][Web of Science][Medline]
Scott SH and Kalaska JF. Reaching movements with similar hand paths but different arm orientations. I. Activity of individual cells in motor cortex. J Neurophysiol 77: 826-852, 1997.
Sergio LE and Kalaska JF. Changes in the temporal pattern of primary motor cortex activity in a directional isometric force versus limb movement task. J Neurophysiol 80: 1577-1583, 1998.
Serruya MD, Hatsopoulos NG, Paninski L, Fellows MR, and Donoghue JP. Brain-machine interface: instant neural control of a movement signal. Nature 416: 141-142, 2002.[CrossRef][Medline]
Taira M, Boline J, Smyrnis N, Georgopoulos A, and Ashe J. On the relations between single cell activity in the motor cortex and the direction and magnitude of three-dimensional static isometric force. Exp Brain Res 109: 367-376, 1996.[Web of Science][Medline]
Taylor DM, Tillery SI, and Schwartz AB. Direct cortical control of 3D neuroprosthetic devices. Science 296: 1829-1832, 2002.
Todorov E. Direct cortical control of muscle activation in voluntary arm movements: a model. Nat Neurosci 3: 391-398, 2000.[CrossRef][Web of Science][Medline]
Todorov E. Cosine tuning minimizes motor errors. Neural Comput 14: 1233-1260, 2002.[CrossRef][Web of Science][Medline]
Warland D, Reinagel P, and Meister M. Decoding visual information from a population of retinal ganglion cells. J Neurophysiol 78: 2336-2350, 1997.
Wessberg J, Stambaugh CR, Kralik JD, Beck PD, Laubach M, Chapin JK, Kim J, Biggs SJ, Srinivasan MA, and Nicolelis MA. Real-time prediction of hand trajectory by ensembles of cortical neurons in primates. Nature 408: 361-365, 2000.[CrossRef][Medline]
This article has been cited by other articles:
|
|
 |
|
  K. Ganguly, L. Secundo, G. Ranade, A. Orsborn, E. F. Chang, D. F. Dimitrov, J. D. Wallis, N. M. Barbaro, R. T. Knight, and J. M. Carmena Cortical Representation of Ipsilateral Arm Movements in Monkey and Man J. Neurosci., October 14, 2009; 29(41): 12948 - 12956. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 A. S. Dickey, A. Suminski, Y. Amit, and N. G. Hatsopoulos Single-Unit Stability Using Chronically Implanted Multielectrode Arrays J Neurophysiol, August 1, 2009; 102(2): 1331 - 1339. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
G. H. Mulliken, S. Musallam, and R. A. Andersen Decoding Trajectories from Posterior Parietal Cortex Ensembles J. Neurosci., November 26, 2008; 28(48): 12913 - 12926. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
E. Stark, A. Globerson, I. Asher, and M. Abeles Correlations between Groups of Premotor Neurons Carry Information about Prehension J. Neurosci., October 15, 2008; 28(42): 10618 - 10630. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. M. Radhakrishnan, S. N. Baker, and A. Jackson Learning a Novel Myoelectric-Controlled Interface Task J Neurophysiol, October 1, 2008; 100(4): 2397 - 2408. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 G. H. Mulliken, S. Musallam, and R. A. Andersen Inaugural Article: Forward estimation of movement state in posterior parietal cortex PNAS, June 17, 2008; 105(24): 8170 - 8177. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Truccolo, G. M. Friehs, J. P. Donoghue, and L. R. Hochberg Primary Motor Cortex Tuning to Intended Movement Kinematics in Humans with Tetraplegia J. Neurosci., January 30, 2008; 28(5): 1163 - 1178. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
C. A. Chestek, A. P. Batista, G. Santhanam, B. M. Yu, A. Afshar, J. P. Cunningham, V. Gilja, S. I. Ryu, M. M. Churchland, and K. V. Shenoy Single-Neuron Stability during Repeated Reaching in Macaque Premotor Cortex J. Neurosci., October 3, 2007; 27(40): 10742 - 10750. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. M. Churchland and K. V. Shenoy Temporal Complexity and Heterogeneity of Single-Neuron Activity in Premotor and Motor Cortex J Neurophysiol, June 1, 2007; 97(6): 4235 - 4257. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Wang, S. S. Chan, D. A. Heldman, and D. W. Moran Motor Cortical Representation of Position and Velocity During Reaching J Neurophysiol, June 1, 2007; 97(6): 4258 - 4270. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Deneve, J.-R. Duhamel, and A. Pouget Optimal Sensorimotor Integration in Recurrent Cortical Networks: A Neural Implementation of Kalman Filters J. Neurosci., May 23, 2007; 27(21): 5744 - 5756. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
N. G. Hatsopoulos, Q. Xu, and Y. Amit Encoding of Movement Fragments in the Motor Cortex J. Neurosci., May 9, 2007; 27(19): 5105 - 5114. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. M. Yu, C. Kemere, G. Santhanam, A. Afshar, S. I. Ryu, T. H. Meng, M. Sahani, and K. V. Shenoy Mixture of Trajectory Models for Neural Decoding of Goal-Directed Movements J Neurophysiol, May 1, 2007; 97(5): 3763 - 3780. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
T. N. Aflalo and M. S. A. Graziano Relationship between Unconstrained Arm Movements and Single-Neuron Firing in the Macaque Motor Cortex J. Neurosci., March 14, 2007; 27(11): 2760 - 2780. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. M. Morrow, L. R. Jordan, and L. E. Miller Direct Comparison of the Task-Dependent Discharge of M1 in Hand Space and Muscle Space J Neurophysiol, February 1, 2007; 97(2): 1786 - 1798. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. R. Townsend, L. Paninski, and R. N. Lemon Linear Encoding of Muscle Activity in Primary Motor Cortex and Cerebellum J Neurophysiol, November 1, 2006; 96(5): 2578 - 2592. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. B. Averbeck and D. Lee Effects of Noise Correlations on Information Encoding and Decoding J Neurophysiol, June 1, 2006; 95(6): 3633 - 3644. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. J. Lehmkuhle, R. A. Normann, and E. M. Maynard Trial-by-Trial Discrimination of Three Enantiomer Pairs by Neural Ensembles in Mammalian Olfactory Bulb J Neurophysiol, March 1, 2006; 95(3): 1369 - 1379. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
E. Stark, R. Drori, and M. Abeles Partial Cross-Correlation Analysis Resolves Ambiguity in the Encoding of Multiple Movement Features J Neurophysiol, March 1, 2006; 95(3): 1966 - 1975. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
N. G. Hatsopoulos Encoding in the Motor Cortex: Was Evarts Right After All? Focus on "Motor Cortex Neural Correlates of Output Kinematics and Kinetics During Isometric-Force and Arm-Reaching Tasks" J Neurophysiol, October 1, 2005; 94(4): 2261 - 2262. [Full Text] [PDF]
|
|
|
|
|
|
L. E. Sergio, C. Hamel-Paquet, and J. F. Kalaska Motor Cortex Neural Correlates of Output Kinematics and Kinetics During Isometric-Force and Arm-Reaching Tasks J Neurophysiol, October 1, 2005; 94(4): 2353 - 2378. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. A. Lebedev, J. M. Carmena, J. E. O'Doherty, M. Zacksenhouse, C. S. Henriquez, J. C. Principe, and M. A. L. Nicolelis Cortical Ensemble Adaptation to Represent Velocity of an Artificial Actuator Controlled by a Brain-Machine Interface J. Neurosci., May 11, 2005; 25(19): 4681 - 4693. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Truccolo, U. T. Eden, M. R. Fellows, J. P. Donoghue, and E. N. Brown A Point Process Framework for Relating Neural Spiking Activity to Spiking History, Neural Ensemble, and Extrinsic Covariate Effects J Neurophysiol, February 1, 2005; 93(2): 1074 - 1089. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
B. B. Averbeck, M. V. Chafee, D. A. Crowe, and A. P. Georgopoulos Parietal Representation of Hand Velocity in a Copy Task J Neurophysiol, January 1, 2005; 93(1): 508 - 518. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 G. M. Friehs, V. A. Zerris, C. L. Ojakangas, M. R. Fellows, and J. P. Donoghue Brain-Machine and Brain-Computer Interfaces Stroke, November 1, 2004; 35(11_suppl_1): 2702 - 2705. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
L. Paninski, S. Shoham, M. R. Fellows, N. G. Hatsopoulos, and J. P. Donoghue Superlinear Population Encoding of Dynamic Hand Trajectory in Primary Motor Cortex J. Neurosci., September 29, 2004; 24(39): 8551 - 8561. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
TOP
ABSTRACT
INTRODUCTION
on this log-log scale); the Poisson line fits the data well at rates up to about 10 Hz; variance is slightly sub-Poisson at higher rates.
