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

This Article Abstract Full Text (PDF) All Versions of this Article: 96/5/2578    most recent 01086.2005v1 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 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 (5) Citing Articles via Google Scholar Google Scholar Articles by Townsend, B. R. Articles by Lemon, R. N. Search for Related Content PubMed PubMed Citation Articles by Townsend, B. R. Articles by Lemon, R. N.

Linear Encoding of Muscle Activity in Primary Motor Cortex and Cerebellum

¹Sobell Department of Motor Neuroscience and Movement Disorders, Institute of Neurology, and ²Gatsby Computational Neuroscience Unit, University College London, London, United Kingdom; and ³Department of Statistics, Columbia University, New York, New York

Submitted 17 October 2005; accepted in final form 18 June 2006

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The activity of neurons in primary motor cortex (M1) and cerebellum is known to correlate with extrinsic movement parameters, including hand position and velocity. Relatively few studies have addressed the encoding of intrinsic parameters, such as muscle activity. Here we applied a generalized regression analysis to describe the relationship of neurons in M1 and cerebellar dentate nucleus to electromyographic (EMG) activity from hand and forearm muscles, during performance of precision grip by macaque monkeys. We showed that cells in both M1 and dentate encode muscle activity in a linear fashion, and that EMG signals provide predictions of neural discharge that are equally accurate to those from kinematic information under these task conditions. Neural activity in M1 was significantly more correlated with both EMG and kinematic signals than was activity in dentate nucleus. Furthermore, the analysis enabled us to look at the temporal properties of muscle encoding. Cells were broadly tuned to muscle activity as a function of the lag between spiking and EMG and there was considerable heterogeneity in the optimal delay among individual neurons. However, a single lag (40 ms) was generally sufficient to provide good fits. Finally, incorporating spike history effects in our model offered no advantage in predicting novel spike trains, reinforcing the simple nature of the muscle encoding that we observed here.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Since the pioneering work of Evarts (1968) many studies have investigated the problem of the coding of motor command signals in primary motor cortex (M1), by relating the responses of neurons recorded in this area to the modulation of different extrinsic movement parameters. A key finding was that M1 neurons are tuned to the direction of movement during a reaching task (Georgopoulos et al. 1986). This has been extended to incorporate representations of kinematic parameters, such as position, velocity, and acceleration (Fu et al. 1995; Moran and Schwartz 1999; Paninski et al. 2004a,c; Scott and Kalaska 1997). Tuning for static and dynamic force has also been extensively demonstrated, in wrist movement and precision grip tasks (Cheney and Fetz 1980; Hepp-Reymond et al. 1999; Porter and Lemon 1993).

Comparatively few studies have addressed the encoding of intrinsic movement parameters such as patterns of muscle activity (Holdefer and Miller 2002; Kakei et al. 1999; Morrow and Miller 2003). This is despite much evidence of functional modulation of connectivity between M1 corticomotoneuronal (CM) cells and muscles used during task performance (Bennett and Lemon 1996; Fetz and Cheney 1980; Jackson et al. 2003). Furthermore, theoretical studies suggest that "muscle space" may be the fundamental coordinate system of M1, given that apparent encoding of kinematic parameters might arise from position- and velocity-dependent compensations made during the direct control of muscle activation (Mussa-Ivaldi 1988; Todorov 2000). It remains unclear which, if any, of these variables are explicitly encoded by M1.

The contribution made by the cerebellum toward movement control has been studied using similar methods. Just as for M1, tuning of cerebellar neurons, both in the cortex and the deep nuclei, has been observed for several parameters including direction of reaching (Fortier et al. 1989), position and/or velocity (Fu et al. 1997; Mano and Yamamoto 1980), force (Smith and Bourbonnais 1981), and muscle activity (Wetts et al. 1985) across a variety of tasks.

Thus it would be of interest to determine the relative importance of different parameters in the encoding of movement and to compare this encoding in M1 and cerebellum. To try to tackle these issues, we applied a basic linear–nonlinear (LN) analysis (Chichilnisky 2001; Simoncelli et al. 2004) to describe the encoding of muscle activity and kinematic parameters by single neurons in M1 and dentate nucleus, in the awake behaving macaque monkey performing a precision grip task. This type of analysis is ideal for relating neural discharge to the complex time course of the electromyographic (EMG) signal and has already been applied to describe the encoding of hand position and velocity by neurons in macaque M1 (Paninski et al. 2004c). We demonstrate that the activity of single neurons in both regions can be accurately predicted using simultaneously recorded EMGs from up to nine hand and forearm muscles. This tuning of M1 and dentate neurons to muscle activity is linear in nature and contrasts with a more nonlinear tuning for kinematic information, whereas the relation of neuronal discharge to both types of parameters is on average weaker in the dentate nucleus compared with M1, in agreement with previous studies of the cerebellar nuclei and cortex (Fu et al. 1997; Thach 1978; van Kan et al. 1993). Overall, our results are consistent with a system that controls muscle activity in a linear fashion (Ethier et al. 2006; Todorov 2000), and thus linear methods are sufficient to describe this encoding, in M1 and cerebellar dentate nucleus.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Experimental procedures

