HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Supplemental Figures All Versions of this Article: 97/6/4235    most recent 00095.2007v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by Churchland, M. M. Articles by Shenoy, K. V. Search for Related Content PubMed PubMed Citation Articles by Churchland, M. M. Articles by Shenoy, K. V.

Temporal Complexity and Heterogeneity of Single-Neuron Activity in Premotor and Motor Cortex

Neurosciences Program and Department of Electrical Engineering, Stanford University, Stanford, California

Submitted 29 January 2007; accepted in final form 15 March 2007

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The relationship between neural activity in motor cortex and movement is highly debated. Although many studies have examined the spatial tuning (e.g., for direction) of cortical responses, less attention has been paid to the temporal properties of individual neuron responses. We developed a novel task, employing two instructed speeds, that allows meaningful averaging of neural responses across reaches with nearly identical velocity profiles. Doing so preserves fine temporal structure and reveals considerable complexity and heterogeneity of response patterns in primary motor and premotor cortex. Tuning for direction was prominent, but the preferred direction was frequently inconstant with respect to time, instructed-speed, and/or reach distance. Response patterns were often temporally complex and multiphasic, and varied with direction and instructed speed in idiosyncratic ways. A wide variety of patterns was observed, and it was not uncommon for a neuron to exhibit a pattern shared by no other neuron in our dataset. Response patterns of individual neurons rarely, if ever, matched those of individual muscles. Indeed, the set of recorded responses spanned a much higher dimensional space than would be expected for a model in which neural responses relate to a moderate number of factors—dynamic, kinematic, or otherwise. Complex responses may provide a basis-set representing many parameters. Alternately, it may be necessary to discard the notion that responses exist to "represent" movement parameters. It has been argued that complex and heterogeneous responses are expected of a recurrent network that produces temporally patterned outputs, and the present results would seem to support this view.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
A central question regarding the neural control of movement has been how neural responses relate to movement parameters. On the one hand, it is not clear that there need be any straightforward relationship. Even for a network that directly generates movement, individual neuron responses may bear no straightforward relationship to any movement parameter, even those they "control" (Fetz 1992; Robinson 1992; Todorov 2000). On the other hand, a large body of work has stressed the simple relationship between the measured responses of motor cortex neurons and kinematic parameters such as reach direction and speed. In this view, neural responses "code" movement parameters and therefore do have a straightforward relationship with them. Support comes from studies using a center-outreaching task, during which cortical responses typically show cosine tuning for reach direction (Georgopoulos et al. 1982; Schwartz et al. 1988), perhaps implying that movement direction is coded (although see Mussa-Ivaldi 1988; Sanger 1994; Zhang and Sejnowski 1999). Perhaps surprisingly, most such studies have focused on the spatial, rather than temporal, structure of responses. For example, many studies have examined the spatial reference frame (hand-centered, joint-centered, muscle-centered) of direction tuning (Ajemian et al. 2000; Caminiti et al. 1990; Kakei et al. 1999; Scott and Kalaska 1997; Wu and Hatsopoulos 2006). Temporal structure has been addressed less directly, by constructing time-varying population vectors (Georgopoulos et al. 1988; Schwartz 1993, 1994) and/or by regressing neural activity against time-varying movement parameters (Ashe and Georgopoulos 1994; Fu et al. 1995; see Paninski et al. 2004b for a nonlinear approach). Population vector–based decode methods can be remarkably successful (Moran and Schwartz 1999a; Schwartz and Moran 1999), but combine potentially heterogeneous responses, and do little to elucidate the temporal profiles of single neuron activity. Regression coefficients also yield a very indirect view of temporal structure, and results depend on including the right movement parameters and on their degree of correlation with one another (Stark et al. 2006; Todorov 2000).

The controversy surrounding this topic (Scott 2000a,b) stresses the necessity of viewing, as directly as possible, temporal patterns of activity at the level of single neurons. This has been done in the context of the center-out reaching task (Moran and Schwartz 1999b), but rarely extensively. An exception is Sergio et al. (2005), in which multiphasic patterns of single-neuron activity were observed when a weighted manipulandum was used. Under these conditions, the population vector diverged sharply from the measured kinematics, in contrast to prior studies of unencumbered reaching movements. A key question, which we address here, is whether similar multiphasic responses are displayed when movements are executed in the absence of imposed loads.

Observing temporal structure requires satisfying some methodological constraints. First, one wishes to obtain reasonably many nearly identical repeats of the same movement. Averaging responses across differing movements can produce an unrepresentative mean, and a paucity of trials will ensure insufficient power to resolve temporal detail. Second, it is desirable that movement kinematics be dissociated (in the temporal domain) from the patterns of muscle activity. One can then ask which of the two more closely resembles neural responses. Third, it is desirable to have a task in which movements with different time-courses are executed to the same endpoint/goal. One can then compare how neural responses vary with movement time-course. Fourth, averages across neurons can be quite unrepresentative, such that it is essential to view large quantities of data at the level of single-neuron responses.

