Stable Ensemble Performance with Single-Neuron Variability during Reaching Movements in Primates

Introduction

Numerous studies have shown that cortical neurons modulate their firing in association with spatial and temporal aspects of behavior (Fetz, 1992; Johnson et al., 2001). In recent years, the idea that neuronal modulations can be used to control a brain-machine interface (BMI) has received significant attention (Wessberg et al., 2000; Nicolelis, 2001; Taylor et al., 2002; Carmena et al., 2003; Nicolelis, 2003; Musallam et al., 2004; Lebedev et al., 2005). There are, however, different opinions in the literature about the size of the neuronal sample needed to operate a BMI efficiently. One principal challenge in building useful BMIs is coping with variability of neuronal firing (Zohary et al., 1994; Abbott and Dayan, 1999). Cortical neurons do not represent any motor parameter precisely (Ashe and Georgopoulos, 1994; Sergio and Kalaska, 1998; Carmena et al., 2003). Instead, neuronal firing fluctuates from one behavioral trial to another (Lee et al., 1998; Shadlen and Newsome, 1998; Cohen and Nicolelis, 2004). This variability in neuronal firing and its correlation with behavioral parameters has recently been explained by a model of neural computation based on perturbations (Mass et al., 2002). This model suggests that stable internal states are not required for a stable output in highly dimensional dynamic systems like the brain.

Here, we report significant variability in firing properties of individual neurons while two monkeys, chronically implanted with multielectrode arrays in frontal and parietal cortical areas, performed a continuous arm movement task. The activity of a single neuron was typically correlated with several motor parameters. However, the strength of a given correlation was generally nonstationary within a 30-60 min recording session. Moreover, although the correlation between the firing of single neurons and movement parameters was variable, stable predictions of arm movements could still be obtained from the activity of the whole neuronal ensemble recorded in each session. Thus, although significant variability was observed in the firing properties of individual neurons, the ensemble performance in the BMI control remained stable.

This result adds to the recent literature that demonstrates redundant encoding by neuronal ensembles (Laubach et al. 2000; Schneidman et al., 2003; Cohen and Nicolelis, 2004; Paz and Vaadia, 2004; Narayanan et al., 2005) and shows that brain networks can achieve the same behavioral goals through different patterns and relative contribution of individual neuron activity. The fact that the brain may use different ensembles of neurons to perform a given task may impact the design and implementation of a robust brain-machine interface that can maintain a stable performance over time.

Materials and Methods

Behavioral training and electrophysiology. Two adult female monkeys (Macaca mulatta) were used in this study. All procedures conformed to the Guide for the Care and Use of Laboratory Animals and were approved by the Duke University Institutional Animal Care and Use Committee. Before the surgical implantation of recording arrays, the animals completed a period of chair training. The arrays contained 16-64 microwires each (separation between adjacent microwires, 300 μm) and were implanted in frontal and parietal cortical areas. Implanted areas specifically included the dorsal premotor cortex (PMd), supplementary motor area (SMA), and the primary motor cortex (M1) in both hemispheres. In monkey 1, an additional implant was placed in the primary somatosensory cortex (S1). In monkey 2, the medial intraparietal area (MIP) of posterior parietal cortex (PP) was also implanted. Both monkeys performed the task with their left arms. After recovery from the surgical procedure, the animals were engaged in behavioral training. Neuronal recordings commenced at the same time. A 512 Multichannel Acquisition Processor (Plexon, Dallas, TX) simultaneously recorded single-neuron and multiunit activity during each recording session. Spike sorting was performed as described previously (Nicolelis et al., 2003). In this study, only single units that maintained stable waveforms across the analyzed recording sessions were used. Monkeys were seated in a primate chair facing a computer monitor. They performed a motor task using a handheld pole. The position of the monkey's hand was obtained from an infrared marker located on top of the pole. The marker was monitored by an optical tracking system (Optotrak; Northern Digital, Waterloo, Ontario, Canada). In the task, the monkeys were shown a small disk (the “cursor”) and a larger disk (the “target”). This is illustrated in Figure 1 Ai. Each trial began with a target presented at a random position on the screen. The monkey had to use the pole to put the cursor over the target for at least 150 ms. If the monkey crossed the target too fast, the target disappeared and the trial was not rewarded. The monkeys had 5 s to hit the target, with a delay of 0.5 s in between trials. They received juice rewards for correct performance.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Experimental design and examples of stable and variable neurons. A, Schematics of the behavioral task and neuronal encoding analysis. The monkey controls a handheld pole while seated in front of a computer screen (i). The goal of the task is to move the cursor (small disk) to a visual target (large disk) that appears at random locations on the screen for each trial. Multidimensional linear regression is applied off-line for linearly mapping neuronal ensemble activity with motor output (ii). Independent models are generated for hand position and hand velocity coordinates. For each sample (bin size, 100 ms), 10 lags (i.e., 1 s total) before movement are assigned as free parameters to the linear model. B, Examples of stable and variable neurons during a 22 min recording session of continuous movement (first 10 min used for fitting the first model are not shown). Single-unit waveforms from PMd (top) and M1 (bottom) neurons recorded during the same session (i). Numbers in vertical axes depict the mean V_pp (in microvolts) of the action potentials. Interspike intervals (ISIs) of the depicted neurons are shown (ii). Evolution of the IRFs obtained by regressing a single unit with each motor parameter are shown (iii). Models were generated with data points from a 10 min sliding window in 30 s increments. Each model consisted of 10 free parameters from the corresponding time lags. The color bar denotes IRF values for each model. Note the different IRF patterns across time between a stable (top) and variable (bottom) neuron. Evolution of the single-unit ranking, single-unit R score, and ensemble R score across the recording session are shown (iv). For plotting purposes, neuron ranking is normalized from 0 (worst) to 1 (best) from the 94-neuron optimized ensemble. Notice the abrupt changes in both ranking and R score of the variable (bottom) neuron.

