Abstract
Previous studies have suggested that complex movements can be elicited by electrical stimulation of the motor cortex. Most recording studies in the motor cortex, however, have investigated the encoding of time-independent features of movement such as direction, velocity, position, or force. Here, we show that single motor cortical neurons encode temporally evolving movement trajectories and not simply instantaneous movement parameters. We explicitly characterize the preferred trajectories of individual neurons using a simple exponential encoding model and demonstrate that temporally extended trajectories not only capture the tuning of motor cortical neurons more accurately, but can be used to decode the instantaneous movement direction with less error. These findings suggest that single motor cortical neurons encode whole movement fragments, which are temporally extensive and can be quite complex.
Introduction
Almost 100 years ago, it was postulated that the motor cortex should be viewed as a synthetic organ for complex motor actions (Leyton and Sherrington, 1917). According to this view, individual motor cortical neurons encode elementary, time-dependent movement fragments, which are combined via the anatomical connectivity within the motor cortex to generate the rich variety of complex motor actions that are ubiquitous in every day life. This view implies that the motor cortex constitutes a sort of language of motor actions where individual motor cortical neurons encode the movement primitives of the language, whereas the grammar of action is defined by the manner in which these neurons combine their activities to generate more complex motor actions.
Despite this early viewpoint that individual motor cortical neurons represent temporally evolving movements, much of the research over the past 40 years since the advent of behavioral electrophysiology has focused on examining time-independent movement parameters, including position (Georgopoulos et al., 1984; Paninski et al., 2004a), velocity (Moran and Schwartz, 1999), acceleration, direction (Georgopoulos et al., 1982), movement distance, joint torques, and hand force (Evarts, 1968; Smith et al., 1975; Hepp-Reymond et al., 1978; Cheney and Fetz, 1980; Kalaska et al., 1989; Taira et al., 1996; Cabel et al., 2001), or some combinations of those parameters (Fu et al., 1993; Kurata, 1993; Ashe and Georgopoulos, 1994; Fu et al., 1995; Stark et al., 2006). Common to most of these studies is the implication that single motor cortical neurons encode a time-independent, instantaneous variable at some fixed time delay.
The dominant viewpoint has been that motor cortical neurons are tuned to particular movement directions (Georgopoulos et al., 1982). However, several studies have shown that directional tuning measured on a fine time scale can vary substantially during the performance of a simple reaching task (Mason et al., 1998; Sergio and Kalaska, 1998; Sergio et al., 2005). These results suggest that motor cortical neurons encode temporally extensive representations of movement that may be quite complex. In addition, previous experiments have demonstrated that electrical stimulation (Graziano et al., 2002) even of a single neuron in motor cortex (Brecht et al., 2004) can elicit complex, time-evolving movements. Here, we explicitly characterize the shape and duration of time-dependent trajectories that appear to be encoded by single motor cortical neurons.
Materials and Methods
Behavioral tasks
Three macaque monkeys (Macaca mulatta) were operantly trained to move a cursor appearing above the monkey's hand location to targets projected onto a horizontal, reflective surface in front of the monkey. The monkey's arm rested on cushioned arm troughs secured to links of a two-joint exoskeletal robotic arm (KINARM system; BKIN Technologies, Kingston, Ontario, Canada) (Scott, 1999) underneath the projection surface. The shoulder joint was abducted 90° such that shoulder and elbow flexion and extension movements were made in the horizontal plane. The shoulder and elbow joint angles were sampled at 500 Hz by the motor encoders of the robotic arm. The x and y positions of the hand were computed using the forward kinematics equations.
In the random target-pursuit (RTP) task, a sequence of seven targets appeared on the projection surface. At any one time, a single target appeared at a random location in the workspace, and the monkey was required to move to it. As soon as the cursor reached the target, the target disappeared and a new target appeared in a new, random location. After reaching the seventh target, the monkey was rewarded with a drop of water.
The loading RTP task was performed by one animal in which viscous loads were applied to both the shoulder and elbow joints (with viscosities of 0.4 or 0.5 N-m-s/rad) on randomly intermingled trials along with unloaded trials. These viscous loads were applied by two motors on the joints of the exoskeletal robot on which the monkey's arm rested, and opposed the motion of the joints with torques proportional to the joint angular velocities.
The center-out task involved movements from a center target to one of eight peripherally positioned targets (5–7 cm distance). On each trial, one of the eight peripherally positioned targets was pseudorandomly selected. The task consisted of two epochs: (1) a fixed instruction period of 200 ms during which the monkey was required to hold its hand over the center target, and one of the eight final targets appeared, and (2) a “go” period during which the target began to blink, informing the monkey to begin moving to the peripherally positioned target. This task was performed immediately after the performance of the RTP task (one data set in one animal) and collected in a single data set so that we could insure that we were recording the same units for both tasks and could make direct comparisons in the temporal shifts of preferred direction.
Electrophysiology
A silicon-based electrode array composed of 100 electrodes (1.0 mm electrode length; 400 μm interelectrode separation) was implanted in the arm area of the primary motor cortex (MI) of each monkey. During a recording session, signals from up to 96 electrodes were amplified (gain, 5000), bandpass filtered between 0.3 Hz and 7.5 kHz, and recorded digitally (14-bit) at 30 kHz per channel using a Cerebus acquisition system (Cyberkinetics, Foxborough, MA). Only waveforms (1.6 ms in duration; 48 sample time points per waveform) that crossed a threshold were stored and spike-sorted using Offline Sorter (Plexon, Dallas, TX). Interspike interval histograms were computed to verify single-unit isolation by ensuring that <0.05% of waveforms possessed an interspike interval <1.6 ms. Signal-to-noise ratios were defined as the difference in mean peak-to-trough voltage divided by twice the mean SD. The mean SD was computed by measuring the SD of the spike waveform over all acquired spikes at each of the 48 sample time points of the waveform and then averaging. All isolated single units used in this study possessed signal-to-noise ratios of 4:1 or higher. A total of 13 data sets (four data sets for animal BO, four data sets for animal RJ, and five data sets from animal RS) were analyzed, where a data set is defined as all simultaneously recorded neural data collected in one recording session. Each data set contained between 34 and 145 simultaneously recorded units from MI. A total of 859 single-unit samples were recorded from MI over all 13 data sets. We use the term “samples” because the recordings were made from chronically implanted arrays in each monkey. Therefore, data collected over recording sessions in each monkey are not necessarily from completely different ensembles of units. Ensembles consisted of “randomly” selected units from MI except for a possible bias for neurons with large cell bodies that would generate higher signal-to-noise ratios. All of the surgical and behavioral procedures were approved by the University of Chicago Institutional Animal Care and Use Committee and conform to the principles outlined in the Guide for the Care and Use of Laboratory Animals (National Institutes of Health publication no. 86–23, revised 1985).
Analysis
PD shifts in the center-out task.
Temporal variations in preferred direction (PD) in the center-out task were computed by aligning the spike trains with respect to movement onset, binning the spikes in 50 ms bins, and fitting the binned spike counts as a function of movement direction with a cosine function. For each neuron, an optimal 400 ms window was selected over which the shifts in PD were compared with those as a function of lead/lag time computed in the RTP task. PDs occurring early in the 400 ms window (e.g., before movement onset) are comparable with PDs computed at lag times in the RTP task (i.e., times when the instantaneous movement direction lags the measured firing rate or, expressed differently, times when the measured firing rate occurs before the instantaneous movement direction). Likewise, PDs occurring late in the 400 ms window are comparable with PDs computed at lead times (i.e., times when the instantaneous movement direction leads or occurs before the measured firing rate). The start of the optimal 400 ms window was found that maximized the mean correlation coefficient between the shifting PDs in the center-out with those as a function of lag/lead time computed in the RTP task. This optimization was done for each neuron separately.
Encoding model.
To build the encoding model for each neuron, we used a generalized linear model that estimated a set of coefficients that maximized the likelihood of observing a spike or not within a small spike sampling window given that the monkey's hand generated a particular velocity trajectory segment. Movement segments (x and y positions sampled every 2 ms) of duration τlead + τlag were extracted every 50 ms throughout successfully completed trials of the RTP task resulting in tens of thousands of sample movement segments in a typical data set (see Fig. 1B). A spike sampling window (δt = 10 ms) was assigned to each sample segment centered at time to (supplemental information and supplemental Fig. 1, available at www.jneurosci.org as supplemental material). For each neuron, if a spike occurred within the spike sampling window, a 1 was associated with the sample movement segment; otherwise, a zero was associated with that sample movement segment.
The x and y trajectories of each sample movement segment were smoothed by replacing each sample point with the time averaged position (±50 ms about the sample point) and then numerically differentiated (segments with peak velocities >100 cm/s were excluded from the model) to generate velocity trajectory vectors of duration τlead + τlag:
where vx(t0 + nΔt) and vy(t0 + nΔt) are the x and y components of the hand velocity, respectively, at time t0 + nΔt where Δt = 2 ms. These vectors were normalized by dividing by the L2-norm to create normalized velocity trajectories:
The full extended model used in this study also included terms reflecting the average speed, v̄t0, and average position, (x̄t0, ȳt0), of the trajectory:
where k⇑ is referred to as the preferred velocity trajectory of the neuron. The preferred path or “pathlet” of the neuron is obtained by integrating k⇑ in time as if it were a vector of velocities. Because k⇑ includes coefficients corresponding to both x and y components of the velocity trajectory, the integration was performed on the x and y components of k⇑ separately.
To reduce the dimensionality of v̂t0 (e.g., 400 dimensions for a 400 ms trajectory length; 200 dimensions for x and y, respectively) we performed a principal component analysis on all sample velocity trajectories (before normalization). This dimensionality reduction was performed to reduce the computational load of estimating the optimal linear coefficients in the generalized linear model. The first 10 principal components accounted for over 90% of the variance, so we assume that k⇑ is a linear combination of these top 10 components (i.e., k⇑ = Uβ⇑, where β⇑ is a 10 dimensional parameter vector (referred to as the transformed preferred velocity trajectory) and U is a matrix of 10 columns each corresponding to a principal component. The 14 free parameters, β⇑ (10 dimensional), a, b, c, γ, were estimated by maximum likelihood using standard methods for generalized linear models. The exponent on the right side of the equation above is denoted the linear predictor. When comparing any of the exponential models to the data, we computed the linear predictor for each sample trajectory. This set of scalar values was then ordered and divided into equal-sized quantiles, or bins. The probability of spiking within each bin was computed from data and shown as a blue point in a graph of the probability versus linear predictor.
Our modeling approach falls within the point process framework used by Truccolo et al. (2005), where the conditional firing intensity, λ(t0 | v̂t0) = P(spike(t) | v̂t0)/δt, depends only on the external motion covariates. In their work, the emphasis was on the dependence of the firing intensity on the firing history of the neuron and other neurons whereas the movement covariate was taken at one single time point. Here, we explore the dependence of the firing intensity on movement segments taken at multiple time points.
Receiver-operating characteristic curves.
To assess the degree to which the encoding model fit the data, receiver operating characteristic (ROC) curves were generated to quantify the relationship between the probability of correctly predicting a spike (hit probability) versus the probability of incorrectly predicting a spike when it was absent (false-positive probability). The curve was generated by systematically varying a threshold for predicting a spike or not. For each threshold value, a spike was assumed to occur if the spiking probability exceeded the threshold. To verify that the model was not over-fitting the data, all ROC analyses were performed on test data that were not used to build the model. We assigned a random 90% of the data (i.e., the training data) to build the encoding model, and the remaining 10% of the data were assigned as test data.
Relationship with previous encoding models.
The complete encoding model that we propose is novel because it combines the normalized velocity trajectory, mean speed, and mean hand position in a multiplicative manner because these terms add in the exponent. However, the model can be reduced to a cosine tuning model of direction (Georgopoulos et al., 1982) and linear tuning model of speed (Moran and Schwartz, 1999) by ignoring position and making two simplifying assumptions. By linearizing the exponential and considering a velocity trajectory composed of the instantaneous velocity at only one time point, v̂t0 = [vx(t0 + Δt), vy(t0 + Δt)], where vx(t0 + Δt) is the instantaneous x velocity and vy(t0 + Δt) is the instantaneous y velocity at time point t0 + Δt:
where v̄t0 is the average speed, θ is the instantaneous direction of the hand and θpd is the direction of k⇑, which corresponds to the preferred direction of the neuron.
Likewise, Cartesian planar tuning of position (Georgopoulos et al., 1984) can be shown by ignoring the velocity trajectory and the speed and linearizing the exponential with respect to the position
Torque encoding model.
We applied the same encoding framework to generate the preferred torque trajectories. Instead of velocity trajectories, we extracted sample torque trajectory vectors:
where qsh(t0 + nΔt) and qelb(t0 + nΔt) are the shoulder and elbow components of the joint torque, respectively, at time t0 + nΔt where Δt = 2 ms. These sample torque vectors were normalized in the same manner as in the kinematic model. The average speed term v̄t0, used in the kinematic model, was replaced with an average torque magnitude term, q̄t0.
Results
Evidence for trajectory encoding
As other studies have found (Mason et al., 1998; Sergio and Kalaska, 1998; Sergio et al., 2005), we observed systematic shifts in PD among multiple single units in MI recorded simultaneously (n = 68 units) while monkeys performed an instructed-delay center-out task (Fig. 1A). The PDs were measured in small, nonoverlapping time intervals (50 ms in duration) relative to movement onset and could exhibit shifts 90° or larger over a 400 ms period of time.
Evidence for trajectory encoding in the motor cortex. A, Temporal evolution of preferred directions (in 50 ms bins) for four MI neurons relative to movement onset computed in an instructed-delay center-out task to one of eight targets. B, Left, Example of a single successful behavioral trial in the RTP task starting with the blue target, proceeding through the green targets, and ending at the red target. A sample trajectory (400 ms long) is also shown (blue and red lines), which was used to compute the preferred direction shifts as a function of lead/lag times relative to the measured firing rate of a neuron. This sample trajectory was also used to build the exponential encoding model. The crosshairs represent the mean position of the sample trajectory. Right, Fifty example movement segments extracted from an RTP task data set. C, Temporal evolution of preferred directions for the same four neurons as in A as a function of lead/lag time relative to the measured firing rate computed in the RTP task. The red and blue arrows represent preferred directions computed for movements lagging or leading the measured firing rate of the neuron (50 ms bins), respectively. D, The preferred directional paths computed in the RTP task (red and blue arrows) and in the center-out task (black arrows) for the four neurons by vectorally adding the preferred direction vectors from right to left. A–C, Preferred directions that fit a cosine function (r2 > 0.25 for 16 directions in RTP task and r2 > 0.5 for 8 directions in center-out task) are shown in solid black. Those that do not fit a cosine function appropriately are shown as dotted lines. It should be noted that a threshold r2 of 0.25 for 16 directions and a threshold r2 of 0.5 for eight directions are statistically significant but somewhat lower than the thresholds used by other studies (Georgopoulos et al., 1982).
Shifts in PDs were also observed when monkeys performed an RTP task with their arm. Instead of repeating a limited set of movements in one of eight prescribed directions, we had the animal generate a rich variety of trajectories and paths with different spatial shapes, velocities, and positions akin to white-noise experiments conducted in sensory systems to probe the tuning properties of sensory neurons (Fig. 1B). Because the monkey's arm was continuously moving in a variety of directions, it was not possible to compute the preferred directions with respect to movement onset. Instead, preferred directions were computed at different time leads and lags relative to the measured firing rate of the neuron. For each time point t0 (discretized in 50 ms intervals) in the experiment, the firing rate of the neuron was measured (in 50 ms bins) and a movement trajectory was extracted at multiple time leads (blue line) and lags (red line) before and after t0 (Fig. 1B). At each time lead or lag, τ, of the extracted trajectory, the instantaneous movement direction of the hand was computed and discretized in 16 direction bins. The firing rate of the neuron measured at t0 was then associated with that direction bin. This was repeated over all time points in the experiment from which an average firing rate was computed for each direction bin. The preferred direction could then be estimated at that time lead or lag using a standard cosine tuning function (Georgopoulos et al., 1982).
Comparing the PD shifts in the random target-pursuit task (Fig. 1C) with those in the center-out task for the same neurons (using the optimal time window) (see Materials and Methods), it was evident that they closely matched each other. Over a population of 68 simultaneously recorded neurons, the mean correlation between the PD shifts in the two tasks was 0.88 (SD 0.23). By treating each PD as a velocity vector of unit length, preferred paths for each neuron could be computed by temporally integrating the PDs, which often revealed complex shapes (Fig. 1D). The fact that systematic shifts in PD were observed in two different tasks and that these shifts closely matched each other suggested that MI neurons may be more accurately described as encoding time-dependent trajectories.
The trajectory-encoding model
We sought to develop a tuning model that captured the encoding of complete trajectories instead of treating each temporal lead or lag independently, as we did in the last section. Given the binary nature of the occurrence of a spike, we applied a class of generalized linear models that estimated the probability of a spike emitted by a neuron (within a small spike sampling window) given that the monkey's hand generated a particular velocity trajectory. Consistent with previous work using a different analytic approach (Paninski et al., 2004b; Shoham et al., 2005), we found that an exponential relationship between the spiking probability and the inner product between the preferred velocity trajectory, k⇑, and the normalized velocity trajectory of the hand, v̂t0, fit the following data:
where each velocity trajectory extends over a range of times before and after the spike sampling time, t0. γ is an offset parameter of the model. The normalized velocity trajectory v̂t0 is a sequence of velocities (i.e., directions and relative speeds) over a predefined time duration (see Materials and Methods). We call k⇑ the preferred velocity trajectory because the inner product in the exponent is maximized when the v̂t0 vector is aligned with it. The absolute magnitude of each of the components of k⇑ are difficult to interpret in isolation. However, their magnitudes in relation to each other represent the pattern of variation in direction and relative speed over time that maximizes the probability of firing. Normalization of the velocity trajectory captures the shape of the velocity trajectory but ignores overall speed and position, which are known to influence the firing of motor cortical neurons (Georgopoulos et al., 1984; Moran and Schwartz, 1999). Therefore, the model in Equation 7 was extended by including terms reflecting the average speed and average position of the trajectory sample (see Materials and Methods).
Examples from the same four MI neurons as in Figure 1A–C demonstrate the fits of the exponential model to the empirical spiking probability of the neurons (Fig. 2A). For each neuron, the preferred velocity trajectories k⇑ exhibited time-dependent variations in both the x and y dimensions (Fig. 2B). By temporally integrating the x and y components of k⇑, the preferred paths or pathlets for each neuron were generated (Fig. 2C). The vector k⇑ is a set of dimensionless coefficients because it multiplies the normalized velocity trajectory, which is also unitless, in the exponent of Equation 7. However, it can be treated as if it had units of velocity because it is proportional to the normalized velocity trajectory that maximizes the probability of a spike occurring (as described above). Its time integral can, therefore, be interpreted as having units of position. The blue segment of the pathlet occurs before the spike sampling time, whereas the red segment follows the spike sampling time. For each neuron, the overall shape of the pathlet (Fig. 2C) resembled that of its preferred direction path to some degree (Fig. 1D). The differences in shape are partly attributable to the fact that the pathlet was created by temporally integrating relative velocity (i.e., direction and relative speed) components, whereas the preferred direction paths shown in Figure 1D integrate only direction information.
The trajectory-encoding model. A, The fitted model probability (red line) of a spike occurring in the spike sampling window as well as the empirical spiking probability (blue points) as a function of the linear predictor for the same four MI neurons as in Figure 1, A, C, and D. The linear predictor is the value of the exponent in the model. B, The x (solid line) and y (dashed line, above) components of the preferred velocity trajectories, k⇑, of the neurons in A. The y components of the preferred velocity trajectories are offset by 0.7 vertically for clarity. The solid and dashed horizontal lines indicate zero for x and y velocity parameters, respectively. Error bars represent SEs. C, The corresponding preferred paths or pathlets generated by temporally integrating the fitted k⇑. The blue portion of each preferred velocity trajectory and pathlet occurs before the spike sampling time whereas the red portion occurs after the spike sampling time. D, Scatter plots of the log tangential velocity versus the log radius of curvature of the four MI pathlets. The best fit lines (red lines) have slopes of 0.32, 0.24, 0.25, and 0.21, respectively. E, The corresponding preferred joint angle paths generated by the same encoding model expressing trajectories in a joint velocity coordinate system.
By examining the spatiotemporal properties of these estimated pathlets, we found that they shared a common feature that characterizes movement segments. Based on psychophysical work, it has been shown that complex movements can be decomposed into segments, each of which follows a power law by which the tangential velocity varies according to the one-third power of the radius of curvature (Viviani and Terzuolo, 1982; Lacquaniti et al., 1983; Viviani and Cenzato, 1985). Many of the pathlets closely followed the one-third power law (Fig. 2D). The mean power coefficient relating the tangential velocity to the radius of curvature among all pathlets (with ROC areas ≥0.6; n = 398 of 859 single units; for explanation of the ROC area, see Predicting spike occurrences from the model) was 0.29 (SD 0.06).
We also applied the same trajectory model in a joint velocity coordinate system by estimating the spiking probability of a neuron given shoulder and elbow angular velocity trajectories. The exponential fit the data as well as with Cartesian velocity trajectories (for a quantitative comparison, see Predicting spike occurrences from the model) and generated shoulder/elbow joint angle preferred paths that evolved in time and were consistent with the Cartesian pathlets (Fig. 2E). Because the data shown in Figure 2 are from a monkey using his left hand, a rightward-oriented pathlet in Cartesian coordinates (Fig. 2C, first panel) is associated with an increasing shoulder angle (Fig. 2E, first panel), whereas a leftward-oriented pathlet (Fig. 2C, second panel) is associated with a decreasing shoulder angle (Fig. 2E, second panel).
We examined the temporal stability of the pathlet shapes by calculating the pathlet of each neuron using either the first or second half of one experimental data set (n = 54 units) (Fig. 3A). The similarity of the pathlets across the two halves of the data were assessed by computing pairwise correlation coefficients between their corresponding preferred trajectories. The mean correlation coefficient was 0.75 (SD 0.21) compared with a mean correlation coefficient of −0.05 between randomly chosen preferred trajectories from different neurons.
Temporal and spatial stability of pathlets. A, Pathlets were generated using data from the first half (blue and red lines) and the second half (cyan and magenta lines) of one data set containing 54 neurons in one monkey. Only 15 neurons whose ROC areas are ≥0.6 are shown. B, Pathlets were generated using trajectories performed on the left side (blue and red lines) and right side (cyan and magenta lines) of the workspace. Only 12 neurons whose ROC areas are >0.6 are shown. To offset the reduced temporal and spatial data set sizes, sample movement segments were extracted every 25 ms for the temporal and spatial stability analyses.
We also examined the spatial invariance of the pathlet representations by partitioning the workspace into two halves and calculating the pathlet of each neuron over the left and right sides of the workspace (Fig. 3B). The mean correlation coefficient was 0.56 (SD 0.21) compared with a mean correlation coefficient of 0.06 between randomly chosen preferred trajectories from different neurons.
Predicting spike occurrences from the model
To assess the goodness of fit of the model to the data, we used the model (built on training data) to predict the occurrence of spikes given sample trajectory segments in test data. An ROC curve was generated that related the probability of correctly predicting a spike (hit probability) versus the probability of incorrectly predicting a spike when it was absent (false-positive probability) (Fig. 4A). The area of the ROC curve ranging from 0.5 to 1 quantifies the goodness of fit of the model. By systematically varying the trajectory length used to build the model in one data set (n = 54 units), we found that a 400-ms-long trajectory (−100 to +300 ms with respect to the spike sampling time) maximized the mean ROC area (Fig. 4B). Even when considering a short pathlet between +100 and +150 ms with respect to the spike sampling time (approaching a time-independent preferred direction representation), which is often treated as the optimal lag between a motor cortical neuron's firing and velocity (Moran and Schwartz, 1999; Paninski et al., 2004a), a 400 ms trajectory more accurately predicted the occurrence of a spike (or lack of a spike) under the model (p < 0.01, sign test) (Fig. 4C). Using a 400 ms trajectory, the mean ROC curve area was 0.60 (SD 0.12) with a maximum value of 0.96 over all 859 neuron samples recorded in MI during the RTP task. We also performed the same analysis in a joint-velocity coordinate system and found that a temporally extensive shoulder and elbow joint velocity trajectory of 400 ms (−100 to +300 ms with respect to the spike sampling time) maximized the mean ROC area over the same data set. There was no statistically significant difference in ROC area when comparing 400 ms trajectories in Cartesian-velocity versus joint-velocity coordinate systems (p > 0.05, sign test).
Predicting spike occurrences from the encoding model. A, ROC curves for the encoding models of 68 simultaneously recorded MI neurons in one data set (light gray lines). The four MI neurons shown in Figures 1 and 2 are shown in magenta, blue, green, and red, respectively. Each curve is computed from test data (10% of the data) not used to build the model (the remaining 90% of the data). B, Predictive power of the encoding model as a function of the trajectory length used for one data set. The predictive power is shown as one minus the area under the cross-validated ROC curve (averaged over all neurons in the data set) so that low values represent stronger predictive power. The length of the cyan and magenta lines represent the lead and lag times, respectively, for each trajectory length chosen. The bold blue and red lines represent the lead and lag times, respectively, of the 400 ms trajectory length (−100 to +300 ms) that generated the maximal predictive power. C, A scatter plot relating the areas under the ROC curves for a trajectory length of 400 ms (−100 to +300 ms with respect to the spike sampling time) versus a trajectory length of 50 ms (+100 to +150 ms) over all 54 MI neurons in one data set. Each data point represents the cross-validated ROC curve area for a given neuron. Notice that most points fall above the diagonal indicating stronger predictive power of the 400 ms trajectory compared with the short (+100 to +150 ms) trajectory. D, Direction decoding performance (mean absolute angular error) as a function of the duration of the set of lead and lag PDs used to predict instantaneous movement direction. A modified version of the population vector decoding algorithm was used such that, at each time point, each neuron voted multiple times using its preferred directions at times before (cyan) or after (magenta) the time point weighted by its firing rates at those lead and lag times. The bold blue and red lines represent the set of PDs with a duration of 400 ms (−100 to +300 ms).
To ensure that the estimated preferred trajectories were not an artifact of the autocorrelation of the hand movement, we also simulated a population of neurons whose preferred directions matched those of a population of 68 actual neurons recorded in one data set. This was accomplished by first determining the preferred directions of the actual neurons using a 100–150 ms delay between the firing of the actual neurons and the instantaneous direction of the monkey's hand. We then generated simulated Poisson spike trains by assuming that the firing probability of each of the simulated neurons was dependent on the similarity between the preferred direction and the instantaneous direction of the monkey's hand (supplemental information, available at www.jneurosci.org as supplemental material). We then applied the encoding model and estimated preferred trajectories for each of the simulated neurons. The mean correlation (over all 68 neurons) between the actual and simulated preferred trajectories was 0.26 (SD 0.42), which is low, suggesting that the actual pathlets were not solely caused by behavioral autocorrelations. We also used these simulated neurons to compare the predictive power of a simulated pathlet of 400 ms in length (i.e., −100 to +300 ms) with a simulated pathlet of 50 ms in length (i.e., +100 to +150 ms) (supplemental Fig. 2, available at www.jneurosci.org as supplemental material). The ROC areas were significantly larger using the 50 ms pathlet compared with the 400 ms pathlet (p < 0.05, paired t test). This should be contrasted with the same comparison made with the actual neurons (Fig. 4C).
Decoding
If time-dependent trajectory encoding more accurately characterizes the tuning of MI neurons, then it should be possible to decode movement more effectively using time-dependent representations. To directly compare our decoding results with those using the prevalent population vector decoding algorithm using a time-independent PD for each neuron, we characterized the tuning of each neuron with a set of temporally shifting PDs (as in Fig. 1C) instead of using the trajectory-encoding model described in Equation 7. Using one data set (n = 54 units) recorded using the RTP task (animal RJ), we set aside 90% of the data to find the best fit set of PDs for each neuron (i.e., the training data) and then decoded the instantaneous hand direction using the remaining 10% of the data (i.e., the test data). For each 50 ms time bin in the test data, we predicted the instantaneous hand direction by computing a population vector based on multiple PDs for each neuron in the population. That is, each neuron voted multiple times using its PDs at multiple time leads and lags weighted by its firing rates at those leads and lags. We then compared the predicted and actual hand direction for each time bin in the test data and computed the average absolute error. As with the trajectory-encoding model, we systematically varied the number of lead and lag PDs used to decode the hand direction (Fig. 4D). The average absolute direction error decreased with increasing PD set size but eventually leveled off using a set of PDs spanning a period of ∼400 ms (100 ms before to 300 ms after the measured firing rate). This is consistent with our findings using the trajectory-encoding model where we found that a 400-ms-long trajectory (100 ms before to 300 ms after the spike sampling time) maximized the ROC curve area. The direction error was over 10° smaller using a 400 ms PD set compared with a standard population vector approach using a single PD measured between 100 and 150 ms after the firing rate. To verify that the superior decoding performance using the extended PD set was not attributable to the fact that each neuron had its own optimal lag, we applied the standard population vector approach in which each neuron voted once based on its PD measured at the neuron's optimal lag. We first found the optimal time lag for each neuron separately by applying the decoding algorithm using only that neuron and found the lag time that minimized the average absolute direction error. We then applied the standard population vector approach using all the neurons with their respective optimal lags and corresponding PDs. Even in this case, the mean absolute error was 45° (SE 0.71°), which was larger than the error using the extended PD set.
Topographic organization of pathlet representations
We examined evidence of any spatial structure among the pathlet representations of the neurons recorded within and across electrode sites on the arrays implanted on the precentral gyrus (Fig. 5A). We compared the shapes of pathlet pairs estimated from neurons (all recorded units with ROC areas ≥0.6) at different interelectrode distances by correlating their preferred velocity trajectory vectors, k⇑. As others have found when examining directional tuning (Amirikian and Georgopoulos, 2003; Ben-Shaul et al., 2003), we found a significantly higher correlation between neurons recorded on the same electrode or neighboring electrodes (mean correlation of 0.08, SD 0.47) compared with neurons separated by 0.5–1 mm (mean correlation of 0.03, SD 0.45; p < 0.05, two-sample t test). The mean correlation between neurons recorded on the same electrode was 0.14 (SD 0.46). Aside from this finding, no simple topographic structure was evident. A map of pathlets based on a population of simultaneously recorded single units in MI demonstrated a variety of different shapes and orientations over the surface of the cortex (Fig. 5B).
Multielectrode array placement and the spatial distribution of pathlets. A, The placement of the electrode arrays in the MI of the three monkeys based on digital photographs taken during the surgical implantation. CS, Central sulcus; Arc, arcuate sulcus; pcd, precentral dimple; spur, spur of the arcuate sulcus. B, The map of pathlets computed over a set of simultaneously recorded MI neurons in one data set. The relative position of each pathlet is determined by the relative position of the electrode on the electrode array that recorded the pathlet of the neuron. Although we recorded from 96 of the 100 electrodes on the array, not all electrodes recorded single units. Moreover, only pathlets whose ROC areas were ≥0.6 are plotted. On some electrodes, more than one unit was isolated. However, for purposes of clarity, the pathlet from only one unit (arbitrarily chosen) is plotted at each electrode site.
Torque trajectory encoding
Because complex movements characterized kinematically depend on the forces or torques applied to the joints, we applied the same encoding model to estimate the spiking probability of a neuron given a temporally extended torque trajectory. Shoulder and elbow torques were estimated using the inverse dynamics equations that consider the inertial properties of the monkey's arm and the robotic device (Scott, 1999). We again found that the exponential model fit the data well (Fig. 6A), yielding preferred torque trajectories (Fig. 6B, shoulder torque, solid lines; elbow torques, dashed lines). By varying the torque trajectory length used in the encoding model, we again found that a 400 ms (−100 to +300 ms with respect to the spike sampling time) pathlet optimally predicted the occurrence of a spike. The mean ROC area over a simultaneously recorded population of 54 units was 0.55 (SD 0.13). The preferred torque trajectories revealed temporal variations such as monophasic fluctuations (Fig. 6B, two top panels) or multiphasic fluctuations often associated with more complex pathlets (Fig. 6B, bottom).
Torque trajectory encoding. A, The fitted torque model probability (red line) of a spike occurring in the spike sampling window as well as the empirical spiking probability (blue points) as a function of the linear predictor for three neurons computed using the same encoding model as with velocity, but using 400-ms-long torque trajectories (−100 to +300 ms with respect to the spike sampling time) instead. B, The associated preferred torque trajectories (shoulder torque, solid lines; elbow torque, dashed lines). The elbow components of the preferred velocity trajectories are offset by 0.1 vertically for clarity. The solid and dashed horizontal lines indicate zero for shoulder and elbow torque parameters, respectively. Error bars represent SEs. C, The associated pathlets for the same three neurons using the velocity trajectory-encoding model.
Velocity versus torque trajectory encoding
It is difficult to compare kinematic versus kinetic encoding under a fixed dynamical context because there is a one-to-one mapping between the kinematics of the hand and the torques generated at the joints. Therefore, we performed an additional experiment with one animal in which a viscous load was applied to both joints on certain trials randomly intermingled with unloaded trials using the same RTP task. We then applied the velocity and torque trajectory-encoding model to loaded and unloaded trials separately. We assumed that a neuron that generated a preferred velocity trajectory that remained relatively unchanged under both unloaded and loaded conditions (compared with its preferred torque trajectory) could be considered more of a “kinematic” encoding neuron. Likewise, a cell that generated a preferred torque trajectory that remained relatively unchanged (compared with its preferred velocity trajectory) under both conditions could be considered more of a “kinetic” encoding neuron. Although we observed both kinematic as well as kinetic encoding neurons, there was a strong bias toward kinematic encoding. Six example neurons exhibit relatively invariant preferred velocity trajectories (Fig. 7A) and corresponding pathlets (Fig. 7B) whereas their preferred torque trajectories were altered more substantially under the two loading conditions (Fig. 7C). Over a population of neurons (two data sets in one animal; 221 single unit samples), there were considerably more neurons whose preferred velocity trajectories were more strongly correlated under the two conditions compared with their preferred torque trajectories (p < 0.001, paired t test on Fisher z-transformed correlations) (Fig. 7D).
Comparison of velocity and torque trajectory encoding in the loading RTP task. A, The x (solid lines) and y (dashed lines) components of the preferred velocity trajectories in unloaded trials (blue and red) and loaded trials (cyan and magenta) for six MI neurons. The y components of the preferred velocity trajectories are offset by 0.5 vertically for clarity. B, The associated pathlets in unloaded (blue and red) and loaded (cyan and magenta) trials. C, The shoulder (solid lines) and elbow (dashed lines) preferred torque trajectories for the same six neurons as in A in unloaded (blue and red) and loaded (cyan and magenta) trials. The elbow components of the preferred torque trajectories are offset by 0.5 vertically for clarity. D, A population scatter plot relating the correlation coefficients between the preferred velocity trajectories in the unloaded and loaded trials and the correlation coefficients between preferred torque trajectories in the two conditions. Each point represents one MI neuron among 70 neurons with ROC areas ≥0.6 recorded over a total of 221 single unit samples in two data sets in one monkey.
Discussion
We have provided three important pieces of evidence to support the hypothesis that single motor cortical neurons encode time-dependent trajectories. First, we have demonstrated that a temporally extended preferred trajectory (expressed in either Cartesian or joint-velocity coordinates) can more accurately predict the occurrence of a spike (or lack of a spike) than a shorter trajectory that approaches a time-independent preferred direction. The improved predictive power of the model using a temporally extended trajectory compared with shorter pathlet lengths is not because a larger number of parameters was used in the model. Because we reduce the dimensionality of the trajectory using principal component analysis, all models have 14 free parameters regardless of the trajectory length. Moreover, we applied the models to predict spike occurrences on test data that were not used to build the model.
Second, we have shown that temporally extensive trajectory encoding can be used to decode the instantaneous movement direction more accurately. A 400-ms-long set of preferred directions reduced the mean angular direction error by over 10° compared with a shorter set of PDs spanning the typically assumed delay between neural activity and movement (i.e., +100 to +150 ms).
Last, trajectory encoding may explain the temporal shifts in preferred direction that are often observed in a standard center-out task when measured on a finer time scale (Mason et al., 1998; Sergio and Kalaska, 1998; Sergio et al., 2005). This is a phenomenon that has been underappreciated but calls into question the idea that motor cortical neurons possess a fixed preferred movement direction on a fine time scale. The shifts in preferred direction during reaching may be explained in part by the fact that force direction shifts from the early, ballistic portion of a reach movement to the late, decelerating phase when braking occurs (Sergio et al., 2005). Shifts in preferred direction would thus occur if neurons were encoding force direction. However, this does not explain the shifts in preferred direction that occur before movement in an instructed delay task, as Mason et al. (1998) have shown, and as is evident in Figure 1A.
The hypothesis that single motor cortical neurons encode trajectories has been examined previously (Hocherman and Wise, 1991). In that study, the authors observed that motor cortical neurons fire preferentially for certain types of hand paths even when the initial and final positions of the paths were identical. In the present study, we extended these initial observations by developing an explicit trajectory-encoding model or tuning function that characterized the shape of the preferred trajectories.
The notion that the spiking probability of a motor cortical neuron would be related to the velocity of the hand 300 ms in the future might seem somewhat counter-intuitive. Spike-triggered EMGs of corticomotoneuronal (CM) cells as well as electrical stimulus-triggered EMGs indicate a short-lasting effect with a short latency of ∼10 ms (Fetz and Cheney, 1980; Park et al., 2004). However, during motor behavior the latency between MI modulation and EMG onset or movement onset can be quite variable and can extend to several hundred milliseconds (Evarts, 1974; Porter and Lewis, 1975; Thach, 1978; Holdefer and Miller, 2002). Even CM cells exhibit a broad range of latencies that can extend to several hundred milliseconds between their firing onset and EMG onset of the muscle they facilitate during movement (Cheney and Fetz, 1980). Therefore, the long temporal relationship between the spiking probability of a motor cortical neuron and arm motion is not surprising in the context of actual motor behavior.
The experimental evidence based on our loading experiments suggests that the fundamental representations of motion encoded by the motor cortex may be better characterized as complex kinematic trajectories instead of complex kinetic trajectories. This conclusion appears inconsistent with the results of Sergio and colleagues (Sergio and Kalaska, 1998; Sergio et al., 2005), who also demonstrated systematic shifts in preferred direction in a reaching task. They argued that the shifts in preferred direction were a consequence of variations in force production by showing that the decoded output of the population vector often, but not always, corresponded more closely to the dynamics of hand force. One explanation for this discrepancy may be the fact that we recorded rostrally on the precentral gyrus whereas Sergio and colleagues (Sergio and Kalaska, 1998; Sergio et al., 2005) recorded caudally from the anterior bank of the central sulcus. There is evidence indicating that neurons in the rostral portion of the motor cortex may be more related to kinematic variables such as velocity and movement direction (Ashe and Georgopoulos, 1994; Crammond and Kalaska, 1996). Therefore, our conclusion regarding a kinematic trajectory representation may apply primarily to the rostral portion of motor cortex. However, regardless of whether a kinematic or kinetic representation is more appropriate, we demonstrated that complex, time-dependent trajectory encoding (in a kinematic or kinetic domain) captures the firing properties of single motor cortical neurons more accurately.
Natural movement trajectories reside in a relatively small dimensional space because of physical and biological constraints such as the one-third power law. This is supported by our principal component analysis, which demonstrated that a high-dimensional trajectory can be well characterized with 10 dimensions. However, in contrast to an efficient algorithmic description in terms of coefficients with respect to an orthogonal basis, it appears that motor cortical neurons represent an extended set of motion templates. The preferred trajectories of millions of motor cortical neurons that are likely to be modulating could easily form such a template set for trajectory representation.
Future work will explore how individual preferred trajectories are combined when multiple neurons fire together. The exponential nature of the encoding model provides a very simple and elegant prediction by which individual preferred trajectories are combined when more than one neuron fires simultaneously. If the neurons are firing independently (conditioned on a given trajectory performed by the animal), the probability that multiple neurons fire simultaneously will equal the product of their individual firing probabilities and, therefore, their preferred trajectories will add.
Trajectory encoding has a conceptual appeal in that it represents a unified framework for understanding why motor cortical neurons appear to encode so many different kinematic features including direction, velocity, speed, and position. As we have characterized it, a trajectory embodies a time-dependent velocity profile that implicitly carries both direction and relative velocity information. At the same time, our model incorporates global position and speed dependencies, implying that neurons prefer trajectories in a certain location moving at an overall speed. Instead of encoding any single abstract Newtonian movement parameter, the motor cortex represents movement fragments that embody many parameters extended in time.
Footnotes
-
This work was supported by a grant from the Whitehall Foundation and National Institute of Neurological Disorders and Stroke Grants N01-NS-2-2345 and R01 NS45853-01A2. We thank Greg Ojakangas for deriving and implementing the equations of motion of the exoskeletal robot together with the animal's arm. We also appreciate Wei Wu's assistance in certain computational analyses. We also thank Matthew Fellows, Elise Gunderson, Zach Haga, and Dawn Paulsen for their help with the surgical implantation of the arrays, training of monkeys, and data collection.
- Correspondence should be addressed to Nicholas Hatsopoulos, Department of Organismal Biology and Anatomy, University of Chicago, Chicago, IL 60637. nicho{at}uchicago.edu