To these ends, we recorded from M1 and dorsal premotor cortex (PMd) of monkeys heavily trained on a task that enforces peak reach speed (Churchland et al. 2006a,b). Repeated reaches of the same type had very similar hand paths, durations, and velocity profiles. There was also a clear dissociation between kinematics and time-varying muscle activity. Finally, we could evoke reaches to the same target but with different time-courses. Moderately large numbers of trials per neuron were collected (mean of 426, up to 1,027), so that average responses would be meaningful at the level of single neurons. Inspection of such responses reveals that neurons in M1 and PMd frequently display temporally complex patterns of activity during movement. Responses were also surprisingly heterogeneous: often the pattern displayed by a given neuron was observed for no other neuron in the dataset. As a result, neural responses spanned a much larger space than would be expected if they represented a modest number of movement parameters. Thus the complexity observed by Sergio et al. (2005) is not unique to the use of a weighted manipulandum. Rather, those features are even more striking for rapid unencumbered reaches.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Animal protocols were approved by the Stanford University Institutional Animal Care and Use Committee. Our basic methods and the details of this task have been described previously (Churchland et al. 2006b,c). Briefly, two adult male monkeys (Macaca mulatta; 10 kg) sat in a customized chair with a head restraint and performed a reaching task on a fronto-parallel screen. Figure 1, A and B, shows the task structure and typical behavior. Trials began with the appearance of a 12-mm-diam central spot. After this spot was touched and held for 400–500 ms (randomized), the target appeared. During the subsequent delay-period, the target jittered randomly (2 mm SD). Trials were aborted if the hand moved during this time. At the end of the delay period, target jitter ceased and the central spot disappeared, providing the go cue. At this point, the monkey was required to reach to the target. Reaches were rewarded if they were accurate, with reaction times between 150 and 500 ms. A juice reward was delivered after the target was held for 300 ms, with the next trial beginning a few hundred milliseconds later. Delay period durations varied randomly from 400 to 800 ms, and most experiments also included occasional short-delay catch trials. These details were important when studying delay-period responses (Churchland et al. 2006c) but are largely irrelevant to this study, in which only data after the go cue are analyzed. Monkeys were trained to reach at different speeds depending on the color of the target. Briefly, green targets instructed slow reaches, whereas red targets instructed fast reaches. Success was based on peak speed falling within a window, which depended on target color and distance.

View larger version (23K): [in this window] [in a new window]   FIG. 1. Illustration of the task, behavior, and neural recordings. A: monkeys sat in a primate chair 26 cm from a display. Movements began and ended with the hand touching the display. The hand was typically a few millimeters from the screen while in motion. White trace shows the reach trajectory for a typical trial. Gray boxes plot location and size of touch-spot and target, using the same scale as the reach trajectory. B: timeline for the same trial. Horizontal hand (black) and target (gray) position are plotted at top. The target jittered on first appearing. Cessation of jitter provided the go cue, at which time the central spot was also extinguished. Gray trace at bottom plots hand velocity (computed in direction of target), superimposed on the voltage recorded from the medial deltoid (black, arbitrary vertical scale). Traces end at time of reward. Labels T, G, and M indicate target onset, go cue, and measured movement onset, respectively. Data are for monkey A from a session dedicated to EMG recordings. C: recording sites (1 dot per neuron) for monkey A (gray dots) and monkey B (black dots). A small amount (0–0.3 mm) of random displacement has been added to penetration locations to make it clear when multiple recordings were made starting from the same surface location. The large circle outlines the limits of the implanted cylinder. Lines give locations of 1) spur of arcuate sulcus, 2) precentral dimple, and 3) central sulcus. For monkey A (gray lines), these were measured postmortem. For monkey B (black lines), they are inferred from an MRI scan. Recordings from M1 were often made in the central sulcus, well below penetration entry point. D: responses of 1 example neuron (B68). Response rasters (1 tick per spike) are shown for the 17 trials in which a "slow" reach was made to a 12-cm-distant target at 145°. All data are aligned on time of movement onset. The 17 largely overlapping gray traces plot hand velocity for each trial. Gray dots show time of the go cue for each trial. For this condition, mean reaction time was 328 ms (a typical value), measured from the go cue until movement onset. As the latter was measured when hand velocity reached 15% of its peak, this value slightly overestimates the true reaction time. Black trace at top shows mean firing rate as a function of time, computed from the data below.

All trials were presented using a randomized-block design, with any failed trials re-presented at a random time within the current block before proceeding to the next block. We collected two types of datasets. The direction series used two target distances (7 and 12 cm for monkey A; 5 and 12 cm for monkey B) and seven target directions (5, 50, 95, 140, 185, 230, and 320° for monkey A; 10, 55, 100, 145, 190, 235, and 325° for monkey B). Seven (rather than 8) directions were used because there was always one target location that the monkey could not see through his arm. The entire target array was rotated so that the missing (roughly downward) direction was aligned with the position obscured by the arm. The distance series (monkey B only) used five distances (3, 4.2, 6, 8.5, and 12 cm) in the preferred and antipreferred direction of each neuron, estimated using an initial cursory direction series. If the antipreferred direction was near the location obscured by the arm, a compromise was made, and the antipreferred target was moved slightly to one side or the other. Some neurons did not have a clear preferred direction, and we simply chose an axis that produced modulation.

Neural and EMG recordings

Neural recordings were made using single electrodes as described previously (Churchland et al. 2006b,c). All recordings are of single units (manually isolated during the experiment) that were judged, based on stored waveforms and interspike intervals, to have no significant contamination from nearby units. Recordings were made from M1 and caudal PMd. Surface locations of all penetrations are shown in Fig. 1C. Note, however, that this figure does not adequately represent the fact that a number of recordings that began on the precentral gyrus progressed deep into the gray matter of the central sulcus below.

EMG activity was recorded using acutely placed hookwire electrodes during separate dedicated sessions as described in Churchland et al. (2006c). Recordings were made from the deltoid, biceps brachii, triceps brachii, trapezius, latissimus dorsi, pectoralis, and brachioradialis. Multiple recordings were often made from different sites in the same muscle, for a total of 17 recordings across both monkeys. When averaging across trials, EMG traces were differentiated to remove any baseline, rectified, and smoothed (25 ms SD Gaussian). The median (across repeated reaches of the same type) was taken as a function of time. The median was typically very similar to the mean but was less influenced by outliers.

Criteria for inclusion of neural data