Behavioral task. Three purpose-bred adult monkeys (M. mulatta) were used for this study: two females (M36 and M38) and one male (M41), weighing 5.0, 6.0, and 6.0 kg, respectively. The animals performed a variant of the precision grip task (see Baker et al. 2001), which involved squeezing two levers between thumb and index finger. Levers were mounted onto the spindles of motors (Phantom, SensAble Technologies), computer controlled to produce a variety of position-dependent forces. One complete trial involved moving both levers into a target displacement window and maintaining this position for 1 s before releasing (Fig. 1A). For monkey M36, this target window was between 6 and 8 mm from baseline for both levers. Monkeys M38 and M41 were not as accurate in their performance of the task, so slightly larger windows between 5 and 8 mm were used. Three auditory cues were presented to the animal: the first indicated that both digits were in the target displacement window, the second that they had been held there for the required duration, and the third was accompanied with a fruit reward once the levers had been released and returned to baseline.

View larger version (19K): [in this window] [in a new window]   FIG. 1. Example trial performed by monkey M38 selected according to the off-line criteria described in METHODS. A: thumb (thin line) and finger (thick line) lever position traces. Dashed lines indicate target position window for levers. B: lever velocity profiles. Dashed line indicates 30 mm/s velocity threshold from which movement period is defined.

Off-line, two time periods were defined (Fig. 1A). The movement period was defined as the time during which either finger or thumb velocity was >30 mm/s. The hold period began once both finger and thumb positions were within the target window and lasted for 1 s. These periods were used for off-line trial selection criteria applied to data sets that were used for the comparison of fits by the model to EMG and kinematic data. Here it was important to ensure that trial performance was as homogeneous as possible because greater variability in task performance could weaken the correlation of kinematic parameters with neural discharge. In all animals, roughly 40% of trials performed passed these criteria and were available for analysis.

The criteria used were as follows: First, the movement period (defined by finger or thumb velocity being >30 mm/s) had to last <1 s (Fig. 1B). This criterion rejected trials where the finger or thumb levers did not move swiftly into the target displacement window on the first attempt. Second, the length of time that either finger or thumb levers were kept in their hold windows was >0.5 s but <1.25 s. This criterion rejected trials where the hold was not achieved directly and where the levers did not rapidly return to baseline at the end of the hold period. This latter point was important because, although successful completion of the 1-s hold was signaled by an auditory signal (see above text), sometimes the monkey kept the levers within the target window for a short period after instead of immediately terminating the trial. An example of an accepted trial is illustrated in Fig. 1A.

Recording. Details of surgical procedures and the Eckhorn multiple-electrode recording system (Thomas Recording, Marburg, Germany) were previously described (Baker et al. 1999, 2001) for recordings from primary motor cortex (M1). All procedures were performed in accordance with appropriate UK Home Office regulations. Data were recorded directly to hard disk by two A2D cards (PCI-6071E, National Instruments).

EMG recordings.    Monkeys M36 and M38 were implanted with subcutaneous EMG patch electrodes (Miller et al. 1993) sutured directly onto the surface of intrinsic hand and forearm muscles, in the hand used to perform the task. The electrode leads ran subcutaneously to a connector on the monkey’s back. The muscles implanted and their abbreviations in this paper are detailed in Table 1. EMGs were recorded bipolarly with gains of 1,000–5,000, high-pass filtered at 30 Hz (NL824, Digitimer), and were sampled at 5,000 Hz. This was downsampled to 500 Hz for purposes of the analysis described below.

The intrinsic hand muscles are particularly important for the control of skilled finger movements, including the generation of finely graded force during precision grip in humans and monkeys (Maier and Hepp-Reymond 1995; Porter and Lemon 1993). We therefore recorded EMGs from the intrinsic hand muscles 1DI and two muscles of the thenar eminence (flexor pollicis brevis and adductor pollicis brevis) that, because of their close proximity and the small size of the macaque thumb, were sampled by a single implanted electrode and are referred to as the "thenar" EMG. The long flexors of the fingers (FDP and FDS) were recorded from because they too are primarily involved in the generation of pinch force by the index finger. The extensors of the wrist (ECU, ECR-L), extensors of the fingers (EDC), and flexors of the wrist (FCU) were implanted because these would be active during the removal and insertion of the hand and digits from/into the manipulandum. These muscles also act to stabilize joint torques and maintain equilibrium during precision grip (Maier and Hepp-Reymond 1995; Schieber and Santanello 2004). Coordinated activity in all these muscles is required to perform the precision grip and thus we recorded EMGs from multiple muscles simultaneously rather than one at a time.

To assess the extent of correlation between EMG signals we calculated the mean absolute correlation coefficient between each possible pairwise combination of EMGs, across all data sets in monkey M38 (nine EMGs = 36 possible pairs) and M36 (seven EMGs = 21 possible pairs). The mean level of correlation was 0.31 (SD 0.10) for M38 and mean 0.54 (SD 0.03) for M36, indicating quite low background levels of correlated activity. However, in both animals particular pairs of EMGs showed strong correlations, such as ECU x 1DI in M38 (0.67) and thenar x 1DI in M36 (0.85), consistent with the similar actions that these muscles exert on the hand. In contrast, other pairs showed particularly weak correlations, such as AbPL x FDP in M38 (0.05) and EDC x FDS in M36 (0.3), in agreement with the different actions that these muscles exert on the hand. Overall there was a mixture of correlation strengths between pairs of muscles, with no clear tendency toward very strong or weak correlations across all pairs. Note that the level of physical cross talk in these EMG recordings was very low (see Brochier et al. 2004).

Task data.    Digital events (trial start and end times, end of hold period) were recorded, together with lever position signals sampled at 500 Hz.