Linear model. Hand position, velocity, and gripping force were modeled as a weighted linear combination of neuronal activity using a multidimensional linear regression or Wiener filter (Haykin, 2002) described by the following equation: (1)

In this equation, x(t - u) is an input vector of neuronal firing rates preceding time t by time-lag u, y(t) is a vector of kinematic variables (e.g., hand position and velocity) at time t, a(u) is a vector of weights for time-lag u, b is a vector of y-intercepts, and ϵ(t) are the residual errors. The lags can be positive (i.e., past) or negative (future). Here, we considered only lags in the past. This equation can be recast in matrix form as follows: (2)

where each row in each matrix corresponds to a unit of time and each column is a data vector. Matrix X contains lagged data, and thus has a column for each lag multiplied by each neural channel (e.g., 100 channels by 10 lags to a total of 1000 columns). The y-intercept is incorporated by prepending a column of ones to matrix X. Matrix A is solved by the following equation: (3)

Neuronal firing rates were sampled using 100 ms bins. Ten bins preceding time t were used for training the model and consequent predictions. Thus, in Equation 1, m and n equal 0 and 9, respectively. Models were trained with 10 min of data and tested on the subsequent records.

Data analysis and variability measurements. Data from 10 recording sessions for each monkey were used in this study. For each session, a linear model was generated for each of the motor parameters investigated in this study (i.e., X and Y coordinates of hand position and velocity) (Fig. 1 Aii). Models were fitted using training data sets of 10 min with a 30 s step sliding window. All of the models generated in a given session were tested against a data set of 5 min drawn from the end of the session. For assessing neuronal variability, firing rates of individual neurons were regressed independently and subsequently tested for prediction of each motor parameter. As with ensemble regression, regressions for individual neurons used 10 time lags, resulting in 10 free parameters as inputs to the model. Dividing by the SD before conducting linear regression normalized firing rates of each neuron. Prediction quality was evaluated as the correlation coefficient R between the 5 min test data set and the predictions from the model. The values of R were then ranked from the highest (first rank) to the lowest. In addition to R-based ranking, neuronal sensitivity was used to describe the contribution of neurons to ensemble predictions. Neuronal sensitivity was calculated as the average of the absolute value of the impulse response functions (IRFs) obtained from the regression for the whole ensemble (Sanchez et al., 2004). Both methods produced similar results (r = 0.88). Significant variability in the contribution of individual neurons to prediction of motor parameters was detected as the presence of fluctuations in ranking across a recording session that exceeded 3σ. If the contribution rank of neuron i at time t exceeded 3σ from the mean of the distribution of the remaining ranking points in that session, the neuron was classified as variable. The fact that this selection criterion is not sensitive to smooth ranking transitions across the session, but rather is only affected by abrupt decreases or increases in rank value, makes it very conservative. The same criterion was applied to measure variability in the mean firing rate of individual neurons.

Cumulative and residual neuron ranked curves. Ten recording sessions were used for this analysis. Prediction capacity of each neuron in the ensemble during each session was ranked as described above. For the cumulative ranked curves, models were trained (10 min epochs) and tested (5 min epochs) for each motor parameter. The size of the population was increased in unitary increments, starting from one (best) single unit and ending with the total size of the ensemble. For the residual ranked curves, the procedure was the inverse (i.e., starting with the whole ensemble and removing one neuron at a time starting with the neurons with the highest ranks and ending with the neurons with the lowest ranks). The final cumulative and residual curves represented the average for all sessions.