We recorded responses from 189 neurons: 138 for the direction series task and 51 for the distance series task. Essentially all recorded neurons exhibited responses either during the delay-period or around the time of the movement. We wished to concentrate our analysis on activity around the time of the movement and to exclude neurons whose principal response was during the delay. To do so, we created two indices of neural response: a delay index and a peri-movement index. The delay index was equal to the SD, across target conditions, of the mean delay-period firing rate. For the peri-movement index, we also computed the SD of firing rates across target conditions, but as a function of time from 100 ms before movement onset until 300 ms after (this was done after filtering). The peri-movement index was the maximum value across time. To be included in our analysis, we insisted that 1) the peri-movement index be at least twice the delay index and 2) that peri-movement activity be reliably task-related (ANOVA, P < 0.0001 for a main effect of target condition) at the time the peri-movement index was measured (i.e., at the time when firing rates varied most across conditions). A total of 137/189 neurons satisfied these criteria and were included in our analysis. These criteria do not insist on any particular form of tuning. For example, a neuron with bimodal direction tuning would be included in our analysis even though it would not have satisfied criteria that insisted on cosine tuning. The above criteria were designed in part to exclude neurons that exhibited primarily delay-period activity, with little additional peri-movement activity. Such neurons often showed strong changes in firing rate around the time of the movement, but only because firing rates were returning to baseline levels. We did not wish such neurons to be classified as having peri-movement responses. Such neurons are successfully excluded by the first criterion, because the peri-movement index typically did not exceed the delay index. In principal, our criteria would also exclude neurons with temporally modulated responses but no tuning for target direction, distance, or speed. However, no convincing examples of such responses were present in our dataset.

Roughly defining cylinder zero (approximately the middle of the precentral dimple) as the PMd/M1 border, 55% (45%) of analyzed neurons were recorded from M1 (PMd). Note that, although the total number of recorded neurons (189) is somewhat smaller than is typical for studies of M1/PMd, the total quantity of data is not. Rather than recording from as many neurons as possible, we chose to instead characterize the activity of each neuron as completely as possible, using a moderately large number of conditions and trials per condition. The mean number of trials per neuron was 426. For comparison, Moran and Schwartz (1999b) recorded approximately the same number of trials, but from more neurons (1,066) with fewer trials per neuron (40).

Analysis criteria

Trials were included in the analysis if the target was hit accurately and held until the time of reward (98% of saved trials). Trials aborted because the hand moved during the delay (or never moved) were not saved, but these comprised at most a few percent of all trials. Failures with regard to peak reach speed were more common (6–8%) but usually involved peak reach speed being only slightly too fast or too slow. We saw no reason to exclude such trials from analysis, as they formed a continuum with normal behavior and presumably did not reflect a lack of effort on the part of the monkey, but rather the challenging nature of the task. In principle, not rejecting such trials might slightly reduce the difference in response between the two instructed speeds. In practice, this seems not to be a concern: large differences were typical.

Plots of mean firing rate as a function of time were made by convolving spike trains with a Gaussian filter (25 ms SD) and averaging across trials (Fig. 1D). Means were computed after aligning individual trial data to the onset of movement, defined as the time when hand velocity reached 15% of its peak. Mean hand velocity, measured in the direction of the target, was similarly computed with single-trial data aligned to the onset of movement.

Some analyses and presentations were restricted to neurons whose responses had the highest signal-to-noise ratios (SNRs). For each neuron, we measured the maximum and minimum firing rates, FR_max and FR_min, each taken across both time and conditions. We also computed SE_max, the maximum (across time and conditions) standard error of the mean. The SNR was defined as (FR_max –FR_min)/(SE_max + ), where is a small value (3 spikes/s) that prevents division by zero and to some degree biases the SNR in favor of neurons with high firing rates. As would be expected, neurons with high SNRs tended to be those with high firing rates, high trial counts, and strong tuning across time and conditions. A number of other reasonable definitions are possible; our intent was simply to use a sensible objective measure.

Computing the preferred direction

The preferred direction (PD) of each neuron was computed as a function of time. For each time, the mean firing rate was plotted against target direction (7 data points). Data were fit with a cosine (free parameters were phase, amplitude, and DC offset) whose peak was taken as the PD. PDs were computed separately for the two target distances and instructed speeds. Note that, for a cosine fit, there is no bias created by the lack of the eighth (downward) direction, which was omitted because the monkey's outstretched arm obscured that target location. A bootstrap procedure was used to compute the sampling distribution: the expected PD distribution, were it remeasured, given the observed measurement variability. For each target direction, we resampled (with replacement) the original distribution of firing rates and computed a new PD. This process was repeated (100 or 1,000 times depending on the analysis), and 95% CIs were computed from the resulting distribution.

Linear regression and kinematic models

The responses of individual neurons were regressed against two linear models. The velocity-tuned model assumes that firing rates are sensitive to reach velocity and speed

(1)

where x is hand position, b_o is the baseline firing rate, b₁ and b₂ capture the horizontal and vertical sensitivities to hand velocity (and together determine the preferred direction), and b₃ captures any nondirectional sensitivity to speed. The complex model assumes that firing rates are also sensitive to hand position and acceleration

(2)

It is frequently assumed that neurons can have variable leads, or even lags, relative to the parameters they represent. Thus for each neuron we tested a range of t, from 200 (firing rate leading kinematics) to –100 ms (firing rate lagging kinematics), in increments of 10 ms, and took the maximum R². We also constructed versions of these models where the firing rate was not allowed to drop below zero. In these cases, fits to the neural data were found by numerical optimization, which was initialized using the parameters obtained by the regression.

We computed an approximate correction for the R² values to account for the fact that some of the data variance is caused by sampling noise and is not expected to be accounted for by the model. This was done only for the complex model. For each neuron, we used the model fits as the underlying firing rate and generated simulated trials in which spikes were generated according to a Poisson process. We generated the same number of simulated trials as we had recorded real trials for that neuron and computed the mean simulated firing rate (after filtering, just as was done for the real data). We refit the model to these noisy simulated responses and computed the variance unaccounted for, which we termed the expected residual variance. We recomputed the R² of the model fits, factoring out the expected residual variance.

(3)

where _total² is the total data variance to be accounted for, _residual² is the actual residual variance after fitting, and _expected² is the estimated expected residual variance caused by sampling error. When applied to simulated recordings, the above correction was successful. Simulated recordings had Poisson spiking statistics, with the underlying firing rate determined by the complex model. We simulated a population of 108 neurons for which trial-counts, baseline firing rates, and firing rate modulations matched the 108 recorded neurons. When fit with the complex model, the mean R² was 0.82, and the mean corrected R² was 1.02, very close to the ideal value of unity.

Principal component analysis