Cortical recordings.    Recordings were made in the hand area of M1, contralateral to the performing hand. Right M1 was recorded in animal M36 and left M1 in M38 and M41. M1 chamber center coordinates were about A10 L17 in all three animals. At least five glass-insulated platinum electrodes (impedance 1–3 M, 4 x 4 grid with interelectrode spacing of 300 µm) were independently lowered into the cortex to search for cells. Pyramidal tract neurons (PTNs) were identified by their antidromic response to stimulation in the pyramid (latencies, 0.9–4 ms; thresholds, 20–200 µA) and collision testing (Lemon 1984). All other cells in M1 that did not respond to stimulation were labeled as unidentified neurons (UIDs). In M38 and M41 at least two electrodes were inserted into the cerebellum in each session, ipsilateral to the performing hand, to make recordings from dentate nucleus simultaneously with M1. Electrodes were introduced by fine sharpened guide tubes that were advanced through the cortical dura, 5 mm below it, before the electrode was advanced. The total depth of penetration was about 25–30 mm. Cerebellar chamber center coordinates were A8 and L6 in M38 and A9.5 L6.5 in M41. During a session, we selected cells for recording that showed clear task-related modulation in their firing rate and whose interspike-interval histograms did not contain counts in bins at short intervals (<2 ms), which is evidence that the recorded spikes we discriminated came from a single neuron. Neuronal activity was recorded as the analog activity, filtered between 1 and 10 kHz, and sampled at 25 kHz. We were typically able to make stable recordings from single neurons for around 30 min. Off-line, single units were discriminated using principal component analysis on the spike waveform and cluster cutting (Eggermont 1990).

Analysis

Instantaneous firing rate estimation. For each discriminated unit, an estimate of the instantaneous firing rate (IFR) throughout the recording period was first calculated, according to techniques described in Pauluis and Baker (2000). Briefly, this method uses the reciprocal of the interspike interval as a first approximation to the IFR, explicitly detecting significant changes in firing rate while smoothing the periods in between these times. For this analysis, IFR profiles were calculated with a sampling resolution of 500 Hz and smoothed with a Gaussian kernel of width 10 ms.

Encoding model. To study how cells in M1 and dentate encode muscle activity, we described the dependency of cell firing rate on EMG activity from multiple muscles using a linear–nonlinear (LN) model (Fig. 2). This analysis used standard techniques based on spike-triggered regression. Full descriptions of these procedures can be found in previous work (Chichilnisky 2001; Paninski et al. 2004c; Shoham et al. 2005; Simoncelli et al. 2004). Applications of these methods to EMG data are outlined as follows.

View larger version (24K): [in this window] [in a new window]   FIG. 2. Predicting primary motor cortex (M1) and dentate instantaneous firing rate (IFR) from muscle activity. A: sample-rectified electromyographic (EMG) signals from 5 hand and forearm muscles. B: muscle activity is linearly filtered by the muscle weight vector to give the filtered signal · . C: · is transformed by the nonlinearity f to produce a nonlinear prediction of the IFR for that neuron. Observed IFR shown for comparison. Figure displays a randomly chosen section of the data. Analysis is for a single M1 pyramidal tract neuron (PTN).

Regression was used to estimate the cell’s weighting of the EMG activity from each muscle , as ). The first term is the inverse of the correlation matrix of w, which is computed to remove any correlations of the EMG input vector with itself. This is a key requirement of spike-triggered regression techniques (the probability distribution for values of the input vector is assumed to be radially symmetric or "white"; Chichilnisky 2001). The second term is the cross-correlation of the EMG input vector with the spike train, which reduces to a conditional expectation (E) because of the binary nature of the spike train.

The term , as the spike-triggered average (STA), gives the average value of the normalized EMG activity in each muscle at = 40 ms after the spike. can be conceptualized as describing how the cell weights activity in each of the muscles at this lag (Fig. 2B, top).

The relationship between spiking of the cell and EMG activity in the set of muscles is then captured by the term ( · ), which is a linearly filtered version of the concatenated muscle activity in carried out by the cell using its weight vector (Fig. 2B, bottom). To determine the value of the filtered signal at time t, the EMG level in each muscle at lag after time t is weighted according to entries of the fixed vector . These values are then summed to give the filtered signal. Thus a muscle with a positive weight will increase the amplitude of the filtered signal at time t + ; a muscle with zero weight will make no contribution; and a muscle with a negative weight will decrease the amplitude of the filtered signal at time t + .

Cells in M1 can show nonlinearity in the relation of their discharge to movement parameters such as hand position and velocity (Paninski et al. 2004c). Therefore we included a nonlinear term f in the description of muscle encoding given above. In the LN model, the filtered signal ( · ) controls cell firing rate through this nonlinearity f, which is written f( · ). We did not assume a particular type of nonlinearity beforehand, but instead estimated f separately for each cell from ( · ) using an intuitive, nonparametric binning process. This procedure consists of finding, for any possible value u of the filtered signal · , all times {t}_u at which · was found to be approximately equal to u. The conditional firing rate f(u) is then given by the fraction of time bins {t}_u that contained a spike. In this procedure, f(u) has an underdetermined scale factor [because a scale factor in the argument u-axis can always be absorbed by rescaling f itself (Chichilnisky 2001)]. Therefore we standardized u by linearly mapping the 1st and 99th quantiles of the observed distributions of u to –1 and +1, respectively, in all plots.

