Abstract
Monkeys can learn to directly control the movements of an artificial actuator by using a brain-machine interface (BMI) driven by the activity of a sample of cortical neurons. Eventually, they can do so without moving their limbs. Neuronal adaptations underlying the transition from control of the limb to control of the actuator are poorly understood. Here, we show that rapid modifications in neuronal representation of velocity of the hand and actuator occur in multiple cortical areas during the operation of a BMI. Initially, monkeys controlled the actuator by moving a hand-held pole. During this period, the BMI was trained to predict the actuator velocity. As the monkeys started using their cortical activity to control the actuator, the activity of individual neurons and neuronal populations became less representative of the animal's hand movements while representing the movements of the actuator. As a result of this adaptation, the animals could eventually stop moving their hands yet continue to control the actuator. These results show that, during BMI control, cortical ensembles represent behaviorally significant motor parameters, even if these are not associated with movements of the animal's own limb.
Introduction
Ever since Evarts (1966) pioneered single-unit recordings from the cortex of awake, behaving monkeys, numerous studies have scrutinized a fundamental observation: cortical neurons modulate their firing during voluntary movements and thereby encode a variety of motor parameters (Georgopoulos et al., 1986, 1992; Ashe and Georgopoulos, 1994; Fu et al., 1995; Sergio and Kalaska, 1998). This encoding incorporates motor, sensory, and cognitive signals and appears to be highly complex (Cheney and Fetz, 1980; Alexander and Crutcher, 1990; Lebedev et al., 1994; Shen and Alexander, 1997; Donchin et al., 1998; Kakei et al., 1999) (for review, see Georgopoulos, 2000; Johnson et al., 2001; Andersen and Buneo, 2002; Paz et al., 2004). Individual cortical neurons do not have a one-to-one relationship to any single motor parameter (Ashe and Georgopoulos, 1994; Sergio and Kalaska, 1998; Carmena et al., 2003), and their firing patterns exhibit considerable variability (Lee et al., 1998; Shadlen and Newsome, 1998; Cohen and Nicolelis, 2004; Wessberg and Nicolelis, 2004). Precise motor control is achieved through the action of large neuronal ensembles (Wessberg et al., 2000; Carmena et al., 2003). With the advent of electrophysiological techniques permitting simultaneous recordings from hundreds of neurons in multiple brain areas (Nicolelis et al., 1995, 1997, 2003; Kralik et al., 2001; Nicolelis and Ribeiro, 2002), detailed information can be collected about ensemble encoding.
Using ensemble recordings, motor parameters can be predicted from neuronal activity in real time and used to control a brain-machine interface (BMI) (Chapin et al., 1999; Wessberg et al., 2000; Nicolelis, 2001; Serruya et al., 2002; Taylor et al., 2002; Carmena et al., 2003; Wessberg and Nicolelis, 2004; Patil et al., 2004). In addition, BMIs can be driven by planning and motivation signals extracted from cortical activity (Musallam et al., 2004). Recently, we implemented a BMI in which cortical ensemble activity recorded from macaque monkey brain directly controlled reaching and grasping movements performed by a robotic arm (Carmena et al., 2003). The monkeys' proficiency in operating the BMI improved with training, accompanied by changes in neuronal directional tuning. Based on these observations, we suggested that, with prolonged usage, an artificial actuator can be effectively incorporated in the internal representation of the subject's body, commonly referred to as “body schema” (Head and Holmes, 1911; Gurfinkel et al., 1991; Maravita et al., 2003; Nicolelis, 2003).
Notwithstanding initial reports that training to operate a BMI evokes plastic changes in the activity patterns of cortical neurons (Taylor et al., 2002; Carmena et al., 2003; Musallam et al., 2004), a number of questions about cortical representation of artificial actuators remain to be addressed. Perhaps the most intriguing finding is that neuronal ensembles previously modulated in association with limb movements continue to be modulated after the subjects stop moving their limbs, and these modulations represent purposeful movements of the actuator (Chapin et al., 1999; Taylor et al., 2002; Carmena et al., 2003). In the present study, we hypothesized that neuronal adaptation to represent the velocity of an artificial actuator instead of the limb can occur even while the subjects continue to move their limbs. Such adaptations would function to optimize neuronal representation of new behavioral goals (Todorov and Jordan, 2002). Accordingly, we examined ensemble encoding of the hand and actuator movements before and after transition to BMI control. We discovered changes in spatial and temporal representation of movements during the BMI control that appeared to establish the representation of the actuator at the expense of the representation of the limb.
Materials and Methods
Behavioral task and brain-machine interface operations. Two adult female monkeys (Macaca mulatta) were used in this study. All procedures conformed to the National Research Council Guide for the Care and Use of Laboratory Animals (1996) and were approved by the Duke University Animal Care and Use Committee. The monkeys were trained to pursue a visual target on a computer monitor with a cursor, initially set to motion using a hand-held pole (see Fig. 1 A) (Carmena et al., 2003). The pole actually controlled the movements of the robot arm, which was not visible to the monkeys, and the cursor provided visual feedback of the position of the robot (see Fig. 1 A). The target appeared at random locations on the screen, and the monkeys had to place the cursor over it to obtain a fruit-juice reward [task 1 of Carmena et al. (2003)]. Multielectrode arrays consisting of 16-64 microwires spaced at 300 μm were implanted bilaterally in primary motor cortex (M1), dorsal premotor cortex (PMd), supplementary motor area (SMA), primary somatosensory cortex (S1), and posterior parietal cortex (PP). Neuronal activity was sampled and sorted using a multichannel acquisition processor (Plexon, Dallas, TX).
Experimental design and modes of operation. A, Experimental apparatus. The monkey was seated in front of a computer monitor on which visual stimuli were shown. It had to pursue a visual target (large circle) with a cursor (small circle). The monkey controlled the cursor by moving a hand-held pole (pole control). The pole actually controlled a robotic arm invisible to the monkey, and the cursor position on the screen reflected the robot's position. A linear model was trained to predict hand/robot velocity from neuronal ensemble activity recorded from the monkey's cortex. Then, the pole was disconnected, and the robot was directly controlled by the model's output (brain control). B, Schematics of movement trajectory, instantaneous velocity, and neuronal discharges preceding or succeeding an IVM. Neuronal rates were estimated using 100 ms bins placed at different lags relative to the IVM. C, Representative traces of the hand (black) and the robot (red) during pole control. D, Hand and robot traces during brain control. E1, E2, Time-dependent traces of the hand and robot position during pole control. F1, F2, Time-dependent position traces during brain control. G1, G2, Hand and robot velocity during pole control. E1, E2, Velocity during brain control.
As the monkeys performed the task by moving the pole, a linear model was trained to predict hand velocity in real time (mathematical description is given below by Eq. 3) (Wessberg et al., 2000; Carmena et al., 2003). We call this mode of operation pole control. After at least 30 min of operation in pole control, the model parameters were fixed. The robot (and the cursor) was disconnected from the pole and was directly controlled by brain activity. Robot velocity was set to the predicted values of the model. We name this mode of operation brain control, of which we distinguish two types. During brain control with hand movements (BCWH), a mode that immediately followed pole control in most experiments, the monkeys continued to move the pole. This mode of operation continued for 20-30 min, after which the pole was removed from the apparatus. At that time, the monkeys stopped moving their arms, as confirmed by EMG recordings (Carmena et al., 2003), but continued to operate the BMI. This mode of operation is termed brain control without hand movements (BCWOH). BCWOH was first attempted after the monkeys were trained in BCWH for 1 week; the monkeys were able to control the actuator in this mode of operation and obtain rewards on the very first attempt. Each recording session continued for 60-90 min, until the monkeys had had enough juice and became unmotivated to perform the task. In some sessions, to obtain longer records of BCWOH, the pole was removed immediately after pole control.
Tuning to velocity during pole and brain control. Significant correlation between neuronal firing rate and movement velocity has been documented in previous studies of M1 (Georgopoulos et al., 1988; Ashe and Georgopoulos, 1994; Moran and Schwartz, 1999; Carmena et al., 2003) and other cortical areas (Carmena et al., 2003; Averbeck et al., 2005). This finding was key for the experiment reported here, because the BMI relied on predictions of velocity. The linear model predicted velocity in real time, and the velocity of the robot was set to the output of the model. In the present analysis, we describe the relationship between neuronal firing rate and movement by examining the firing rate of each neuron as a function of velocity and of the time lag between the firing rate and instantaneous velocity measurement (IVM) (see Fig. 1 B). IVMs were sampled every 100 ms. Neuronal activity was binned in 100 ms nonoverlapping windows. Binned firing rates were then normalized by subtracting the mean and dividing by the SD. Lags in the range of 1 s before to 1 s after the IVM were examined. For each lag τ, neuronal firing rate n was fitted to velocity using the following multiple linear regression: 
(1)
where t is time, n(t + τ) is neuronal firing rate at time τ from velocity measurement (negative τ corresponds to activity preceding velocity measurement, and positive τ activity succeeding velocity measurement), Vx(t) and Vy(t) are x- and y-components of velocity, a(τ), b(τ), and c(τ) are regression coefficients, and ϵ(t, τ) is the residual error. We characterized the degree of correlation between the neuronal firing rate and the velocity as the velocity tuning index (VTI), which was calculated as the square root of R2 for the regression (Eq. 1). This measure describes correlation between neuronal firing rate and movement velocity (both direction and amplitude taken into account) and represents it in normalized form. The better the neuronal firing rate tracked movement velocity, the higher was the VTI. The VTI did not depend on the absolute value of firing rate (e.g., highly active versus less active neurons) because of normalization. Preferred direction (PD; i.e., movement direction for which the neuron fired the most) was calculated as PD(τ) = arctan b(τ)/a(τ).
Frequency distributions of velocity (Vx and Vy) for different operations for monkey 1 (A-D) and monkey 2 (E-H). A, E, Pole control, hand velocity. B, F, Brain control with hand movements, hand velocity. C, G, Brain control with hand movements, robot velocity. D, H, Brain control without hand movements, robot velocity.
To examine whether the PDs of individual neurons in the population were preserved after the monkey started to perform in a new mode of operation (e.g., pole control vs BCWH), we used a correspondence index defined as follows: 
(2)
where statistical values are calculated for the neuronal samples of the PDs (measured in degrees) before and after the transition, α and β, respectively, the bar denotes the mean, and 90° is the mathematical expectation of 
in case they are uncorrelated. C was scaled from 0 to 1. Values of C approaching zero indicate no correspondence between the PDs, whereas values of C approaching unity indicate an exact correspondence.
To examine the contribution of correlation between neurons and tuning parameters, we conducted a shuffle test that destroyed the correlation between neurons by shifting spike trains of different neurons with respect to each other by a random interval in the range of 0 to 200 s. This procedure preserved the sequence of spikes of each neuron but removed correlations between the neurons. After shuffling, robot movements were recalculated using the same model with the fixed weights that had been used for the real-time brain control. The VTI of each individual neuron for shuffled data was a measure of tuning that would be expected simply because that neuron had a weighted contribution to the robot control signal. Higher values of VTIs for unshuffled data indicated that correlated firing between neurons improved tuning characteristics of individual neurons (i.e., neuronal sub-groups acted as “teams” producing movements in particular directions).
A separate analysis tested velocity tuning properties for instances in which the monkeys' hands or the robotic actuator moved in straight-line trajectories. Epochs of straight-line movements of at least 300 ms duration were selected. Firing rates were measured using a 300 ms window placed at different time lags relative to the selected epochs. Regression of Equation 1 was used to determine the VTI. In addition, directional tuning depth was calculated (Carmena et al., 2003). Directional tuning depth was equal to the difference in the average firing rate between the directions of maximum and minimum firing, divided by the SD of the firing rate.
Off-line predictions of hand velocity. Predictions of velocity were made from a weighted linear combination of neuronal firing rates using the following: 
(3)
where n(t + τ) is a vector of neuronal firing rates, at time t and time lag τ (negative lags correspond to past events), Vx is the x-component of velocity at time t, w(τ) is a vector of weights for each neuron at time lag τ, b is the y-intercept, and ϵ(t) is the residual error. An analogous equation was used to estimate Vy. These equations were solved using linear least-squares regression. Ten time lags preceding velocity measurement with a bin size of 100 ms were used for predictions. In a separate analysis, estimations of velocity were calculated using a single 100 ms bin placed at different lags relative to velocity. Estimation quality was quantified by the correlation coefficient (R).
Random neuron dropping was implemented as follows: 10 min of neuronal population data were used to fit the model and find the weights, w(τ). These weights were then used to make predictions of velocity on a different 10 min period. A single neuron was then randomly removed from the population, the model retrained, and new predictions generated on the same data. This process was repeated until no neurons remained. Average R as a function of number of neurons was obtained by repeating this procedure 100 times.
Results
Hand and robot velocity during different modes of operation
The linear model was trained during pole control. However, after the mode of operation was switched to brain control, robot movements did not precisely match movements of the hand. This is illustrated by Figure 1C-H, which shows representative traces of the hand and the robot during pole and brain control. During pole control, movements of the hand and the robot were in close correspondence (Fig. 1C), which can also be seen in time-dependent plots of the x- and y-components of position (Fig. 1E1,E2) and velocity (Fig. 1G1,G2). During BCWH, the hand and the robot moved differently (Fig. 1D), albeit in a correlated manner (Fig. 1F,H).
Figure 2 shows the frequency distribution of hand and robot velocities (Vx and Vy) during different modes of operation in two monkeys. The correlation coefficient between the robot and hand velocities was 0.995 ± 0.001 (mean ± SE) during pole control versus 0.61 ± 0.02 during BCWH in monkey 1 and 0.980 ± 0.001 versus 0.50 ± 0.04 in monkey 2. The difference between these values was highly significant statistically (p < 0.001; Wilcoxon rank-sum test). Ranges of hand velocities were similar during pole (Fig. 2A,E) and brain (Fig. 2B,F) control. However, a comparison of robot velocities (Fig. 2C,G) with the hand velocities that accompanied them during brain control (Fig. 2B,F) showed that lower robot velocities were less frequent than those of the hand. Absolute velocity of the hand (5.65 ± 0.03 cm/s in monkey 1 and 4.51 ± 0.04 cm/s in monkey 2) was statistically different (p < 0.001; Wilcoxon rank-sum test) from absolute velocity of the robot (6.93 ± 0.03 and 7.00 ± 0.03 for monkeys 1 and 2, respectively). Thus, the monkey's hand often moved slowly, but the robot moved faster, driven by neuronal ensemble activity. Ranges of robot velocities were similar during BCWH (Fig. 2C,G) and BCWOH (Fig. 2D,H).
Velocity tuning in M1 neuron tuned during pole control and brain control. A, Color plots of the firing rate of the neuron (color coded; key, bottom left) as a function of Vx and Vy (key, bottom left) for different lags with respect to IVM, different modes of operation (pole control and brain control with and without hand movements), and different velocity parameters (hand or robot movements). B, Velocity tuning index as a function of lag for different types of operation (color coded; key on top). C, Preferred direction as a function of lag (color coded; key on top). deg, Degrees.
Neuronal tuning during different modes of operation
We analyzed neuronal data collected in nine daily recording sessions in monkey 1 and in 20 sessions in monkey 2. Analyzed neuronal samples on average consisted of 177 units per recording session in monkey 1 and 72 units in monkey 2. Examination of neuronal rates at different time lags relative to IVM (Fig. 1B) showed modulations correlated with velocity magnitude and direction (Figs. 3, 4, 5). We refer to correlation between neuronal rate and velocity (Georgopoulos et al., 1988; Ashe and Georgopoulos, 1994; Moran and Schwartz, 1999; Carmena et al., 2003) as velocity tuning. The degree of velocity tuning was quantified as the VTI (see Materials and Methods). This normalized index measured the dependency of neuronal firing rate on both direction and magnitude of velocity. Because of normalization, the VTI was not biased toward highly active neurons. Neurons tuned before the IVM are more likely to have a causal relationship to movements, for example by projecting to spinal motoneurons. Neurons tuned after the IVM may reflect sensory feedback. Pre-IVM and post-IVM tuning can also be related to movement patterns, such as acceleration and deceleration (Sergio and Kalaska, 1998; Johnson and Ebner, 2000; Todorov, 2000). The dependency of firing rate on x- and y-components of velocity is depicted in color plots (Fig. 3A) (Johnson and Ebner, 2000; Lebedev et al., 2003; Paninski et al., 2004). The color encodes firing rate (red, high rate; blue, low rate). The center of each color plot represents the firing rate for zero velocity; x- and y-axes correspond to x- and y-components of velocity, respectively. These color plots are useful to evaluate the velocity tuning of a neuron at a glance. They contain information about both directional tuning of the neurons (direction is represented by a vector originating in the center of the plot) and their sensitivity to velocity magnitude (velocity magnitude is measured as the distance from the center of the plot). The color plots were constructed for different lags relative to the IVM and for different modes of operation (compare Figs. 3, 4, 5). Two sets of color-plots were calculated for BCWH, because during this mode of operation, hand movements did not match robot movements. One set relates firing rate to velocity of hand movements (Fig. 3A, second row from top). The other was calculated relative to robot velocity (Fig. 3A, third row from top). The color plots for BCWOH were calculated relative to robot velocity (Fig. 3A, bottom row). In addition, VTIs (Fig. 3B) and PDs (Fig. 3C) were plotted as a function of lag for different modes of operation.
An M1 neuron that was tuned to movement velocity only if the monkey's hand moved. Conventions are as in Figure 3. deg, Degrees.
Individual neurons exhibited a diversity of tuning patterns that depended on the mode of operation. Figures 3, 4, 5 illustrate three representative neurons recorded in M1. Their tuning characteristics were averaged over five consecutive recording sessions. The neuron shown in Figure 3 was tuned during both pole and brain control. The VTI of this neuron peaked at ∼200 ms before the IVM (Fig. 3B), and for all modes of operation, it was better tuned before the IVM than after it. Peak VTI values (Fig. 3B) for different modes of operation were significantly different (p < 0.001; Kruskal-Wallis ANOVA; α = 0.05 for Tukey's multiple comparison). The highest VTI values were observed during pole control (Fig. 3B, blue curve). During BCWH, tuning relative to hand movements (Fig. 3B, red curve) was significantly diminished. However, tuning relative to robot movements (Fig. 3B, black curve) was comparable with tuning during pole control (no statistical difference). Tuning during BCWOH was measured relative to robot movements. It was significantly lower than tuning relative to robot movements during BCWH but higher than tuning relative to hand movements during BCWH. The PD of this neuron was stable and similar for different modes of operation for the lags preceding the IVM (Fig. 3C). For the points after the IVM, the PD rotated for all modes of operation except BCWOH. Thus, after transition to brain control, this neuron retained many features of its original velocity tuning. However, it became less tuned to hand movements.
Figure 4 illustrates a neuron whose tuning mostly reflected hand movements. Changes in the VTI during different modes of operation were statistically significant (p < 0.001; Kruskal-Wallis ANOVA; α = 0.05 for Tukey's multiple comparison). This neuron was tuned during pole control and BCWH. However, velocity tuning vanished during BCWOH (Fig. 4A,B). Similar to the neuron illustrated in Figure 3, the highest VTIs of this neuron were observed for lags preceding the IVM, which may have been related to the causal relationship between the neuronal activity and hand movements. However, this neuron was also tuned after the IVM, especially during pole control. The PD depended on the lag relative to the IVM. It rotated counterclockwise (Fig. 4C) and completely reversed for lags succeeding the IVM. Comparable PD reversals have been reported previously in arm movement tasks and attributed to a representation of accelerative and decelerative forces (Sergio and Kalaska, 1998; Todorov, 2000). Note, in addition, that post-IVM tuning was markedly decreased when tuning was assessed relative to robot movements during brain control. Prevalence of pre-IVM tuning was also observed for the whole neuronal ensemble (see below). Because of the strong dependency of tuning on the presence of hand movements, this neuron is likely to have been critically involved in generation of descending motor commands (e.g., pyramidal tract neuron).
Pairwise comparison of pole control with different types of brain control (which was somewhat different from grouping the data from different sessions together, because different types of brain control were not tested in every session) showed that in the majority of neurons (77 ± 2%; mean ± SE), the peak VTI for hand movements decreased after the transition to BCWH. During BCWH, the peak VTI relative to robot velocity exceeded the peak VTI relative to hand velocity in 84 ± 2% of neurons. Each of these differences in proportions was highly significant statistically (p < 0.001; Wilcoxon signed rank test). Importantly, the peak VTI of a substantial number of neurons (38 ± 2%) was higher during BCWOH than during pole control. A neuron of the latter type (Fig. 5) had a peculiar tuning pattern: it was much better tuned during BCWOH than other conditions (p < 0.001; Kruskal-Wallis ANOVA; α = 0.05 for Tukey's multiple comparison).
Population averages for all recording sessions confirmed that neurons typically became less tuned to hand movements and better tuned to robot movements (Table 1). Population averages of peak VTIs were 0.074 ± 0.002 for pole control, 0.056 ± 0.002 for BCWH relative to hand, 0.077 ± 0.003 for BCWH relative to robot, and 0.070 ± 0.001 for BCWOH relative to robot (mean ± SE) in monkey 1 and 0.104 ± 0.005, 0.069 ± 0.004, 0.118 ± 0.002, and 0.126 ± 0.007 in monkey 2. Changes in the peak VTI were statistically significant (p < 0.001; Kruskal-Wallis ANOVA). Pairwise comparison (Tukey's multiple comparison; α = 0.05) showed that the peak VTI during BCWH relative to hand was significantly less than the peak VTI during pole control and BCWH and BCWOH relative to robot. This decrease in tuning to hand movements and better tuning to the actuator instead was observed in the majority of sampled areas. Statistical comparison of peak VTIs for individual cortical areas is presented in Table 1.
Velocity tuning index (mean ± SE) by cortical area and condition
An M1 neuron with enhanced tuning to robot velocity during brain control without hand movements. Conventions are as in Figure 3. deg, Degrees.
Tuning patterns for the whole ensemble
Color plots of Figure 6 show VTIs as a function of lag for the whole neuronal ensemble. They present VTI curves (Figs. 3B, 4B, 5B) in color-coded format, stacked one above the other, so that each line corresponds to a neuron. It can be seen that for the ensemble as a whole, tuning to hand movements decreased during BCWH (Fig. 6A2) compared with tuning during pole control (Fig. 6A1) and tuning to robot movements during BCWH (Figs. 6A3, 4). Individual neurons exhibited a diversity of VTI changes. Averaging VTIs for the whole ensemble (Fig. 6B) confirmed a decrease in tuning to hand movements during BCWH (Fig. 6B, compare blue and red curves). VTIs for robot movements (Fig. 6B, black and green curves) were clearly higher before than after the IVM. This occurred for two reasons. First, only the bins preceding the IVM contributed to the on-line predictions of velocity. Second, as a result of the delay of the robotic system (∼100 ms), the curves were shifted to the left. VTIs for hand movements during pole and brain control peaked closer to IVM compared with robot movements. Average VTIs for the whole ensemble were similar to those for M1 (Fig. 6B,C), because in the sampled neuronal population, M1 neurons were tuned most strongly.
The occurrence of tuning of the neuronal ensemble to robot velocity during BCWOH is not surprising, because the weighted sum of neuronal firing rates was made to control the robot. However, neuronal tuning to robot movements during BCWOH was not a simple consequence of using the linear model (as it would be if just one neuron controlled the BMI). Rather, correlated firing of the whole ensemble played a role in the tuning properties of individual neurons. This was shown using the shuffle test. The average VTI decreased after the correlation between the neurons was destroyed by the shuffling procedure (see Materials and Methods), as seen in Figure 6, B and C (dotted green lines vs solid green lines) and Table 1. This decrease was statistically significant for both monkeys (p < 0.001; Wilcoxon signed rank test; means and SDs are presented in Table 1). Thus, tuning of any individual neuron was not related only to the weights assigned to that neuron by the linear model but also depended on the firing of the rest of the ensemble. We have shown previously (Carmena et al., 2003) that firing of individual neurons was correlated, and this correlation increased in brain control. Because of the enhancement in correlated activity, VTIs of individual neurons increased as well.
Preferred directions plotted as a function of lag relative to the IVM (Fig. 7A) showed that, generally, there was no fixed PD for each neuron. Rather, PDs rotated with lag changes similarly to the examples in Figure 3 and 4. This rotation is consistent with previous results of Sergio and Kalaska (1998), Todorov (2000), and Johnson and Ebner (2000). After the mode of operation was switched to brain control, if the monkeys continued to move their arms, the pattern of PDs generally resembled that of pole control (Fig. 7A1-A3). Correspondence index analysis indicated that the best correspondence between the PDs measured with respect to hand movements during pole control and BCWH occurred at the time of the IVM (Fig. 7B, red line), whereas PDs measured with respect to robot movements were in best correspondence to pole control ∼100 ms before the IVM (Fig. 7B, black line). This is in agreement with average VTI curves (Fig. 6B). Statistical analysis of peak values of correspondence indices showed better correspondence between the PDs during pole control and BCWH measured with respect to robot movements compared with the values measured with respect to hand movements, which was statistically significant for monkey 1 (0.38 ± 0.01 vs 0.33 ± 0.01; mean ± SE; p < 0.05; Wilcoxon signed rank test) and did not reach statistical significance for monkey 2 (0.66 ± 0.03 vs 0.59 ± 0.04; p > 0.05). The PDs during BCWOH (Fig. 7A4) were confined to a narrower angular range and rotated much less with changes in the lag than those observed during pole control (Fig. 7A1) and BCWH (Fig. 7A2,A3) (Carmena et al., 2003). Correspondence index analysis, however, indicated a relatively good match between the PDs during BCWOH and pole control for the 100 ms lag preceding the IVM (Fig. 7B, green line). Peak values of the correspondence between the PDs during BCWOH and pole control for monkey 1 (0.24 ± 0.01) were significantly less than the correspondence index between pole control and BCWH (p < 0.05; Wilcoxon signed rank test), and the values for monkey 2 (0.52 ± 0.05) were not significantly different (p > 0.05).
Neuronal tuning for selected hand movements
Changes in neuronal tuning to hand movements that occurred after the monkeys started to perform in brain control could be related, among other factors, to differences in the pattern of hand movements. To test whether these changes would still be present if the same hand movements were compared, we selected epochs during which the monkeys' hands or the robot moved in straight lines for at least 300 ms (Fig. 8A). For these epochs, we constructed classical tuning curves, which represented average firing rate as a function of movement direction (Fig. 8B). Directional tuning depth (Fig. 8C) and VTI (Fig. 8D), calculated as functions of the lag between the 300 ms window in which firing rate was measured, showed the same effect as demonstrated with other methods: both directional tuning depth and VTI with respect to hand movements were less during BCWH than during pole control. Average peak tuning depth was 0.55 ± 0.01 and 0.71 ± 0.01 in monkeys 1 and 2, respectively, during pole control, and it was 0.44 ± 0.01 and 0.55 ± 0.01 during brain control. Average peak VTI was 0.203 ± 0.005 and 0.277 ± 0.005 (monkeys 1 and 2) during pole control versus 0.161 ± 0.004 and 0.207 ± 0.005 during brain control. These differences between tuning to hand movements in pole control and brain control were statistically significant (p < 0.001; Wilcoxon signed rank test).
Off-line predictions of hand velocity during pole and brain control
Similar to the decrease in VTIs, the quality of off-line predictions of hand velocity decreased during BCWH (Fig. 9A). In this analysis, the linear model was trained to predict hand velocity using either pole-control or brain-control data as training epochs. Ten minute (Fig. 9A) or 5 min (Table 2) data segments were used for training. Predictions were then tested on a new 10 or 5 min epoch (splitting the data into 5 min epochs was useful for statistical analysis, although prediction quality was less compared with 10 min epochs, which was also used for training the model during online predictions). Prediction quality dropped during brain control (Fig. 9A). This decrease occurred quite abruptly after the operation was switched to brain control (Fig. 9B). Average R values for predictions of hand velocity using 5 min segments of data during pole control were 0.64 ± 0.02 and 0.48 ± 0.01 for monkey 1 and 2, respectively, whereas they were 0.40 ± 0.01 and 0.30 ± 0.01 for predictions of hand velocity made during BCWH (Table 2). This change in R values was statistically significant (p < 0.001; Kruskal-Wallis ANOVA; α = 0.05 for Tukey's multiple comparison).
Prediction of hand velocity during pole control and brain control with hand movements using 5 min epochs for training the model and testing predictions
To inspect changes in contribution to velocity prediction for different lags between the velocity and neuronal activity, velocity predictions were performed using a single 100 ms sampling window. The window was placed at different lags with respect to the IVM. This analysis showed a decrease in prediction for all tested lags (Fig. 9C).
It is important to note that the same drop in prediction quality for hand velocity as we show here in the off-line analyses occurred during the experimental sessions as well after the mode of operation was switched to brain control. However, such deterioration in prediction quality only meant that the model became less predictive of the velocity of the hand. This decrease in predictions of hand movement did not impede the behavioral performance, which was determined by the robot movements rather than hand movements. Behavioral performance improved with training in brain control (Carmena et al., 2003), suggesting that cortical representation of the robot was optimized at the expense of representation of the animal's own limb.
Changes in velocity tuning for the whole ensemble. A, VTI as a function of lag relative to IVM for the ensemble recorded in monkey 1. Each row corresponds to a neuron. VTI values are shown for pole control (1), brain control with hand movements, relative to hand (2) and to robot (3), and brain control without hand movements, relative to robot (4). ips, Ipsilateral. B, Average VTIs as a function of lag for different conditions (color coded; key on right). Average VTI for shuffled data during brain control without hand movements is shown by the green dotted line. Averages are shown for the whole ensemble. C, Average VTIs as a function of lag for M1 only.
Discussion
The principal finding of this study is that once cortical ensemble activity is switched to represent the movements of the artificial actuator, it is less representative of the movements of the animal's own limb. Although the original rate of the sample of recorded cortical neurons was initially modulated in association with hand movements, after these neurons became involved in brain control, they became less correlated with hand movements. Moreover, neuronal tuning to the movements of the actuator was enhanced because of increased correlation between the neurons. This result accords with the vast literature showing that neuronal modulations in cortical motor areas do not exclusively reflect limb movements and that cognitive signals are abundant in these areas (for review, see Wise et al., 1997; Georgopoulos, 2000; Johnson et al., 2001; Andersen and Buneo, 2002; Paz et al., 2004). For instance, neuroimaging studies in human subjects have shown an activation of cortical motor areas in the absence of limb movements (e.g., when subjects imagine that they produce movements) (Leonardo et al., 1995; Porro et al., 1996; Ehrsson et al., 2003) or when the illusion of limb movement is induced by muscle vibration (Naito et al., 1999). Moreover, Graziano et al. (Graziano, 1999; Graziano et al., 2000) demonstrated that if the vision of the monkey's own arm was occluded, PMd and PP neurons reflected the position of a visible, realistic false arm, but at the same time, they were not sensitive to unrealistic substitutes for the arm. A recent neuroimaging study confirmed this result by demonstrating similar activations that occurred in human premotor cortex when the subjects experienced a rubber hand illusion (Ehrsson et al., 2004). Thus, cortical representation of the limb is highly flexible and susceptible to illusions and mental imagery.
Changes in preferred directions for the whole ensemble. A, PD as a function of lag relative to IVM for the ensemble recorded in monkey 1. Each row corresponds to a neuron. PDs are shown for pole control (1), brain control with hand movements, relative to hand (2) and to robot (3), and brain control without hand movements, relative to robot (4). ips, Ipsilateral. B, Correspondence index describing the similarity between PD distribution during pole control and brain control. deg, Degrees.
Neurophysiological studies have shown that modulations in cortical areas that occur during movement planning and execution often reflect abstract parameters of movements rather than low-level details. Cortical neurons can encode movement direction regardless of muscle patterns (Kakei et al., 1999, 2001), represent target of movement and movement of the hand-controlled visual marker toward it rather than actual hand trajectory (Alexander and Crutcher, 1990; Shen and Alexander, 1997a,b; Ochiai et al., 2002), encode multiple spatial variables (Lebedev et al., 2004; Cisek and Kalaska, 2005), reflect orientation of selective spatial attention (Lebedev and Wise, 2001), and represent misperceived movements of visual targets (Lebedev et al., 2001; Lebedev and Wise, 2002; Schwartz et al., 2004). In addition, frontal cortex neurons can flexibly adapt to new visuomotor (Wise et al., 1998; Paz et al., 2003; Paz and Vaadia, 2004) and dynamic (Li et al., 2001; Padoa-Schioppa et al., 2004) conditions. Based on the results presented here, we suggest that neuronal mechanisms underlying these adaptive properties can also subserve cortical ensemble adaptation during the BMI operation.
Operating a BMI is not unlike using mechanical tools, and the idea that body schema can extend along a wielded tool was formulated almost 100 years ago (Head and Holmes, 1911). This suggestion was substantiated much later by neurophysiological data showing that cortical neurons responsive to both somatosensory and visual stimuli extended their receptive fields along a rake that a monkey used to retrieve distant objects (Iriki et al., 1996). Effects consistent with the idea of remapping body schema after tool usage have been also observed in human psychophysics experiments (for review, see Maravita et al., 2003). A recent neuroimaging study (Maruishi et al., 2004) reported specific activations of the right ventral premotor cortex during manipulation of a myoelectric prosthetic hand. The authors interpreted this result as the recognition by the brain of the prosthetic hand as a “high performance alternative to a real hand.” Perhaps a long-term operation of a BMI-controlled actuator would lead to even more vivid perception of the actuator being a body part rather than merely a tool or an alternative to the real hand. In support of this suggestion, perceived voluntary movements of phantom limbs in amputees were associated with primary sensorimotor cortex activation (Roux et al., 2003).
When interpreting neuronal patterns during BMI control as reflecting a new representation of the artificial actuator, however, one should be careful, because the BMI design effectively makes any modulation of neuronal activity translate into the movements of the actuator. One piece of evidence that the system adapts to control the actuator comes from the nonrandomness of the actuator movements and behavioral improvements as the monkey trains (Taylor et al., 2002; Carmena et al., 2003). The second piece of evidence is coordinated activity of the neurons in the ensemble. Extending our previous demonstration of increased correlation between the neurons during brain control (Carmena et al., 2003), we showed here using the shuffle test that neuronal correlation enhanced tuning to robot velocity in individual neurons. Moreover, tuning of the same neurons to hand velocity decreased, and this decrease occurred even for the relatively stereotypical subset of hand movements having a straight-line trajectory (Fig. 8). Together, these observations suggest that the neuronal ensemble adapted to represent the actuator at the expense of the representation of the animal's own limb.
Analysis of neuronal tuning for straight trajectories of the hand. A, Straight trajectories of the hand selected for eight movement directions during pole control and brain control with hand movements. B, Directional tuning curves calculated for the traces above by averaging firing rates for different movement directions within a 300 ms window leading the movement by 100 ms. Each horizontal line represents a curve for a particular neuron. ips, Ipsilateral; deg, degrees. C, Directional tuning depth as a function of lag between the firing rate window and movement trace. D, Velocity tuning index.
Tuning patterns in the entire ensemble (Fig. 5A) during pole control and brain control with and without hand movements showed that, in general, most of the neurons that encoded movement velocity during pole control continued to do so during the BMI operation. How then could the monkeys operate the BMI and at the same time inhibit movements of their hands? We suggest that subtle rather than dramatic differences in the way the activity of individual neurons was combined at the ensemble level determined this novel behavioral outcome. For example, although the overall modulation pattern was similar across behavioral states, some neurons (Fig. 4) did not modulate when the monkeys' hands stopped moving, whereas others (Fig. 3) were strongly modulated. It should be noted that the neuronal ensemble that was recorded and used to control the BMI constituted only a tiny part of the entire network normally involved in controlling arm movements. It is possible that these selected neurons were preferentially modulated to control the actuator during BCWOH, whereas the ones that were not recorded from continued to operate as normally. Alternatively, the ensemble controlling the BMI could be assisted by a neuronal network to which it was connected. Additional studies will be needed to resolve between these possibilities.
Our finding that after learning to control an artificial actuator using a BMI, neuronal ensembles represented the actuator velocity better than the velocity of the animal's own arm appears to support the theory of optimal feedback control (Todorov and Jordan, 2002; Scott, 2004). This theory describes the motor system as a stochastic feedback controller that optimizes only those motor parameters that are necessary to achieve task goals. In our study, the task goals were achieved by the robot. Thus, in terms of optimal feedback control theory, cortical ensembles should adapt to represent the goal-related variable although it is related to the robot instead of the limb.
The present observations also suggest that the quality of predictions in open-loop mode that computational neuroscientists are striving to achieve may be relatively unimportant for closed-loop BMI control. Once the neuronal ensemble starts to control a BMI, imperfections in the model may be compensated by neuronal adaptations that improve control of the actuator movements. Whereas they improve the BMI performance, these adaptations may actually worsen the predictions of hand movements, the effect that we observed here. These considerations concur with the results of Taylor et al. (2002), who showed differences in tuning properties of M1 neurons during hand movements and during a BMI control that used a coadaptive algorithm for adjusting the model weights in the absence of hand movements.
To clarify the dynamics of cortical plasticity associated with BMI control and to elucidate the exact mechanism of cortical adaptation to represent an artificial actuator will require additional studies. It is very likely that the features of such adaptations will depend on concrete requirements for the BMI operations. For example, it is possible that if the subject is required to subserve the BMI control with limb movements, the system may come to an optimal solution in which both the limb and the actuator are represented. Alternatively, if limb movement is not required, the system may gradually evolve to a different state (Taylor et al., 2002; Carmena et al., 2003). It is also possible that certain neuronal adaptations could work as temporary solutions to maintain performance. Previously, we reported (Carmena et al., 2003) that, as a consequence of the increase in correlation between the neurons during BCWOH, PDs of many neurons became more similar (Fig. 6A). This effect may appear maladaptive, because PD diversity is likely to improve directional encoding (Georgopoulos et al., 1986, 1988). However, our present analysis using the shuffling procedure showed that tuning of individual neurons was increased because of the correlation between neurons. Thus, the neuronal ensemble controlling the BMI may have been optimizing its performance by finding a trade-off between the magnitude of tuning and the diversity of PDs. It is possible that increased correlation was a temporary effect, and with prolonged and intensive training to control the BMI, the neuronal ensemble would eventually return to the original distribution of PDs and neuronal correlations. Even during learning of relatively simple motor skills, highly distributed cortical networks undergo transient changes before settling to the level of automated performance, as shown by human neuroimaging studies (Jueptner et al., 1997; Floyer-Lea and Matthews, 2004). Interestingly, epochs of high attention to motor performance (and increased attention is very likely needed for successful learning to control a BMI) have been found to be associated with increases of synchrony between neurons (Murthy and Fetz, 1992, 1996a,b; Riehle et al., 1997, 2000; Hatsopoulos et al., 1998; Lebedev and Wise, 2000; Baker et al., 2001).
Off-line predictions of hand velocity under different conditions (key on top). A, Neuron-dropping curves depicting hand prediction quality as a function of neuronal sample size. B, Temporal dynamics of prediction quality for the entire ensemble. C, Prediction quality estimated using a single 100 ms bin placed at different lags with respect to IVM.
In conclusion, we propose that a wide-spread process of cortical plasticity, resulting from BMI operation, may lead to the establishment of a representation of an artificial actuator into multiple cortical areas. This raises the interesting hypothesis that the distinct levels of primate proficiency in tool manufacturing and utilization may be defined by the ability of these animals to incorporate the properties of these tools as extensions of the subject's own body.
Footnotes
This work was supported by grants from the Defense Advanced Research Projects Agency, the James S. McDonnel Foundation, the National Institutes of Health, the National Science Foundation (M.A.L.N.) and the Christopher Reeve Paralysis Foundation (J.M.C.). We thank Dragan F. Dimitrov for conducting the neurosurgery, Gary Lehew for manufacturing multielectrode arrays and engineering the experimental set-up, Roy E. Crist and Ming Qian for programming experimental software, Parag G. Patil, Erin Phelps, Laura M. Oliveira, Kristen Shanklin Dewey, and Aaron J. Sandler for surgical assistance and animal care, Sung-Phil Kim for evaluating additional adaptive filters, and Susan A. Halkiotis and Nathan Fitzsimmons for manuscript preparation.
Correspondence should be addressed to Dr. Miguel Nicolelis, Department of Neurobiology, Duke University Medical Center, Bryan Research Building, Box 3209, Durham, NC 27710. E-mail: nicoleli{at}neuro.duke.edu.
Copyright © 2005 Society for Neuroscience 0270-6474/05/254681-13$15.00/0