Principal component analysis (PCA) was applied to the neural responses collected using the direction series, 108 of which passed the criteria for inclusion. The response pattern of each neuron was considered to be a single observation. We considered times from 150 ms before until 300 ms after movement onset, with the mean firing rate sampled every 10 ms, for a total of 46 data points. Each neuron thus contributed a single 1,288-dimensional observation (46 time-points and 28 target conditions). Each of these vectors had its mean subtracted and was normalized so that the maximum value minus the minimum value was one. PCA was performed on the population of 108 such observations, resulting in 107 PCs. Each of these is itself a 1,288-dimensional vector, embodying a component response pattern. For a given set of PCs (e.g., the 1st 10), the proportion of variance accounted for was computed as the cumulative sum of the eigenvalues associated with those PCs.

PCA was also applied to simulated populations of neural responses, generated according to different models. Simulated populations were generated by drawing random parameters (e.g., the PD) for each neuron. However, the baseline firing rate of each neuron and its degree of modulation were set to approximate those of one of the recorded neurons. Simulated populations were always of 108 neurons, with each matched to one of the recorded neurons. The underlying firing rate of each neuron determined the mean of a Poisson spiking process. We collected simulated trials, matching the number of these to the number collected for the matched recorded neuron. Mean firing rates were computed by filtering and averaging, just as was done for the recorded responses.

For each PC, we computed an "internal consistency" metric. This was based on the correlation between the subcomponents of the PC at different target locations. For each target location, we took the correlation between its subcomponent and the subcomponent at the three nearest neighbors (the other distance in the same direction and the 2 other directions at the same distance). We define the internal consistency as the mean correlation. The decision to use the three nearest neighbors is reasonable if somewhat arbitrary. Virtually identical results were obtained if we used the nearest five neighbors.

Linear fits to EMG and kinematic responses

For each muscle, we fit the pattern of EMG activity using a linear sum of neural activity, in which each neuron was assigned a weight. The same weights were used to fit EMG for all target conditions; different weights were used for different muscles. To allow for varying conduction delays, each neuron was allowed to contribute at one or more delays: 20, 40, and 60 ms (effectively, the response of each neuron was convolved with a simple linear filter). This method was considered preferable to attempting to assign one delay to each neuron, because it allowed fits to be made by linear regression (rather than numerical optimization). This is not only faster, but aids comparison of the relative fit quality: linear regression has a single solution, whereas numerical optimization is not guaranteed to find the global minimum. Fits were made only for the direction series, and for a given EMG recording were made using only neural activity recorded from that same monkey.

Fits were made to EMG activity from 170 ms before movement onset until 350 ms after movement onset. Data from all 28 conditions were concatenated to produce a single vector, y, of EMG values (14,588 data points long). The same was done for the neural data, but for each neuron, the data were sampled three times: starting 20, 40, and 60 ms before the time when the EMG activity was taken. Thus for a sample of 100 neurons, the resulting matrix, X, would be 14,588 by 301 (the 2nd dimension is 301, rather than 300, because of the addition of a column of ones to allow a nonzero intercept). A linear regression was used to find the vector of weights, b, such that Xb provided the best fit to y.

For analyses regarding generalization, an initial regression used data from a subset of targets (training). The resulting weights, b, were used to provide a prediction for the remaining targets (generalization). For each muscle recording, this procedure was repeated four times, using different subsets of targets for training and generalization. For each of these repeats, training excluded one of the four direction/speed combinations (e.g., far/fast), and generalization was measured for that combination (7 targets total).

We compared fit and generalization performance to that when EMG data were "mirrored" relative to the neural data. For this procedure, we treated EMG responses for rightward reaches as if they had been recorded for leftward reaches (i.e., they were fit with neural responses during leftward reaches). The full mapping was to take EMG responses for the three responses to the right and exchange them with the three responses to the left. Vertical alignment was not disrupted (e.g., top right mapped to top left), and the direction nearest upward remained unaltered. All these procedures were also applied when fitting kinematic measurements (horizontal and vertical hand velocity and acceleration). To allow the best comparison, kinematic values were extracted from the same datasets as the EMG recordings (note that this yields 4 times as many kinematic recordings as muscle recordings).

Neural network simulation

We simulated a very simple artificial neural network of nine units. Connections were recurrent, with each unit connecting to all others. Each unit's activity ranged from –1 to 1. A 10th "bias" unit, connected to all other units and of constant activity one, was also included. At each time-step, the activity of the network, a, was updated using the following equations

(4)

(5)

(6)

(7)

Equations 4 and 5 derive the synaptic input, s. That input was assumed to be driven by both the magnitude of the weighted presynaptic activity and by its rate of change (sensitivity to the latter being determined by c in Eq. 4). A low-pass filter was applied (Eq. 5). Activity was limited to the range –1 to 1 by a sigmoidal transfer function, T, where T(x) = –1 + 2/(1 + e^–x). A second low-pass filter was applied (Eq. 7) to produce the updated activity, a(t + t). The 10th element of a (the activity of the bias unit) was not updated, but was clamped at 1. Simulations used a 5-ms time step. The network output was the sum of three of the nine units (chosen arbitrarily and not allowed to change during optimization). The connection weights, W, and the initial state of the network, a(t_o), were numerically optimized to produce the desired output: the activity of a given muscle for a given target location.

The parameters controlling the dynamics of individual units, c, ₁, and ₂, were set to 0.04, 47 ms, and 25 ms. These values resulted in simple dynamics at the level of individual units but with the network as a whole being capable of complex behavior. These particular values were not critical—similar results were obtained using different parameters.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Behavior

Behavior on this task has been documented previously (Churchland et al. 2006a,b). Briefly, fast reaches (red targets) had durations of 100–200 ms (depending on distance). Even slow reaches (green targets) were fairly swift, with durations of 200–300 ms. Repeated reaches for the same target/instructed speed were typically very similar. This was true both for reach path and for the velocity profile, allowing the mean firing rate to be computed across nearly identical movements. Figure 1D gives an example (17 slow reaches), for which the mean firing rate (black trace) displays considerable temporal structure. Such structure can be lost when the mean is computed across trials with different velocity profiles. This highlights an important feature of our task: it not only allows comparison of responses across instructed speeds, it also ensures that averages are made across very similar behavioral responses. To allow averages that could reveal high-frequency components of the neural responses, we typically endeavored to collect as many trials as possible. The mean was 16/condition and was occasionally >40/condition. Per neuron, the average was 426 and was occasionally as high as 1,000.