In Fig. 2C it can be seen that f( · ) gives the estimated firing rate of the cell for each value of ( · ). Mapping the signal ( · ) through the nonlinearity f therefore provided an estimate of the instantaneous firing rate of the cell during the experiment. This relationship is encapsulated in Eq. 1, including terms for the current time bin of width dt (here this was 2 ms, small enough so that only one spike was observed per bin) centered at time t

(1)

Cross-validation. We used the following cross-validation procedure to test the accuracy of the cascade model at predicting cell IFR from EMG activity. The analysis was performed on sections of data from successfully performed trials only, and periods of inactivity between the end of each trial and the start of the next one were omitted. For the data set from each recording session, 60% of the trials were designated training data and 40% of trials as test data. We selected trials for each of the two sets in an interleaved fashion, running through the total period of time analyzed from a given recording session (i.e., trials 1, 3, and 5 = train; trials 2 and 4 = test) to ensure that no order effects were introduced. These data sets were separate and nonoverlapping: each trial was allocated to only one set.

For each neuron, the training set was used to fit and the nonlinearity f( · ). First, the signal · was computed using from the train set and from the test set. As described above, this signal captures the relationship between cell discharge and muscle activity. In effect, it is a linear prediction of cell discharge made from the EMG (although scaled between –1 and +1 as detailed above). This prediction was then compared with the cell’s observed IFR for the test set simply by computing the correlation coefficient of the two signals. This provided a measure of how accurate the linear stage of the model was at predicting cell activity from a nonoverlapping data segment.

After this, f( · ) was used to generate the nonlinear prediction of cell IFR using f and fitted from the training set, and from the test set. Again, this nonlinear prediction was compared with the observed IFR for the test set by means of the correlation coefficient, to give a measure of how accurate the nonlinear stage of the model was at predicting cell activity. In addition, all example encoding functions shown in RESULTS and Fig. 2 were also cross-validated by fitting to the training subset and f( · ) to the test subset.

This cross-validation procedure provided a measure of how accurate our analyses were in predicting the cell’s activity rather than simply reproducing observed data. In the analysis, steps were taken to minimize overfitting, whereby predictions of the test data decrease in accuracy when too many regressors are fitted to the training data. These steps are detailed where relevant.

We also carried out the same analysis on kinematic instead of EMG data, forming the input signal from the finger position and velocity signals, at the same lags that were used for the EMG.

Although the analysis studied the combined activity of several concurrent EMGs, it was possible that a small subset of these muscles could dominate the predictions of cell discharge made from the total combined EMG signal. To address this, we first calculated the frequency at which each muscle was allocated the strongest absolute weight by . In M38 all muscles were not weighted equally: instead, two muscles in particular tended to be allocated the strongest weights (EDC and FDP). This may in part be attributed to the strong task-dependent activity of these muscles: EDC showed strong activity during the release phase at the end of each trial because it acted to extend the digits as they were removed from the manipulandum, whereas FDP would be critically involved in generating and maintaining force during the grip (Maier and Hepp-Reymond 1995). In M36, the thenar muscles and FDP tended to show the strongest weighting across the population. Thus some muscles were weighted more strongly than others across the population.

Second, for each cell, nonlinear predictions of IFR were generated separately for each of the nine EMGs from M38 and each of the seven EMGs in M36. These were compared with the observed IFR by calculating the correlation coefficient between the two signals. We then found the frequency at which each muscle provided the best prediction of cell IFR across the population. In M38, thenar, 1DI, and ECU tended to give the best predictions when EMGs were fitted one at a time; in M36, it was 1DI and EDC. Thus for a given cell, it was not possible to determine in advance the best set of EMGs for making the IFR prediction from looking at the pattern of weighting in . For example, a cell that allocated a strong regression weight to muscle 1DI might lead to a poor prediction when cell discharge was fitted to EMG from that muscle alone because the cell’s activity was dependent on the difference in activity between 1DI and other muscles with different weights.

Spike history. Our model was used to predict the time-varying firing rate of the cell, given by f( · ). In addition, we tested the accuracy of our model at explicitly generating spike trains based on fits to EMG data. Spike trains were generated by an iterative method: in a separate test portion of data, the average spike count in each time bin t was calculated according to a Poisson process, with the rate determined by f( · ) (Chichilnisky 2001). Because this method of predicting spikes bin by bin is Poisson in nature, it is inherently noisy, so predictions were repeated 20 times for the same section of data, summed, and averaged. This average firing rate, derived from spike trains generated by the model, was compared with the observed IFR for the same period by means of the correlation coefficient to give a measure of the accuracy of our model in predicting the spiking behavior of each cell.

However, neural responses depend on the spiking history of the cell (Berry and Meister 1998; Keat et al. 2001; Paninski et al. 2004b; Truccolo et al. 2005). Therefore we fitted M1 and dentate neurons to an adjusted LN model, which incorporated some of this response history. This was implemented in two steps. First, we formed a modified input vector _H, formed by concatenating the original input vector and the spike count of the neuron in the previous five bins

(2)

where w₁ is the first element of the original and w_end is the last and r_{i} denotes the observed spike count i time bins ago. The linear filter fitted to this modified input vector is referred to as _H. A single lag was chosen here in contrast to the full filter length (see Multiple filter delays section below) to minimize overfitting.

Second, in a separate test set of data, a spike train was again generated iteratively, using _H and _H

(3)