Results

In this study, we examined the relationship between neuronal activity and hand movements. We used 94 simultaneously recorded units in monkey 1 and 38 in monkey 2, all of which had stable action potentials within all recording sessions. Neuronal variability was quantified in terms of significant fluctuations (p < 0.01) of correlation between neuronal activity and parameters of hand movements. In the neuronal ranking analysis, this correlation was compared for different neurons in the ensemble.

We first examined firing patterns of individual neurons recorded from different cortical areas during several days of the monkeys' performance in the task. Individual neurons were classified in two main groups depending on the degree of variability in their encoding characteristics in each recording session. If the ranking of a neuron fluctuated significantly (p < 0.01) during the session, the neuron was classified as variable; otherwise the neuron was classified as stable. Figure 1B shows representative neurons from each group recorded in the same recording session in monkey 1. The evolution of the properties of each neuron during the course of the session is represented using the IRFs from the linear model obtained by regressing neuronal activation with hand position (y-coordinate) (Fig. 1Biii). Whereas the IRFs of the stable neuron remained very similar across the session, the IRFs of the variable neuron changed dramatically. The result of testing the models (i.e., the sum of products between IRFs and the spike trains of their corresponding neurons) is shown in Figure 1Biv. For the stable neuron, none of the curves representing the ranking and contribution of the neuron (R score) changed significantly over time. The variable neuron (Fig. 1Biv, bottom), however, showed dramatic changes both in contribution to the model (blue line) and in ranking (green line) during the first 3 min of the session, whereas the ensemble encoding remained stable. Notice the large fluctuations in ranking of the neuron; it was the worst performer in the beginning of the session and became the best by the end of the session.

Next, we examined the spatiotemporal patterns of single-neuron variability based on ranking changes across the whole ensemble in monkey 1 (Fig. 2A). Individual neuronal activity was regressed to hand velocity (X component) in a recording session from monkey 1. Figure 2Ai shows neuronal variability in ranking (1, best; 94, worst) across all cortical areas, including M1, which was the best area for encoding movements. Sorting the ensemble by ranking in the first epoch (Fig. 2Aii) and plotting significant variable epochs (p < 0.01) for each neuron (Fig. 2Aiii), shows that only the best four neurons of the ensemble remained relatively stable across the whole session. Up to 10% of the total population could simultaneously exhibit significant variability in a given epoch without affecting ensemble predictions (Fig. 2Aii, top mono-dimensional plot). For a given neuron, up to 16% of the epochs within a recording session were classified as significantly variable (Fig. 2Aii, right mono-dimensional plot). Neurons that provided a high contribution to a parameter in a given epoch were not necessarily the best contributors for predictions of the same motor parameter during the other epochs. This effect was observed for all cortical areas. Variability in a given session was quantified as a percentage of neurons of the ensemble that showed changes of >3 SDs at least once during the session. Figure 2B shows the variability of the rankings that resulted from predictions of all motor parameters and the variability of the mean firing rate of the individual neurons across four consecutive sessions. In monkey 1, variable neurons reached 65% of the ensemble size for ranking and 55% for firing rate. In monkey 2 (data not shown), these values reached 54 and 50%, respectively. Although individual neurons exhibited dramatic changes, overall ensemble performance did not decrease significantly. Ensemble performance was stable because of the compensation of other neurons contributing to the linear-model predictions. The compensation is revealed in Figure 2C, which shows redundancy of information in the ensemble. In the figure, cumulative curves show that the optimal linear mapping between spike trains and motor output can be obtained using the top 50 neurons and, more remarkably, that ∼90% of the accounted variance is obtained with only the first top 10 neurons. In residual curves (i.e., the performance of the remaining ensemble after individually removing neurons in ranked order), the figure shows the rest of the ensemble containing a considerable amount of useful information. For example, after removing the best 50 neurons, 60% of the total ensemble contribution remained in the residual population, suggesting that the majority of the remaining ensemble was not composed of noisy, task-unrelated neurons but rather of highly modulated task-related neurons.

Finally, Figure 2D illustrates the importance of recording from large neuronal ensembles to maximize and stabilize model performance. The figure compares performance, for each epoch, between ensembles of 10 randomly selected neurons (mean R ± SD; 15 samples) and ensembles of the top 10 neurons of a total population of 94 units. The large size of this population increases the chance of yielding neurons with large contributions for a given epoch, resulting in significantly higher and more stable performance than with small populations as represented by the randomly chosen ensembles.