Example responses: direction series

Of the neurons tested using the direction series, 108 of 138 showed robust peri-movement responses (see METHODS) and were included in our analysis. We observed a variety of peri-movement patterns, some simple and some complex. Perhaps most strikingly, there was profound heterogeneity in response patterns across neurons. Figures 2 and 3 show responses of 12 example neurons, selected to show the range of observed neural behavior. Each panel (A–F) corresponds to one neuron; each subpanel to a target location. Mean firing rates for slow and fast reaches (green and red traces, respectively; SE given by trace width) are plotted on top of mean hand velocity traces (gray). The example neurons in Fig. 2, A and B, display simple patterns of activity that are literally textbook (Fig. 14.13 of Bear et al. 1996; Fig. 32.1 of Georgopoulos 1995; Fig. 38–13 of Krakauer and Ghez 2000). Activity is roughly cosine-tuned for direction and scales in magnitude with velocity. These patterns adhere closely to a model in which a neuron's activity is determined by the dot product of its preferred direction with the velocity of the reach (Georgopoulos et al. 1982; Moran and Schwartz 1999b). Conversely, the activity of a population of such neurons could be effectively decoded by multiplying each neuron's activity by its preferred vector and then summing to decode velocity (Georgopoulos et al. 1988; Moran and Schwartz 1999b; Schwartz 1993).

View larger version (25K): [in this window] [in a new window]   FIG. 2. Responses of 6 example neurons tested using the direction series (7 directions, 2 distances, 2 instructed speeds). Each panel plots responses of 1 neuron (its identifier given at center). Each subpanel plots data for 1 target location. Gray traces (1 per instructed speed) plot mean hand velocity. Red and green traces plot mean firing rate for fast and slow reaches. Trace widths show ±SE. Vertical calibration bars indicate 20 spikes/s.

View larger version (28K): [in this window] [in a new window]   FIG. 3. Responses of 6 more example neurons tested using the direction series. Format is the same as in Fig. 2. All vertical calibration bars are 20 spikes/s.

View larger version (27K): [in this window] [in a new window]   FIG. 14. Patterns of EMG activity (thick, light traces) for 4 muscle recordings (A–D). Format is similar to prior figures (e.g., Fig. 2). Thin traces show fits based on neural responses. Fits are based only on neural responses recorded from the same monkey.

Another intriguing response pattern is seen in Fig. 2D. For this neuron, the phase of the responses is advanced for slow versus fast reaches. Less obviously, the phase of the response for fast reaches differs with direction: leading the movement for 230° reaches and lagging for 95° reaches. This effect is unlikely to be caused by sampling noise: SE (width of trace) is modest, and the phase change is evident for both distances. Figure 2E shows that responses could also be temporally complex. In this case, responses are either bi- or triphasic depending on reach direction. For many directions, there was surprisingly little difference in the response for the two reach speeds. Figure 2F shows responses for a neuron whose preferred direction shifted 90° between the two reach speeds (from up for fast to left for slow). Thus the neurons in Fig. 2, C–F, all diverge from the canonical, but in different ways.

Further examples are shown in Fig. 3. The responses in Fig. 3A show a temporally complex pattern, especially for leftward reaches. The responses in Fig. 3B are greater for slow reaches and last considerably longer than the movement itself. The responses in Fig. 3, C and D, show temporally bimodal patterns for some reach directions (55 and 100° for C, 235° for D). The neuron in Fig. 3E roughly obeys straightforward expectations (cosine tuning, larger responses for fast speeds). However, the lag between neural response and movement differs for the two reach speeds (especially for 320°), and for some target locations (230°) the response is larger for the slower reaches. Finally, the responses in Fig. 3F exhibit bimodal direction tuning: large responses are seen for 50 and 230°, with little response for 5 and 185°. This pattern is much more pronounced for fast reaches—tuning for slow reaches is roughly cosine.

Further examples

The responses in Figs. 2 and 3 show that a range of response patterns was observed. Yet those examples do not span the entire range of observed patterns. Frequently, a recorded neuron would exhibit a response pattern that was not observed for any other recorded neuron (even allowing for rotations). Sampling error cannot account for the complexity of the response patterns, their heterogeneity across neurons, or the frequent departures from the canonical. Trial counts were reasonably high, and SE (trace width) was usually quite modest.

Given that response heterogeneity is such a striking feature of the recorded responses, we have taken the unorthodox step of presenting a plethora of further examples. To avoid selection bias, we ranked all neurons by the SNR of their responses (see METHODS). We selected the top 24 neurons (skipping those shown previously: Figs. 2, B and E, and 3, A and E). Figures 4 and 5 show the responses of these 24 neurons. We provide below brief summaries of some (although by no means all) of the notable response features for each of these neurons.

View larger version (41K): [in this window] [in a new window]   FIG. 4. Responses of the 12 neurons tested using the direction series that had the highest signal-to-noise ratios (SNRs; skipping those shown earlier). Format is the same as for Figs. 2 and 3. All vertical calibration bars are 20 spikes/s.

View larger version (39K): [in this window] [in a new window]   FIG. 5. Same as for Fig. 4, showing the response of the 12 neurons with the next highest SNRs. All vertical calibration bars are 20 spikes/s.

We recognize that the presentation of this much raw data is somewhat unorthodox. However, in our view, the best way to truly appreciate the complexity and heterogeneity of the neural responses is to inspect many examples.

Example responses: distance series