Here each element r_{i} of _H was the model’s prediction of the neuronal spike count i time bins ago, instead of the observed spike count used in Eq. 2. Simulating response history in this way, by "feeding back" the model’s output rather than using the cell’s real spike history, ensured that we did not contaminate our model’s predictions with spike trains that had already been observed. Again, predictions of the spike count were repeated 20 times, averaged, and compared with the IFR, which enabled us to compare the prediction accuracies of the spike history and basic LN models.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
Data set

Some 51 data sets from three monkeys were analyzed, constituting a total of 216 neurons. We divided these data as follows. For fits of the encoding model to both EMG and kinematic signals, data from M36 and M38 were used. Cells recorded from these animals were split into three groups: M1 pyramidal tract neurons (PTNs), unidentified M1 neurons (UIDs), and cells from cerebellar dentate nucleus. This gave 22 PTNs and 14 UIDs from monkey M36, together with 29 PTNs, 32 UIDs, and 33 dentate cells from monkey M38. Sixteen of 22 PTNs in M36 exhibited significant postspike facilitation or suppression effects in EMG of one or more hand muscles, which is evidence that they were corticomotoneuronal (CM) cells. Additional data for the analysis of fits to kinematic information came from recordings in monkey M41 (no EMGs were recorded in this animal). This constituted 63 M1 neurons and 23 dentate neurons.

We first describe linearity of muscle activity encoding by these neurons and the temporal properties of this encoding. Next, we compare this encoding to tuning for kinematic parameters. Finally, the results of an attempt to model precise spiking behavior in a subset of this population are discussed.

Nonlinear encoding

We computed the EMG encoding functions f( · ) for a total of 130 cells from M1 and dentate nucleus. The majority were nonlinear in nature; therefore it was of interest to determine whether predictions of each cell’s instantaneous firing rate were more or less accurate when incorporating this nonlinearity in the model. This was tested by comparing the correlation coefficient of each cell’s observed IFR with the linearly filtered signal · , versus the coefficient between observed IFR and the nonlinear predicted activity given by transforming · through f (Fig. 3, A and B). The mean difference between nonlinear and linear prediction accuracies for monkey M38 was small (0.003) and not significantly greater than zero (one-tailed t-test, P < 0.05). For monkey M36, the mean difference was also small (0.01) but significant (one-tailed t-test, P < 0.05) (Fig. 3, C and D).

View larger version (24K): [in this window] [in a new window]   FIG. 3. A: scatterplot of nonlinear vs. linear prediction accuracies for cells from monkey M38. Each point represents the correlation coefficient for a single cell between the observed and predicted IFR, for the 2 types of predictions. Diagonal line indicates unity: cells fall close to this line, indicating that incorporating the nonlinearity f did not significantly increase the prediction accuracy. Point types correspond to the 3 cell types sampled. Predictions were made using maximum filter length (see Multiple filter delays). C: histogram of differences between nonlinear and linear prediction accuracies. Mean difference is very small (0.003). B and D: same analyses for M1 neurons recorded in monkey M36. Mean difference = 0.01.

Figure 4A compares the encoding function f( · ) with the distribution of values of · (Fig. 4B) for a PTN from M1. The small improvement offered by incorporating the nonlinearity into our model of the cell—despite the fact that the shape of the encoding function is clearly nonlinear—can be explained by the restricted range of · . Most values of · are confined to the linear portion of the full encoding function and thus the behavior of the cell in relation to muscle activity is effectively linear.

View larger version (17K): [in this window] [in a new window]   FIG. 4. A: example nonlinear encoding function f( · ) = p(spike| · ) for a PTN from M1. This function gives the conditional firing rate (y-axis) for each value of the filtered signal · (x-axis), after applying the nonparametric binning process outlined in Analysis. B: although the form of f( · ) is nonlinear, examining the distribution of values of · observed for this neuron indicates that most of the signal was restricted to a narrow range corresponding to the linear portion of the encoding function. Therefore the dependency of the firing rate of this cell on the EMG signal filtered by its weight vector was effectively linear.

Comparison with corticomotoneuronal cells

The greater accuracy of fits to PTNs versus UIDs could be explained by the fact that these cells are more directly involved in the control of muscle activity. To test this hypothesis further we analyzed fits made to nine of the subset of 16 PTNs recorded in monkey M36 that were shown to be CM cells, using the maximum filter length for highest fit accuracy (see Multiple filter delays below). This group of nine cells showed only postspike facilitation (not suppression) of one or more of the seven EMGs recorded in this monkey. Methods for identification of genuine postspike effects are described in Jackson et al. (2003). Nonlinear prediction accuracies were compared with those from a separate collection of nine PTNs from M36, which showed no significant postspike effects in any of the EMGs. Of this group, four PTNs came from the original group of 22 PTNs in M36 and an additional five PTNs from extra data sets were included purely for this comparison and are not analyzed elsewhere. For all cells, models were fitted at a single = 40 ms using all seven EMGs recorded in this monkey.

Mean prediction accuracies for the two cell types were as follows. For the nine PTNs with no postspike effects, mean accuracy was 0.26 (SD 0.24). Mean prediction accuracy for the nine PTNs showing postspike facilitation was 0.61 (SD 0.17). This difference was significant (one-tailed t-test, P < 0.01). Also for the nine CM cells showing postspike facilitation there was a weak positive relationship between the size of the muscle field and prediction accuracy for the cascade model, which was not significant (R² = 0.41). The more accurate fits of CM cells to EMG in our model is consistent with the fact that the discharge of these cells directly influences the activity of one or more of the muscles analyzed, through the precise time locking of a proportion of CM cell spikes to motoneuron discharge. Similarly, CM cells showing postspike effects in a larger proportion of the muscles analyzed (i.e., cells with larger muscle fields) would be expected to be fitted more accurately.