Discussion

The analysis performed in this paper was motivated by the observation that the neurons that provided the highest contribution to a given parameter during a given epoch did not necessarily provide the highest contribution to the same parameter during other epochs. Additional analysis revealed that, although there was significant variability in neuronal properties, the performance of the whole ensemble was stable. One explanation for the ensemble stability is that the change in performance of one neuron was compensated by other neurons. Another explanation is that the neuronal variability is simply biological noise. The stable activity of the ensemble might then arise because the noise contributed by more variable neurons remains approximately constant across the ensemble, although the contribution of individual neurons changes stochastically. Although the underlying basis of this neuronal variability is unknown, one possibility is that the neuronal population may have a bias for encoding other potential variations of movements that are not exercised in the current paradigm. In this regard, additional experiments with behavioral tasks containing sequential movements in different contexts (Ben-Shaul et al., 2004) will be needed to further investigate this phenomenon.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Stable ensemble performance with single-neuron variability: spatiotemporal patterns, quantification, and implications of redundancy in neuronal ensembles. A, Spatiotemporal patterns of single-neuron variability based on ranking changes during a recording session. A color map shows an ensemble of neurons arranged by cortical area (i). The color bar denotes ranking (1, best; 94, worst) and the range of individual correlation coefficient scores for the whole ensemble. Arrows depict the top four neurons. A color map shows the same ensemble of neurons as in i but arranged by their initial ranking at the first epoch (ii). The mono-dimensional plot on top shows stable ensemble performance. A color map (red color) depicts statistically significant variability (p < 0.01) of a particular neuron at a particular epoch (iii). The ensemble is arranged as in ii. Mono-dimensional plots show the amount of significantly variable neurons at a given epoch (top) and the amount of epochs for a given neuron that were classified as variable within the recording session (right). B, Overall quantification of neuronal variability across all motor parameters and for four consecutive recording sessions. Variability was quantified as statistically significant (p < 0.01) changes in neuron ranking and firing rate. The numbers in parentheses indicate the size of the neuronal population in a particular cortical area. C, Cumulative (solid) and residual (dotted) neuron-ranked curves for position (X, Y) and velocity (VX, VY). Notice the smooth decay of the residual curves showing the high contribution of most of the neurons of the ensemble. D, Performance comparison between ensembles of the best 10 neurons in a given epoch and ensembles of 10 randomly selected neurons, of a total population of n = 94 units. For each epoch, 15 ensembles of randomly selected neurons were estimated (mean R ± STD).

It is tempting to speculate that neuronal variability is beneficial for neuronal computations (for example, fault tolerance in a highly redundant system like the brain). In fact, the redundancy of encoding by neuronal ensembles that we observed concurs with recent reports (Schneidman et al., 2003; Ben-Shaul et al., 2004; Cohen and Nicolelis, 2004; Narayanan et al., 2005). As shown in Figure 2C, the redundant information that remains after removal of the best performing neurons is depicted by the smooth decay of the residual curves. This redundant information allows the ensemble to maintain optimal performance across time and supports our previous contention (Wessberg at al., 2000; Carmena et al., 2003) that reliable task performance in BMIs will require recording from large neuronal ensembles. Specifically, by sampling over a large population of neurons, the effect of individual neuronal variability on the ensemble output will be minimized. Moreover, by recording from many neurons, the BMI will be able to select a subpopulation of neurons with the greatest contribution to achieve optimal task performance in the control of a prosthetic device at a given moment. Consistent with the results reported here, Kurtzer et al. (2005) reported random shifts in load representation during behavioral tasks that involved control of posture and reaching movement. The authors argue that current BMI approaches that assume fixed correlation between neuronal activity and motor output may be suboptimal. Therefore, by recording from large neuronal ensembles, it will be possible to overcome variable representations by fitting separate models for specific components of the behavioral task. It should also be noted that recording from many neurons also protects the performance of the BMI from the loss of neuronal activity to electrode mal-function, modification of electrode tip properties, or cell death. BMIs that rely on only a very small sample of neurons (Serruya et al., 2002; Taylor et al., 2002) are not likely to accommodate the alterations in the neuronal sample over time.

In conclusion, this study showed that, although individual cortical neurons can exhibit significant variability in both firing properties and correlation with a motor task within an ensemble, the ensemble can achieve the same behavioral goals through different patterns and the relative contribution of the activity of individual neurons. This finding reinforces the importance of recording from large neuronal ensembles to achieve reliable BMI performance in the presence of individual neuronal variability.