We also tested a modest number of neurons using a distance series. Of 51 neurons tested, 29 showed robust peri-movement responses and were included in our analysis. The responses of five of these are shown in Fig. 6. As was the case for the direction series, response patterns were often temporally complex and were heterogeneous across neurons. Even the response pattern in Fig. 6A is not as simple as it first appears. The neural response leads the movement for 10° targets and occurs in phase with the movement for 190° targets. Responses are greater for slow reaches for 190° targets but not for 10° targets. The response patterns in Fig. 6, B–E, are even more complex, often with multiple phases, and with changes across distances/directions that would not have been readily predicted. Nominally, these neurons were tested with the two target directions aligned with and against their preferred directions. However, for some neurons, it is unclear that there is a well-defined preferred direction. The direction that evokes the largest response can depend on distance, instructed speed, and when the response is observed. Again, note that trace width gives the SE, and thus even some very fine features that might be mistaken as measurement noise can be seen to be reliable, especially when shared across multiple conditions.

View larger version (37K): [in this window] [in a new window]   FIG. 6. Responses of 5 example neurons tested using the distance series (2 directions by 5 distances). Each panel corresponds to 1 neuron; each subpanel to 1 target location. Red and green traces plot mean firing rate for fast and slow reaches. Trace widths show ±SE. All vertical calibration bars are 20 spikes/s. Gray traces plot mean hand velocity for the 2 instructed speeds.

It is typically reported that firing rates in M1 and PMd correlate with a variety of kinematic parameters, especially velocity and speed. This was also true of the responses we recorded, despite the wide variety of observed response patterns. We regressed each neuron's firing rate against three kinematic parameters: horizontal velocity, vertical velocity, and speed (Moran and Schwartz 1999b). We assumed that firing rates could either lead or lag kinematics and tested leads up to 200 ms and lags up to 100 ms (in 10-ms increments), taking the value that maximized R². The black bars in Fig. 7 plot the distributions, across the 108 neurons recorded using the direction series, of R² values (A) and best leads/lags (B). Mean values were 0.39 (R²) and 88 ms (firing rate leading the movement).

View larger version (14K): [in this window] [in a new window]   FIG. 7. Correlation of firing rate and movement kinematics. A: distribution of R2 values across the 108 neurons tested using the direction series. Firing rate was regressed against either 3 kinematic parameters (black bars for velocity-tuned model) or 7 kinematic parameters (gray bars for complex model). We took the maximum R2 for each neuron, across a range of temporal offsets between 200 (firing rate leading behavior) and –100 ms (firing rate lagging behavior), tested in 10-ms increments. B: distribution of "best" temporal offsets, at which the highest R2 was obtained.

The above results are consistent with prior reports: firing rates correlate with a variety of kinematic parameters, and the sensitivity to velocity/speed seems particularly large: mean R² values rise only modestly, from 0.39 to 0.46, after the addition of position and acceleration to the regression. Nevertheless, the values of R² we obtained were somewhat lower than reported in some prior work. In particular, Moran and Schwartz (1999b) performed an analysis nearly identical to that in Fig. 7A, in which firing rates were regressed against velocity (horizontal and vertical components) and speed. The median R² obtained for that analysis was 0.65, considerably higher than for our parallel analysis (median R² = 0.37 for the velocity-tuned model). Even more strikingly, they obtained an R² of 0.99 when correlating mean hand velocity with the ensemble nondirectional neural response in M1. These discrepancies prompt concern that our dataset may be fundamentally different from theirs. Perhaps the recorded population of neurons was rather different, given that neurons in PMd were included. Or perhaps the tasks differ in some fundamental way. To address these concerns, we repeated their analysis of the ensemble nondirectional response. Their task used a single distance, so we likewise considered data for one distance (12 cm). Their task did not involve instructed speeds, and considerable trial-to-trial variability was present in reach velocity. To emulate this, we collapsed data across the two instructed speeds. The ensemble nondirectional response was thus the mean firing rate across all neurons, directions, and speeds. Mean reach speed was computed in the same way. Figure 8A plots both the ensemble neural response (gray) and mean reach speed (black), and can be compared directly to Figs. 3 and 5 of Moran and Schwartz (1999b). To emulate their original analysis, the speed trace is truncated shortly after reaching its peak, and the neural response is similarly truncated. The rising edge of the neural response leads the rising edge of the speed trace by 150 ms, in close agreement with the value of 145 ms obtained in their original analysis. Indeed, time-shifting the neural response by 145 ms resulted in a very high correlation with reach speed (R² = 0.95), only slightly less than was obtained in the original study (0.99).

View larger version (15K): [in this window] [in a new window]   FIG. 8. Correlation of ensemble nondirectional response with reach speed. Shown are a truncated analysis (A) and a more complete analysis (B). A: replication of analysis in Moran and Schwartz (1999b), similar to that shown in Figs. 3 and 5 of that study. Ensemble nondirectional response (gray trace) and mean reach speed (black trace) are plotted vs. time (0 is movement onset). Assuming a lead of 145 ms (the ideal lead found in Moran and Schwartz 1999b), R2 was 0.95 (maximum R2 for our data was actually 0.96 at a lead of 154 ms, but we use the 145-ms lead to allow comparison with the original analysis). The correlation was computed for a truncated segment of speed trace (from 120 ms before movement onset to 115 ms after) to emulate the original analysis. The ensemble response was similarly truncated. The ensemble neural response and mean speed were computed across all neurons and directions and also across both instructed speeds (but using only the 12-cm distance). B: similar analysis, but where the speed trace is not truncated and the 2 instructed speeds are considered separately. Assuming (as above) the 145-ms lead from the original analysis, the correlation between speed and neural response was weaker (R2 = 0.57) than for the truncated analysis. The correlation increased only modestly after reoptimizing the lead (maximum R2 = 0.66 at a lead of 120 ms).

A critical and very general point is that averaging, whether across neurons or across reaches of different speeds, tends to inflate the apparent correlation between neural responses and movement kinematics. In particular, the average neural response tends to appear more and more Gaussian (and thus more and more like the velocity profile) when it is made across more neurons and/or conditions. The analysis shown in Fig. 8 shows one manifestation of this effect. We also repeated the analysis in Fig. 7, after collapsing data (both firing rates and kinematics) across the two reach speeds. This emulates a typical situation, in which target location, but not reach speed, is experimenter controlled. After this manipulation, the mean R² value for the velocity-tuned model rose from 0.39 to 0.52 (0.55 if we included only the greater distance). Thus the averaging of data across similar but nonidentical behavioral responses can make the relationship between neural activity and kinematics seem rather stronger than it actually is.