Furthermore, we looked at the relationship between the pattern of muscle weighting given by and the pattern of postspike effects in the same muscles, for a subset of 7/16 CM cells from M36. Three of these cells showed a good degree of correspondence between the two, so that muscles with large positive weights in also showed significant postspike facilitation from the cell, and muscles with negative weights showed significant postspike suppression. These neurons are akin to a group of CM cells ("set A") reported by Bennett and Lemon (1996) in which the pattern of postspike effects and cell–muscle covariation would act together to promote a fractionated pattern of muscle activity important for the performance of precision grip.

In contrast, four cells showed a poor overlap between and the pattern of postspike effects. Although the postspike effects of these CM cells would also tend to fractionate activity, this was not reinforced by the weighting of muscle activity in , and so these cells are similar to the "set B" neurons described by Bennett and Lemon (1996).

Spike-triggered covariance

Earlier in the analysis section, it was assumed that the spike rate depends on a single dimension of the input signal : that is, the amplitude of this signal. The relationship of cell discharge to the level of EMG activity in multiple muscles depended on a weight vector or "linear filter," . However, the spike rate of many neurons (e.g., cells in primary visual cortex, V1) is best related to not one but multiple dimensions of the input signal. To account fully for the dimensionality of the signal to which these cells respond therefore requires multiple linear filters, (Adelson and Bergen 1985; Rust et al. 2005). Here, we examined whether single neurons in M1 and dentate nucleus encode multiple dimensions of the EMG signal, by testing whether the relationship of cell spiking to EMG was better described by multiple filters instead of only one. To do so, we applied the following analysis.

The different dimensions along which varies can be thought of as dimensions in a vector space, and each possible value of as a vector in this space. Returning to the basic LN model, the firing rate depended only on the projection of this EMG vector along a single direction in this space, .

However, if in reality the cell responds to multiple dimensions of the signal, then the firing rate would depend not on one but on multiple vectors in the stimulus space, corresponding to the multiple linear filters . These multiple vectors can be estimated using spike-triggered covariance (STC) techniques, full details of which can be found in previous studies (Brenner et al. 2000; Simoncelli et al. 2004).

For each cell this involved computing the covariance matrix C_spike, by taking all stimulus vectors s (i.e., values of ) that were conditional on the occurrence of a spike (at time t_spike) according to

(4)

We next computed C_prior, the covariance matrix formed from all stimulus vectors s at all times t. C_prior is the covariation matrix of the EMG input with itself

(5)

We then calculated C = C_prior – C_spike. The eigenvectors E of the matrix C correspond to linear combinations of the STC vectors or linear filters: the STC vectors were then extracted by computing (C_prior^–1)E. We selected two STC vectors giving two linear filters , , by taking the two eigenvectors that were most different from zero. The LN model was then applied in the usual manner except that now the response of the cell modeled in Eq. 1 was given by the filtering of the EMG input by two filters

(6)

According to Eq. 5, the cell firing rate is dependent on a two-dimensional (2D) nonlinear encoding function f(₁ · , ₂ · ) applied to the two filtered signals ₁ · and ₂ · . We computed 2D encoding functions for each cell from M38. Figure 5A shows an example for a representative cell in M1: contours of firing rate modulation are approximately linear in the region where most stimulus vectors conditional on spiking were distributed (gray points in Fig. 5B), suggesting that the relationship of cell firing to muscle activity is sufficiently captured by the application of a single filter. To test the accuracy of this new model in predicting cell activity we mapped the filtered input signal (₁ · , ₂ · ) through the corresponding 2D encoding function and compared the predicted firing rate to the observed IFR by means of the correlation coefficient as above. For all cells, these predictions were less accurate as a result of overfitting; computing two filters and the resulting 2D encoding functions for each cell involved fitting a greater number of regressors to the training data, which decreased the accuracy of predictions of the test data. This suggests that a single linear filter (i.e., the spike-triggered average) is sufficient to describe the encoding of muscle activity by M1 and dentate neurons.

View larger version (37K): [in this window] [in a new window]   FIG. 5. Spike-triggered covariance (STC) analysis of a PTN from M1. A: example 2-dimensional (2D) nonlinear encoding function f(1 · , 2 · ) = p(spike|1 · , 2 · ) giving the conditional firing rate for each value of the 2 filtered signals 1 · (x-axis) and 2 · (y-axis). Color axis indicates the firing rate (in Hz) conditional on these 2 variables. Contours of firing rate appear to be curved toward the right-hand side of the plot, suggesting that conditional firing rate depends on both filtered signals. Masking is applied to the regions where insufficient data are observed for accurate estimation of the firing rate. B: black points show values of the EMG signal filtered by 2 corresponding filters of the cell 1 · and 2 · plotted against each other. Gray points show values of these filtered signals that were conditional on spikes fired by the cell. Note that although the form of f[1 · (t), 2 · (t)] shows modulation of firing rate in 2 dimensions, these gray points are restricted to a region toward the left of the 2D encoding function in A where contours of firing rate are approximately linear, showing that spiking activity of the cell is mainly modulated by 2 · . Thus a single filter is sufficient to describe the encoding of muscle activity by this neuron.

Prediction accuracies were calculated as described above, for different single delays between spike and EMG (ranging from –160 to +320 ms in 10-ms increments), to assess the temporal evolution of tuning for muscle activity in these neurons. Because linear and nonlinear predictions were found to be equivalent, we refer to firing rate predictions made using f( · ) unless otherwise stated. Different cells had a different value of that gave the best fit (Fig. 6A). For example, the tuning of the cell given by the topmost curve reaches peak prediction accuracy at close to = 0 s, whereas the curve below peaks at about = 0.18. There was no consistent optimum across cells.

View larger version (28K): [in this window] [in a new window]   FIG. 6. Temporal tuning functions for the encoding of muscle activity. A: tuning curves for single cells. Each curve corresponds to the nonlinear prediction accuracy as a function of for a single cell. Asterisk indicates optimal lag. Only a randomly chosen subsample of the population is shown to avoid overcrowding. B: average temporal tuning functions for the population, divided into 3 cell types (see legend). Each trace represents the mean change in nonlinear prediction accuracy as a function of for each cell group. Error bars indicate SE.

Overall, temporal tuning curves were broad (Fig. 6A). These broad curves were not an artifact of the filtering and smoothing that we applied to the EMG data. The total filter length (as measured by applying the same preprocessing to a signal comprising a spike in a single bin) was 78 ms, shorter than the time scale of the correlations we observed. We also confirmed that this broadness was not simply a result of combining multiple EMGs that were acting at different times relative to each other during task performance because the same broad tuning curves were found for correlations with single muscles.

Next, we analyzed the shapes of these temporal tuning curves in more detail. All cells from monkey M36 were included. In monkey M38 less-accurate task performance may have resulted in greater trial-to-trial variability of cell firing rates and EMG, reducing the strength of correlation between cell discharge and muscle activity (mean nonlinear prediction accuracy at 40-ms lag for M1 cells from M38 was 0.38, compared with 0.48 for neurons sampled in M36). To reduce the effect of noise on the shapes of temporal tuning curves we therefore selected the most accurate 50% of the cells in each group from M38.

Dividing cells into three groups—PTNs, UIDs, and dentate—revealed that there were no significant differences either in the peak lag distributions between the three groups (two-tailed Kolmogorov–Smirnov test between each distribution pair, P < 0.05) or in the mean peak lag (two-tailed t-test, P < 0.05). Figure 6B plots the mean change in prediction accuracy with lag for the population. PTNs and dentate units reached maximum prediction accuracy at 10-ms lag, whereas the peak lag for UIDs was 30 ms. The peak lag for PTNs corresponds to values obtained from previous regression analyses (Morrow and Miller 2003) and spike-triggered averaging work (Fetz and Cheney 1980). At all values of , PTNs were on average more strongly correlated with EMG than with UIDs, which showed stronger correlations than dentate units.

Finally, we looked at temporal tuning curves for 16 PTNs from M36 identified as CM cells, once again using the concatenated activity of all seven EMGs recorded from this monkey. Here also, the shapes of temporal tuning curves for single cells were broadly curved. Optimal lags were measured from these curves, which again were heterogeneous in their distribution. The level of heterogeneity for CM cells was not significantly different from that for other cells. For a fair comparison against cells showing no CM connections, we compared the variance in optimal lags measured for CM cells with those for the top 50% of dentate cells in M38, using Levene’s test for homogeneity of variance—this was not significant (P > 0.05). Thus CM cells showed temporal profiles of tuning for muscle activity similar to those for the other cell types in our model.

Multiple filter delays

Prediction accuracies were compared for different lengths of the linear filter using 40-ms increments from 40 to 360 ms after each spike. Again, there was some heterogeneity in the optimal filter length (Fig. 7A). Figure 7B shows the mean change in prediction accuracy with increasing filter length for all cells, once again taking all cells from M36 and the top 50% of each subpopulation of neurons from M38. For all three cell types, the mean correlation coefficient increased smoothly with increasing length but the range of this increase was small. PTNs showed a larger range of increase in correlation coefficient values than that in dentate units. Looking at the total population (without selecting the most accurate cells), at all filter lengths the highest prediction accuracies were achieved using PTNs, followed by UIDs, and then by dentate units. At the maximum filter length, M1 neurons together were nearly twice as accurate as dentate cells.

View larger version (19K): [in this window] [in a new window]   FIG. 7. Prediction accuracy as a function of number of delay samples. A: modulation for single cells. Each trace corresponds to the nonlinear prediction accuracy vs. number of delays for a single cell, where the delay increment is 40 ms. Asterisk indicates optimal number of delays. A random subsample of cells is shown to avoid overcrowding. B: population mean prediction accuracy as a function of number of delay samples, divided into 3 cell types (see legend). Each trace represents the mean change in nonlinear prediction accuracy vs. number of delays. Error bars indicate SE.

Kinematic information

LN models were fitted to neurons from M36 and M38 using finger position and velocity information recorded during task performance, for the same single lag and multiple lag increments described above. We compared these fits to correlations of cell discharge with muscle activity. Before making this comparison of fits to different movement parameters, trials were selected for analysis using the criteria described in METHODS to ensure that trial performance was as homogeneous as possible. In monkey M38, three PTNs were subsequently excluded from analysis because they fired mainly during trials that were rejected.