Model fits

Inspection suggests that the responses of many individual neurons would be poorly fit by the canonical velocity-tuned model (cosine-tuned for velocity with a nondirectional speed component) or even by the more complex model (also sensitive to position and acceleration). The analysis in Fig. 7A supports this impression. Linear regression is equivalent to fitting the data with a linear model (i.e., a model with linear sensitivities to the different kinematic parameters). The resulting R² values, which were typically well below 1 (means of 0.39 and 0.46 for the 2 models), argue that the data are not well explained by either model. However, interpretation is complicated by two factors. First, the linear model has an obvious flaw: its firing rate can drop below zero, whereas real firing rates cannot. Second, some of the variance we are asking the model to account for is merely sampling error, caused by the finite number of trials recorded from neurons with noisy spiking statistics. Thus even an accurate model would not be expected to account for all the data variance. Inspection of individual examples indicates that the above factors cannot be the sole reason that the kinematics-based models fail to provide good fits. For example, the mean firing rates of neuron B68 (Fig. 2E) possess a tight SE (trace widths), and, although they fall to near-zero briefly, are not strongly truncated at zero. Concerns regarding the above factors should thus be minimized, yet R² values of the model fits were still modest for this neuron (0.49 and 0.51 for the velocity-tuned and complex models, respectively).

However, at the level of the population, it seems likely that linear regression underestimates the quality of the fits that could be provided by the kinematic models. To assess how large this underestimation might be, we constructed a version of the complex model where the firing rate could not drop below zero. This (mildly) nonlinear model was fit to the data. R² values for the model fit were slightly higher than for the linear regression: the mean was 0.48, up from 0.46. We applied an approximate correction to account for the fact that some of the variance is caused by sampling noise and is not expected to be accounted for by the model (see METHODS). The mean corrected R² was 0.56 (for the complex model). This is higher than the uncorrected value of 0.48 but only modestly so. Thus most of the variance unaccounted for by the model is related to real aspects of the neural responses rather than sampling noise.

An important note is that the way in which the model fails is different for each neuron. The responses of neuron A46 (Fig. 2C) are poorly captured because they are larger for the slower reaches (corrected R² = 0.37 for the complex model). The responses of neuron A56 (Fig. 2D) are poorly captured because of the phase difference between the instructed speeds (corrected R² = 0.28 for the complex model). The responses of neuron B68 (Fig. 2E) are poorly captured because of their multiple phases (corrected R² = 0.55 for the complex model). The responses of neuron B107 (Fig. 2F) are poorly captured because of the change in preferred direction between instructed speeds (corrected R² = 0.36 for the complex model). Thus the model of kinematic tuning does not have any single obvious failing. Rather, the variability of neural responses makes them difficult to capture with a simple model. This issue is addressed further in the next section.

As a final point of comparison, we fit the recorded EMG responses using the complex model. The mean R² was 0.49 (uncorrected), which was not terribly different from the mean uncorrected R² for the neural data (0.48). Thus the ability of the complex model to partially fit the neural responses does not necessarily imply that the neural responses code those kinematic features. EMG activity presumably does not code kinematics, even if it can be roughly fit by the complex model.

Response heterogeneity and dimensionality

Not only did the recorded response patterns often differ from the canonical, they often differed sharply from one another. This suggests it may be difficult to adequately capture neural responses using any simple model (kinematics-based or otherwise). However, perhaps this is not the case—perhaps there are only a few basic response patterns, with the apparent heterogeneity arising when these are combined in differing proportions for different neurons. This issue can be addressed by studying the dimensionality of the neural responses using PCA. This approach asks whether we can reconstruct the response pattern of a given neuron from a linear sum of basic response patterns (the PCs). As an example, consider an idealized population of neurons that are cosine tuned for velocity. The response of any neuron could be reconstructed as a linear combination of two PCs: one capturing sensitivity to horizontal velocity and one capturing sensitivity to vertical velocity (a neuron tuned to 45° would use both equally). Thus all the variability in the data would be accounted for by only two PCs.

Figure 9 shows a more realistic version of this scenario. Shown in Fig. 9A are the responses (light traces) of one simulated neuron, taken from a population of 108 simulated neurons. Underlying firing rates were generated by the velocity-tuned model: cosine tuning for velocity plus a nondirectional component related to speed. PDs and gains were randomly assigned for each neuron, but the baseline firing rate and overall modulation were set to match that of 1 of the 108 recorded neurons (t, the lead between neural responses and kinematics, was fixed at 100 ms). We recorded the same number of trials as for the matched real neuron (using a Poisson spiking model). Responses were filtered and means taken, as for the actual neural responses. We extracted 107 PCs, rank ordered by the amount of data variance accounted for. Of those 107 PCs, the first 10 capture 82% of the variance. The first 10 PCs do not capture all of the variance, despite the low dimensionality of the underlying model, because the conversion of underlying rates to spike trains introduces sampling noise and a threshold at zero. However, reconstructions (Fig. 9A, dark traces) based on the first 10 PCs were typically reasonable and largely relied on the first few PCs, which captured the underlying sensitivities built into the model. For example, the reconstruction in Fig. 9A relies almost entirely on the first five PCs (see black histogram). These PCs are shown in Fig. 9B (ordered given their importance for this simulated neuron). PCs 1, 2, and 3 capture, respectively, sensitivity to vertical velocity, horizontal velocity, and speed (the 3 sensitivities of the model itself). These simulations show that PCA can extract small numbers of patterns on which diverse responses are based.

View larger version (13K): [in this window] [in a new window]   FIG. 9. A simulated response pattern and its reconstruction using principal component analysis (PCA). A: light-colored traces plot responses of 1 simulated neuron [106, preferred direction (PD) of 113°] from a population of 108 simulated neurons, with firing rates that obeyed the velocity-tuned model. Each subpanel corresponds to 1 target location, and red and green traces plot responses for instructed-fast and instructed-slow conditions. Responses are filtered averages of spike trains from 26 simulated trials. This trial count was matched to 1 of the recorded neurons, as were baseline firing rate and overall modulation. Calibration bar is 250 ms (roughly the duration of a "slow" reach). Its start is aligned to movement onset. Thin, dark traces plot reconstruction based on 1st 10 PCs. Histogram shows absolute value of PC scores (i.e., weights) used for that reconstruction. B: form of the 5 PCs with largest scores, in order of those scores (scores shown by bars) for this simulated neuron. Labels indicate PC number at level of population (e.g., PC1 is the PC that accounts for the most variance across all simulated neurons). The weighted sum of these 5 PCs provides a reconstruction nearly identical to that shown by dark traces in A (although not quite identical, as PCs 6–10 also make a small contribution).

View larger version (17K): [in this window] [in a new window]   FIG. 10. Response pattern of an example neuron and its reconstruction using PCA (same format as Fig. 9). A: light-colored traces plot responses of neuron B68 (same neuron as in Fig. 2E). Responses are shown only over the range of times used by PCA. Thin dark traces plot reconstruction based on the 1st 10 PCs. Histogram shows absolute value of PC scores used for reconstruction. B: form of the 5 PCs with the largest scores (for this neuron) in order of those scores.

View larger version (14K): [in this window] [in a new window]   FIG. 11. Assessing relevance of PCs through the variance they account for (A) and their internal consistency (B). A: black circles plot residual variance (across all neurons) vs. number of PCs used in reconstruction. PCs were given the usual ordering, with PC1 accounting for the most variance and PC107 accounting for the least. Thus the datapoint at "10" indicates residual variance after using the best 10 PCs. To aid viewing, points are shown only for PCs 1–40 (higher PCs account for very little variance). This analysis is also shown for populations of simulated neurons, generated according to different models (gray symbols). B: black circles plot the internal consistency metric for each PC. Gray symbols do the same for the complex and velocity-tuned models.

The above analysis assesses dimensionality in terms of the percent of variance accounted for by a given number of PCs. In doing, so, we assume (probably reasonably) that those PCs that account for little of the measured variance (the higher-numbered PCs with lower eigenvalues) are primarily capturing sampling noise. However, if we consider the form of the PCs themselves, there is another way in which we can attempt to gauge which PCs are capturing "real" variance and which are capturing only sampling noise. Consider the PCs shown in Fig. 9B, which were extracted from the population of simulated responses. Given that we know the model was constructed to be sensitive to vertical velocity, horizontal velocity, and speed, it is easy to recognize that these features are roughly captured by PCs 1, 2, and 3. However, even if we did not know the form of the model, we would suspect that PCs 1–3 capture real response features. Whatever our prior expectation regarding the features for which the model is tuned (kinematic variables, dynamic variables, muscle activity, etc.), we expect that those features do not change arbitrarily quickly, either across time or across target locations. Thus PCs 1–3 seem quite reliable: values change smoothly with time, and subcomponents at nearby target locations look similar to one another. In contrast, PCs 4 and 5 seem somewhat marginal. There is still some smoothness with time (unsurprising, because mean firing rates were filtered), but subcomponents at different target locations resemble one another only weakly. Thus, for the velocity-tuned model only the first three PCs appear reliable, in agreement with the construction of the model itself, which had three basic sensitivities (horizontal velocity, vertical velocity, and speed). However, for the neural responses, even the higher-number PCs seem reliable (e.g., PCs 7 and 9 in Fig. 10B).

To capture the above impressions quantitatively, we computed the internal consistency of each PC: the correlation between subcomponents at neighboring target locations (see METHODS). Figure 11B plots this metric for each PC for the neural responses (black circles) and for the velocity-tuned and complex versions of the model (light squares and dark triangles, populations simulated as above). In all three cases, internal consistency was high for the lower number PCs (those that accounted for large amounts of the data variance). However, for higher number PCs, internal consistency was much greater for the neural responses than for either model. For example, the mean internal consistency of PCs 20–30 was 0.41 for the neural data, 0.04 for the complex model, and 0.00 for the velocity-tuned model. Thus for the neural responses, these high-number PCs seem to be capturing some real component of the responses that varies smoothly across target locations. For the models, these high-number PCs seem to be capturing primarily sampling noise, which is uncorrelated across target locations.

As a reference, the mean internal consistency of the raw neural responses themselves was 0.66. Choosing an arbitrary but reasonable threshold of one half that value, we can roughly assess dimensionality by asking how many PCs have values higher than the threshold. For the velocity-tuned model, three PCs had values exceeding threshold, in agreement with the three basic sensitivities of the underlying model. For the complex model, 10 PCs had values exceeding threshold, also in rough agreement with the underlying model (the model was sensitive to only 7 kinematic parameters, but also had a variable latency that cannot be captured by a single linear dimension). For the neural responses, there were 29 PCs with values exceeding threshold.

The dimensionality of the model responses cannot be made to match that of the neural responses by adding more measurement noise. We have tested this by sampling fewer simulated trials and reducing the model gain, both of which reduce the SNR (data not shown). For the analysis in Fig. 11A, this merely elevates the later baseline for higher-number PCs; the initial rapid decline is preserved. This manipulation also does not increase the dimensionality of the data when assessed as in Fig. 11B. If anything, dimensionality is slightly reduced (as there is now more noise and less signal to be captured).

The analyses in Fig. 11 reveal that the dimensionality of the neural responses is quite high. Not only are the responses of individual neurons heterogeneous (something readily observed by inspecting individual examples), but they cannot be captured as the sum of some small number of component sensitivities. A caveat is that, although dimensionality is high in a linear space, a nonlinear model could in principle capture the range of observed responses using only a small number of free parameters. Whether this could actually be accomplished using a plausible model is unclear, and is certainly beyond the scope of this study. What can be concluded from this analysis is that any model that considers responses to be driven by the sum of sensitivities to individual movement parameters (kinematic, dynamic, or otherwise) will fail to capture the observed heterogeneity of neural responses unless a rather large number of such parameters is included.

A corollary is that, given the high dimensionality of the neural responses, the space spanned by tho