Correlation coefficients for nonlinear versus linear predictions were compared (Fig. 8A). In contrast to the analysis conducted with EMG data, cascade models using finger position and velocity information showed an increase in the strength of the nonlinearity: more cells were above unity (M38: 76% for kinematic fit, compared with 69% for EMG fit; M36: 95% for kinematic fit compared with 70% for EMG). The mean difference in prediction accuracy between nonlinear and linear predictions was greater for fits to kinematic data: in M38 it was 0.03 for kinematic fits compared with 0.003 for EMG; in M36 it was 0.03 for kinematic fits compared with 0.01 for EMG. However, in both monkeys this difference did not reach significance (one-tailed t-test, P > 0.05). The increased effect of the nonlinearity for fits to kinematic variables is consistent with a previous study of M1 neurons during an arm-reaching task, where kinematic encoding was found to be significantly nonlinear, although the effect of this nonlinearity was also small (Paninski et al. 2004c).

View larger version (17K): [in this window] [in a new window]   FIG. 8. Comparing the encoding of muscle activity vs. kinematic information. A: scatterplot of nonlinear vs. linear prediction accuracies using finger position and velocity information. Conventions are the same as for Fig. 3: note that points fall further above the diagonal in the current figure, suggesting a greater increase in prediction accuracy for encoding of kinematic information when incorporating the nonlinearity f. B: scatterplot comparing prediction accuracies for muscle activity (y-axis) vs. kinematic information (x-axis). Each dot represents the correlation coefficient for a single cell between the observed and predicted IFR, obtained using the 2 types of information. Diagonal line indicates unity. Correlation coefficients do not fall significantly either side of diagonal, indicating equal accuracy of the model when fitted to muscle activity or kinematic information. C: scatterplot of prediction accuracies for a combination of muscle activity and kinematic information (y-axis) vs. muscle activity alone (x-axis). Correlation coefficients fall significantly above diagonal, indicating greatest accuracy of the model when fitted to both types of information. All plots show data from M36 and M38.

Finally, fits made using kinematic and EMG information in M36 and M38 were compared across all neurons. There was no significant difference in mean prediction accuracy between kinematic and EMG information (paired t-test, P > 0.05) (Fig. 8B). Combining both kinematic and EMG information increased prediction accuracy above what was observed with either EMG or kinematic data alone (Fig. 8C). In all of these instances, PTNs were always modeled with greater accuracy compared with UIDs and dentate neurons.

Predictions of spiking activity

We compared the ability of the Poisson and spike history (SH) models described above to predict novel spike trains, in the population of cells from monkey M38. Fitting and f to training data that incorporated the spiking history of the neuron improved the dynamic range of the encoding function (Fig. 9A). For all 94 cells, both the absolute and percentage modulation of firing rate measured from the function were significantly higher when cells were fitted with the SH terms compared with EMG alone (one-tailed t-test, P < 0.05). However, despite this apparent improvement, cross-validated predictions of the observed IFR made by mapping the filtered signal · through the encoding function f were significantly less accurate for the SH model, compared with when cells were fitted with the EMG alone (one-tailed t-test, P < 0.05). The increase in dynamic range captured by incorporating the SH terms failed to improve fit quality, resulting from the fact that, as with the basic LN model, the rising part of the curve still corresponded to the tail of the distribution of the filtered signal, where there were few data (Fig. 9B). In turn, fits made using SH terms did not lead to improved performance of the model when predicting the spike count in a novel, cross-validated section of the data (Fig. 10). Because the iterative spike generation process was fairly time consuming, a subset made up of the 10 cells that gave the most accurate nonlinear predictions of the IFR was tested. For each cell, the cross-validated average firing rate computed from 20 repeats of the spike generation process was compared with the observed IFR by means of the correlation coefficient, for the LN and the SH models. The SH model was not significantly more accurate at predicting the spike count (one-tailed t-test, P > 0.05).

View larger version (18K): [in this window] [in a new window]   FIG. 9. Incorporating spike history (SH) effects. A: effect of including SH terms on dynamic range of cell encoding function. Black trace: example encoding function f, for a PTN from M1, for LN model fitted to muscle activity as in Fig. 4. Gray trace: encoding function for same cell fitted to both muscle activity and the previous response history of the cell 10 ms previously. Note that although the dynamic range of modulation increases, the curve shifts right to higher values of the filtered signal · . Error bars indicate SE. B: effect of SH on filtered signal. Black histogram shows distribution of filtered signal values for fit to muscle activity as in Fig. 4. Gray histogram shows distribution of filtered signal values for fit incorporating SH, which shifts to the right. Thus for both fits, much of the change in the firing rate of the cell takes place at the tail of the corresponding distribution, where there are few data. Filtered signal values have been scaled to –1 and +1 as described in METHODS.

View larger version (31K): [in this window] [in a new window]   FIG. 10. Predicting spikes from muscle activity. Figure shows same cell as Fig. 9. A: raster display showing 20 iterations of predictions of spiking activity for an M1 neuron using the basic linear–nonlinear (LN) model. B: equivalent raster display for predictions made using the spike history model described in Analysis. C: averages of the above rasters (green and red traces, respectively) compared with observed IFR. Both types of predictions do about as well as each other at predicting the periods of spiking activity. D: observed spikes shown for comparison. Predictions are illustrated for a small, randomly chosen portion (3 s long) of the total data recorded from this neuron.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES