Abstract
The cerebellum has been implicated in processing motor errors required for on-line control of movement and motor learning. The dominant view is that Purkinje cell complex spike discharge signals motor errors. This study investigated whether errors are encoded in the simple spike discharge of Purkinje cells in monkeys trained to manually track a pseudorandomly moving target. Four task error signals were evaluated based on cursor movement relative to target movement. Linear regression analyses based on firing residuals ensured that the modulation with a specific error parameter was independent of the other error parameters and kinematics. The results demonstrate that simple spike firing in lobules IV–VI is significantly correlated with position, distance, and directional errors. Independent of the error signals, the same Purkinje cells encode kinematics. The strongest error modulation occurs at feedback timing. However, in 72% of cells at least one of the R2 temporal profiles resulting from regressing firing with individual errors exhibit two peak R2 values. For these bimodal profiles, the first peak is at a negative τ (lead) and a second peak at a positive τ (lag), implying that Purkinje cells encode both prediction and feedback about an error. For the majority of the bimodal profiles, the signs of the regression coefficients or preferred directions reverse at the times of the peaks. The sign reversal results in opposing simple spike modulation for the predictive and feedback components. Dual error representations may provide the signals needed to generate sensory prediction errors used to update a forward internal model.
Introduction
Detecting and correcting performance errors is essential to fine-tune motor behavior, including on-line control and motor adaptation (Todorov and Jordan, 2002; Berniker and Kording, 2008; Shadmehr et al., 2010). Early concepts focused on closed loop control in which the ongoing motor commands are updated by sensory feedback (Wolpert and Miall, 1996). However, closed loop control using delayed feedback can be inadequate and even unstable (Wolpert and Miall, 1996; Bhushan and Shadmehr, 1999; Kawato, 1999; Wolpert and Ghahramani, 2000). Forward internal models provide a solution to the problem of long feedback delays by predicting the sensory consequences of a motor command (Robinson, 1975; Miall et al., 1993; Wolpert et al., 1995; Wolpert and Ghahramani, 2000; Tseng et al., 2007; Shadmehr et al., 2010). The difference between the prediction and the sensory feedback, a sensory prediction error, can be used to control movements and facilitate motor adaptation (Wolpert and Miall, 1996; Wolpert and Ghahramani, 2000; Shadmehr et al., 2010).
A long-standing hypothesis is that the cerebellum detects and corrects for errors (Oscarsson, 1980). Error processing is a central component of many hypotheses of cerebellar function, including motor learning (Marr, 1969; Gilbert and Thach, 1977; Ito, 2002) and providing internal models (Wolpert et al., 1998; Kawato, 1999; Shadmehr and Krakauer, 2008; Shadmehr et al., 2010). It has been hypothesized that the cerebellum computes sensory prediction errors (Wolpert and Miall, 1996; Tseng et al., 2007; Shadmehr et al., 2010). Error-related activation occurs in the cerebellum during the learning of visuomotor transformations or force fields (Imamizu et al., 2000; Diedrichsen et al., 2005; Ide and Li, 2011). Cerebellar damage results in deficits in adapting limb and eye movements, and impairs a patient's ability to learn from movement errors (Martin et al., 1996; Maschke et al., 2004; Smith and Shadmehr, 2005; Morton and Bastian, 2006; Xu-Wilson et al., 2009; Criscimagna-Hemminger et al., 2010). Processing errors appears to be integral to cerebellar function.
Error signals in the cerebellum are generally thought to be transmitted by climbing fiber input. Complex spike discharge is associated with errors during eye movements (Barmack and Simpson, 1980; Graf et al., 1988; Kobayashi et al., 1998; Soetedjo and Fuchs, 2006; Medina and Lisberger, 2008). During arm movements, complex spikes modulate with perturbations (Gilbert and Thach, 1977; Wang et al., 1987), adaptation to visuomotor transformations (Ojakangas and Ebner, 1994), and end point errors (Kitazawa et al., 1998). Whether Purkinje cell simple spike discharge encodes errors has received less attention. However, instructive signals in the simple spike firing contribute to cerebellar-dependent learning in the vestibulo-ocular reflex (Ke et al., 2009). Learning-related plasticity in the vestibular/cerebellar nuclei also suggests that simple spike firing encodes errors (Lisberger et al., 1994; De Zeeuw and Yeo, 2005; Ohyama et al., 2006). During arm movements, simple spike activity modulates with trial success or failure (Greger and Norris, 2005) and correlates with several error measures (Roitman et al., 2009). However, in the latter study, the error parameters evaluated were not independent of arm kinematics. The present study examines performance error encoding in the simple spike firing of Purkinje cells in lobules IV–VI. We report that simple spike discharge signals both prediction and feedback about errors that are independent of the kinematics encoded.
Materials and Methods
All animal experimentation was approved by the Institutional Animal Care and Use Committee of the University of Minnesota and conducted in accordance with the guidelines of the National Institutes of Health.
Random tracking.
This study used a previously described random tracking task (Hewitt et al., 2011) and, therefore, the paradigm is only briefly detailed here. Two rhesus monkeys (Macaca mulatta; female 6.3 kg, male 6.2 kg) were trained to track a randomly moving target (3 cm diameter) on a vertical screen (12 × 12 cm, positioned 50 cm in front of the animal) using a robot manipulandum (In Motion 2; Interactive Motion) with 2 degrees of freedom in the horizontal plane. Visual feedback was provided by a “+”-shaped cursor (0.5 × 0.5 cm). During the paradigm, the monkey needed to maintain the cursor within the target. Due to the task difficulty, however, brief excursions (<500 ms) outside the target were allowed. The paradigm started with an initial hold inside a stationary target for a random period of time (1000–2000 ms). The initial target position on the screen was also random. Next, the target moved along a 6–10 s trajectory selected randomly from 100 trajectories defined a priori. Generated using the sum of random sine waves (Imamizu et al., 2000; Paninski et al., 2004), the trajectories were lowpass filtered and selected to avoid sharp turns and large changes in speed. Target speed was controlled by the two-thirds power law, speed = cK−1/3, where K is the path curvature and c is a constant chosen so that the average speed was ∼4 cm/s (Viviani and Terzuolo, 1982; Lacquaniti et al., 1983). The trajectory ended with a final hold period of 1000 ms during which the target was stationary.
Surgical procedures, electrophysiological recordings, and data collection.
Head-restraint hardware and a recording chamber targeting lobules IV–VI of the intermediate and lateral cerebellar zones were chronically implanted over the ipsilateral parietal cortex in each animal using aseptic techniques and full surgical anesthesia (Hewitt et al., 2011). The positions of the electrodes were confirmed by radiographic imaging techniques that combined a CT scan of the skull with magnetic resonance imaging of the cerebellum (Hewitt et al., 2011).
Purkinje cells were recorded extracellularly after the animal's full recovery, using quartz-platinum/tungsten single microelectrodes (1–2 MΩ impedance; Thomas Recording). When searching for a cell, the animal was sitting quietly and not performing any task. Purkinje cells were identified by the presence of complex spikes (Thach, 1968; Ojakangas and Ebner, 1992; Fu et al., 1997; Roitman et al., 2005; Hewitt et al., 2011). However, the recordings focused on optimizing the simple spike discrimination due to the long time period required to complete the paradigm. The spike trains and the raw electrophysiological data were digitized and stored to disk at 1 and 32 kHz, respectively. For analyses, the spike trains were transformed into an instantaneous firing rate with 20 ms binning using the fractional intervals method and then lowpass filtered (12th order Butterworth with 12 Hz cutoff). The fractional intervals method is a common technique used to determine the instantaneous firing rate in equal bins based on the inverse of the interspike intervals. The method accounts for incomplete (i.e., fractional) intervals that occur at the start and end of each bin (Schwartz, 1992; Taira et al., 1996; Reina et al., 2001). For display and analyses, the mean firing rate for each trial was subtracted from the instantaneous firing rate.
Target center and cursor positions were sampled and stored at 200 Hz. All other kinematic and error variables were derived from the position data (Roitman et al., 2005; Pasalar et al., 2006). Similar to the firing data, the kinematic and error parameters were also down sampled to 50 Hz and filtered (lowpass 12th order Butterworth with a 12 Hz cutoff). The analyses evaluating the relation among the simple spike firing and the behavioral variables were restricted to the tracking period.
Behavioral analysis: task errors and kinematic parameters.
Although several error terms were evaluated, four parameters were selected based on an initial assessment of the simple spike modulation and the need for independence between these error terms and hand kinematics (see below). The error parameters were based on the instantaneous relationships between the cursor (c), which mapped the hand movement onto the screen, and the target center (tg) (see Fig. 1A). The two position errors, XE and YE, were defined as the differences between cursor (Xc, Yc) and target (Xtg, Ytg) centers:
The radial error (RE) was defined as the distance between the cursor and target center:
We also defined the position direction error (PDE), that is the direction the hand needs to move to return to the target center. Formally, the PDE is the difference between the cursor movement direction and the direction of the position error vector:
in which VXc and VYc are the cursor velocity components obtained by numerical differentiation of the corresponding position data. A PDE of 180° shows that the hand is moving from its current position directly toward the target center and a PDE of 0° shows the hand is moving away from the target center.
The probability density function for each error was determined across all recording sessions for both monkeys. Probability density functions were determined by binning the area over which the error could occur and calculating the probability for each bin. For XE and YE, we used a 2D grid (6.0 × 6.0 cm2) centered on the target center (0, 0) that was partitioned into 0.15 × 0.15 cm2 bins. RE spanned from 0 to 3.2 cm to accommodate errors well beyond the 1.5 cm target radius. RE was partitioned into 0.05 cm bins that divided the target space into concentric rings of 0.05 cm width. Note that the areas of the rings increase as the distance from the center increases. PDE spanned from 0 to 360° and was divided into 5° bins.
Linear modeling of simple spike firing based on firing residuals.
The key aim of this study is to assess the simple spike firing modulation with the error parameters described above. Quantification of the simple spike firing involved several levels of linear regression models including: (1) examining each error parameter in isolation in a simple linear regression, (2) validating the independence of the model parameters, and (3) incorporating the error parameters into a single, multilinear model. However, our previous study using the same random tracking task demonstrated that the simple spike discharge is also highly modulated by movement kinematics (Hewitt et al., 2011). Specifically, simple spike firing modulates with hand position (P), velocity (V), and speed (S). We refer to this as the PVS model. Several other studies have documented similar encoding of kinematics during arm movements (Marple-Horvat and Stein, 1987; Fortier et al., 1989; Coltz et al., 1999; Roitman et al., 2005). Therefore, it is essential to demonstrate that the simple spike discharge encodes the error parameters independent of the kinematics. It is also essential to demonstrate that encoding of individual error parameters is mutually independent (e.g., modulation with XE is independent of YE, RE, and PDE). We approached this issue of independence by regressing the parameters of interest against the relevant firing residuals.
To evaluate single error parameters, firing residuals (FRs) were first obtained from a multilinear model that included the PVS model terms and the remaining error parameters. The PVS parameters are defined by the cursor movements and include: position (Xc and Yc), velocity (VXc and VYc) and speed ((VXc)2 + VYc)2)1/2. As an example, the firing residuals needed to evaluate XE independently of the other error terms and kinematics are obtained by regressing actual firing (F) to this multilinear model:
The purpose of this regression is to remove the variability in the firing related to the kinematics (PVS) and the complementary error parameters (YE, RE, and PDE). The error term, ε, is the firing residuals and represents the remaining variability in the simple spike discharge. All regressions used to obtain firing residuals (e.g., Eqs. 5–7) were computed using the instantaneous firing and movement data from individual trials (i.e., nonaveraged data). Firing residuals were determined at each of the 20 ms time shifts (τ) between the simple spike firing and the model predictors in a −500 to 500 ms window. We denote the firing residuals from Equation 5 as FRτPVS, YE, RE, PDE(t). In this notation, the subscript refers to the regression parameters used initially to obtain the firing residuals and the superscript refers to the τ value at which the firing residuals were computed. Similarly, we computed the firing residuals FRτPVS, XE, RE, PDE(t), FRτPVS, XE, YE, PDE(t), and FRτPVS, XE, YE, RE(t), needed for the analysis of YE, RE, and PDE, respectively.
We also determined the firing residuals, FRτPVS(t), in which the firing variability related to kinematics was removed (Eq. 6). FRτPVS(t) was obtained from regressing the firing to the PVS model:
Finally, we computed the firing residuals, FRτER(t), in which the firing variability related to the error parameters was removed by regressing the firing to the error parameters:
Analysis of simple spike firing modulation with individual task errors.
Next, using a simple linear regression, the appropriate firing residuals were regressed to the individual error parameters:
The regression of the firing residuals on individual error parameters (Eqs. 8–11) are denoted by: FRτPVS, YE, RE, PDE (XE(t)), FRτPVS, XE, RE, PDE (YE(t)), FRτPVS, XE, YE, PDE (RE(t)), and FRτPVS, XE, YE, RE (PDE(t)). The standard cosine model (Eq. 11) was used to assess directional tuning with PDE (Georgopoulos et al., 1982; Amirikian and Georgopoulos, 2000).
For the regressions fitting firing residuals (Eqs. 8–11 and 17–18 below) or actual firing (Eqs. 12–15 and 19 below) to the error models, the error space (XE, YE, RE, PDE) was partitioned in eight equal bins along each dimension ranging from −3 to 3 cm for XE and YE, 0 to 3.2 cm for RE, and 0 to 360° for PDE. This partitioning resulted in 4096 4D bins. Error parameters and neural data were sorted into these bins and averaged. Therefore, the final regression step used averaging to increase the signal-to-noise ratio as commonly used in cerebellar single unit studies. However, the 4096 bins represent a much finer partitioning than used in most studies using averaging (Marple-Horvat and Stein, 1987; Fortier et al., 1989; Gomi et al., 1998; Coltz et al., 1999; Medina and Lisberger, 2009; Roitman et al., 2009).
All linear model analyses used repeated fittings in which the time series of neural discharge (firing residuals or actual firing) was shifted relative to the model's time series of independent variable(s) (Ashe and Georgopoulos, 1994; Gomi et al., 1998; Medina and Lisberger, 2009; Roitman et al., 2009; Hewitt et al., 2011). The temporal shifts were used to assess the lead/lag (τ) between the neural signals and behavioral parameters (i.e., errors or kinematics). Temporal shifts ranged from −500 to 500 ms in 20 ms increments, with negative τ values representing that the neural signals lead the model regressors. Therefore, for each τ, a new regression based on the firing or firing residuals is computed. As in other studies (Ashe and Georgopoulos, 1994; Gomi et al., 1998; Medina and Lisberger, 2009; Hewitt et al., 2011), τ does not add an additional degree of freedom. From the time-shifted regressions, we generated temporal profiles for the coefficient of determination (R2) and the regression coefficients (βs) as functions of τ. The τ value at which an R2 profile achieved its maxima is referred to as the optimal τ. Significance was tested at each τ using the F statistic with the degrees of freedom determined by the number of observations and the number of predictors. The Type 1 error rate was set at α = 0.05. We imposed an additional threshold of R2 ≥ 0.02 for a significant fit.
The βs from Equations 8–10 were used as a measure of a Purkinje cell's firing sensitivity with XE, YE, and RE, respectively. To compare the sensitivity among different parameters and across cells, standardized regression coefficients (βstd) were also determined (Rencher, 2000). For PDE (Eq. 11) depth of modulation (DM) was defined as the amplitude of the best fit cosine and is given by DM = (βsPDE2 + βcPDE2)1/2. The preferred direction (i.e., the direction of maximal firing) was also determined using the coefficients βsPDE and βcPDE (Eq. 11) (Georgopoulos et al., 1982; Amirikian and Georgopoulos, 2000).
Validation of individual error models.
As noted above, a fundamental issue is whether the simple spike encoding of an error parameter is independent of the kinematics and the other error parameters. Therefore, additional regressions were performed. We first compared the results based on the firing residuals for each error parameter (Eqs. 8–11) with the regression results based on the actual firing (Eqs. 12–15):
To demonstrate that the firing residuals lack contributions from the kinematics and the complementary error terms, “control” regressions were computed for each error term. The control regressions are used to demonstrate the effectiveness of removing information about specific parameters because the firing residuals were based on the nonaveraged data and the final regressions involved averaged data. Two types of controls were determined. The first type is referred to as a PVS control and was based on the firing residuals for an error parameter in which the contributions of the kinematics and complementary error parameters were removed. These firing residuals were fit to the PVS model. For example, the PVS control regression for XE is modeled at each τ as follows:
Therefore, Equation 16 tested whether any kinematic variability remained in the firing residuals used to assess the modulation with respect to XE (Eq. 8). Similar PVS control regressions were performed for each error term.
The second type of control regression assessed for an error parameter whether the firing residuals contained any information related to the complementary error terms. Referred to as the error control regression, for XE it is modeled as follows:
Therefore,Equation 17 evaluated whether any firing variability related to YE, RE, and PDE remained in the firing residuals used to assess the modulation with respect to XE (Eq. 8). Similar error control regressions were performed for each error term.
Because the PVS control regressions fit neural data to the PVS model, the kinematic space (Xc, Yc, VXc, VYc, Sc) was partitioned into five equal bins along each dimension. Xc and Yc ranged from −6 to 6 cm, VXc and VYc ranged from −12 to 12 cm/s, and Sc ranged from 0 to 12 cm/s resulting in 3125 5D bins. Note this is the same binning used in our previous study that evaluated the simple spike firing with the PVS model (Hewitt et al., 2011). This partitioning was used for all regressions fitting firing residuals (Eqs. 16, 20, 21, and 23) or firing (Eq. 22) to the PVS model.
Decoding of errors based on the linear regression results.
Using results from the linear regression analysis, we determined whether the simple spike firing accurately predicts errors at the population level. For each cell, all movement trials were randomly divided into training (80%) versus test (20%) trials. For each cell the linear regression analyses were repeated for parameters XE, YE, and RE (Eqs. 12–14) using only the training trials. For each parameter, only cells with a significant R2 peak at a negative τ value (predictive peak) were included in the population analysis. Using only test trials to compute the predicted errors, we inverted Equations 12–14 at the τ value corresponding to the predictive peak. The predicted errors were then weighted by the ratio of the cell's predictive R2 peak to the average of the predictive R2 peaks. For each error parameter the decoded and observed values were pooled across the cell population and sorted by the observed values into 20 equal bins spanning the parameter space (twice the resolution used for the original regression analyses). Bins were averaged to obtain both the error population code and the observed error values.
Decoding of PDE followed a similar approach, but used a population vector analysis (Georgopoulos et al., 1982; Schwartz, 1993). Again, training trials from significantly tuned cells were used to assess the τ value and the corresponding PDE directional tuning (i.e., preferred PDE unit vector derived from Eq. 15) for each cell. Data from the test trials, were sorted into 15° bins based on the observed PDE and averaged. Matched firing data were shifted to correspond to the value compared with the predictive peak and also averaged. For each bin (j) and cell (i) we computed a weight wi,j = (Di,j − Di)/(Dmax − Di), where Di,j is the mean firing rate of the bin j, Di is the geometric mean for cell i, and Dmax is the maximum of cell i. A population vector could then be calculated for each bin by averaging the preferred PDE unit vectors determined from the training trials weighted by the firing modulation during the test trials (wi,j). The direction of the population vector is the decoded PDE value in each bin. For each error parameter, the decoding analysis was performed 25 times, each time randomly selecting the training and test trials. Final results are the averages of these 25 runs.
Analysis of simple spike firing modulation with a multilinear error model.
An analysis combined the error parameters into a single model, referred to as the multilinear error (ER) model. To determine the error-related modulation independent of kinematics, we used the firing residuals in which kinematic variability was removed, FRτPVS(t) (Eq. 6). FRτPVS(t) was fit to the ER model:
The actual firing was also fit to the ER model:
Note this is the same regression model used to generate the firing residuals, FRτER(t) (Eq. 7).
We also computed a control regression to evaluate whether kinematic firing variability was completely removed from the firing residuals. The control regression, FRτPVS (PVS(t)), is defined as follows:
The ER models used only a single τ (Eqs. 18–20). As addressed before (Hewitt et al., 2011), it is computationally prohibitive to directly calculate a unique τ value for each error parameter in a combined model. Therefore, we used the regressions on the individual parameters (Eqs. 8–11) as approximations for the contributions of each parameter with independent timing τ. To validate this assumption, the sum of the four individual error parameter R2 profiles was compared with the temporal profile of the ER model. The degree of similarity was quantified using the root mean square deviation (RMSD). This analysis also provided an additional test of the mutual independence of the error parameters.
Analysis of the simple spike firing modulation with kinematics.
The final regression analysis assessed whether the previously reported kinematic encoding (Hewitt et al., 2011) is independent of the error parameters evaluated in this study. To evaluate this, the firing residuals in which error variability was removed, FRτER(t) (Eq. 7), were fit to PVS model:
For comparison to the firing residuals results, the firing was also regressed to the PVS model:
Again we note that this is the same regression model used to generate the firing residuals, FRτPVS(t) (Eq. 6).
Finally, we computed a control regression to evaluate whether firing variability with the error parameters was completely removed from the firing residuals, FRτER(t). This regression, FRτER (ER(t)), is defined as follows:
Results
Behavioral analyses
The probability density functions for the four error parameters are illustrated in Figure 1B–D. The highly symmetric probability densities for XE and YE show that the cursor has a maximal probability of being at the target center and there is no preference for any quadrant. The probability density distribution for RE shows a gradual decrease with distance from the center (Fig. 1C). For both position and distance errors, the low levels of probability beyond ±1.5 cm reflect brief excursions outside the target. These excursions illustrate task difficulty and were permitted by the paradigm. Across all recording sessions, the probability of being outside the target was 0.12. PDE measures the divergence between the position error vector and cursor movement direction. The mode of 173° indicates that the monkeys preferentially aim their movements toward the target center (Fig. 1D). Therefore, the probability density functions for each error demonstrate that the monkey's strategy is to track the target center, supporting the assumption that the error parameters are strongly related to task performance.
Task error definitions and behavioral analyses. A, Geometric representation of the position, distance, and direction error parameters. The target (gray disk) moves on a pseudorandom trajectory (blue trace) and the cursor (red trace) tracks the target. B, Probability density as a function of XE and YE. White circle represents the target. C, Probability density distribution as a function of RE. The gray box denotes when the cursor was outside the target. D, Probability density distribution as a function of PDE. The probability density distributions in B–D are based on all trials of random tracking obtained during the Purkinje cell recordings from both monkeys. Only tracking period data are included.
Simple spike modulation with task errors
The neuronal data consists of simple spike activity recorded from 120 Purkinje cells during the random tracking task. This is the same population of cells used to document the kinematic signals encoded (Hewitt et al., 2011). The present study establishes that the firing of these neurons is also simultaneously modulated by the four error signals. Figure 2 illustrates the encoding of both error and kinematic parameters for two example Purkinje cells. The firing of each cell modulates strongly with position errors XE and YE (Fig. 2A,G). The simple spike firing of cell 23 is maximal when the cursor is on the upper right quadrant of the target and cell 41 firing is greatest when the cursor is positioned on the left side of the target. Tuning with RE is also present for both Purkinje cells (Fig. 2B,H). Cell 23 exhibits a linear increase in firing with RE and cell 41 a linear decrease. The simple spike firing of cell 41 has well defined tuning to PDE with the preferred direction at 194° (Fig. 2I). The discharge of cell 23 also exhibits significant directional tuning with its preferred direction at 151° (Fig. 2C).
Purkinje cell simple spike modulation in relation to task errors and hand kinematics. A–C, Simple spike modulation for cell 23 is plotted in relation to task errors: XE and YE (A), RE (B), and PDE (C). Plots of the simple spike firing with each error parameter based on the optimal τ obtained using firing residuals (Eqs. 8–11; Fig. 3). D–F, Simple spike firing for cell 23 is plotted in relation to position, velocity, and speed: XC and YC (D), VXC and VYC (E), and SC (F). Plots are at the optimal τ using the PVS model (Eq. 22). G–L, Similar plots of the simple spike firing with errors and hand kinematics for cell 41. For the plots of firing to RE, PDE, and speed, the averaged firing (± SEM) and the best fit regression line are plotted. The simple spike modulation is plotted relative to the mean firing (see Materials and Methods). White box in the position error plots (XE and YE) indicates bins with no data points. For both cells, the simple spike modulation is significant (p < 0.001) for all parameters except for the speed (F and L, p > 0.05). White circles in A and G represent the target.
Simple spike firing of the same two Purkinje cells is significantly modulated with hand kinematic parameters as previously described (Hewitt et al., 2011). The simple spike firing of cell 23 has strong tuning with position (Fig. 2D) and velocity (Fig. 2E). The firing of cell 41 modulates strongly with velocity (Fig. 2K) and to a lesser extent with position (Fig. 2J). Although the firing of cell 23 tends to increase with speed, neither cell has a statistically significant relationship between firing and speed (Fig. 2F,L). These examples highlight that Purkinje cell simple spike firing modulates with varying combinations of both kinematics and task performance errors.
Fit of firing to individual error parameters
As described (see Materials and Methods), we assessed the correlation between the firing residuals and each error parameter after removing the contributions of the kinematic parameters and the three error terms (Eqs. 8–11). In the two example cells, Figure 3A–H shows the R2 temporal profiles based on firing residuals (red) and actual firing (green). The two methods result in highly similar R2 profiles for each error parameter, both in timing and amplitude. The population data confirm the similarity in peak R2 values (Fig. 3I–L) and indicate that the individual error parameters are independently signaled in the firing discharge.
Fit of simple spike firing to the individual error parameters. A–D, R2 profiles for cell 23 as a function of τ obtained from fitting the actual firing (green; Eqs. 12–15) or the firing residuals (red; Eqs. 8–11) to each error term. For regressions based on firing residuals, the red arrows denote the times of the significant R2 peaks. Also plotted are control regression R2 profiles obtained from fitting the firing residuals for each term to either the PVS model (PVS ctrl, black; Eq. 14) or the complementary error parameters (ER ctrl, blue; Eq. 15). E–H, Similar plots of the R2 profiles for cell 41. I–L, For each parameter and regression, the average maximal R2 value (± SD) across all Purkinje cells is plotted.
The results show that each error parameter makes important contributions to the simple spike modulation and is widely represented in the population. Regression analyses based on the firing residuals reveal that 93% (112/120) of Purkinje cells has a significant fit with XE (average R2 = 0.19 ± 0.17; Fig. 3I) and 111 of 120 cells have a significant fit with YE (average R2 = 0.14 ± 0.14; Fig. 3J). For 108 Purkinje cells, the firing residuals have a significant fit to RE with an average R2 of 0.09 ± 0.09 (Fig. 3K). Finally, the firing residuals of 100 Purkinje cells have a significant fit to PDE, with an average R2 of 0.06 ± 0.05 (Fig. 3L).
Also plotted in Figure 3 are the control regression results that assess whether firing variability related to the kinematics or the complementary error terms was effectively removed from the firing residuals. The R2 temporal profiles for the control regressions are essentially flat for each error parameter (Fig. 3A–H, blue and black traces) demonstrating that there is no kinematic or complementary error information remaining in the firing residuals throughout the −500 to 500 ms time window. The average peak R2 values for the control regressions are below the threshold level (R2 ≤ 0.02) required for significance (Fig. 3I–L). This finding confirms that the R2 values based on the firing residuals cannot be explained by interactions between the kinematics and complementary error parameters. Therefore, the regression analyses based on the firing residuals demonstrate two important results. First, the modulation due to the kinematics and the complementary error parameters is effectively removed from the firing residuals (Eqs. 8–11). Second, the four error parameters are represented in the simple spike discharge and these representations are mutually independent and also independent of hand kinematics.
Bimodal timing
An important and novel feature of the R2 temporal profiles is the presence of two prominent peaks. Bimodal timing is evident in both example cells (Fig. 3A–H). For cell 23, XE, YE, and PDE have two significant R2 peaks and for cell 41, XE and PDE have two significant R2 peaks. For each example, the first peak occurs at a lead and the second at a lag. As the regressions were performed using firing residuals, the bimodal timing is not due to the influence of kinematic parameters or interactions with other error terms. Bimodal timing demonstrates there are two distinct times of peak correlation between the simple spike firing and an error parameter. The peak R2 at a lead suggests simple spikes carry a prediction about an error parameter, while the peak R2 at a lag suggests the simple spikes encode feedback about the same parameter.
To test for the prevalence of bimodal timing across the population of Purkinje cells, the number and timing of local maxima (i.e., times at which the first differential of the R2 profile changed from positive to negative) were determined from the regressions based on firing residuals (Eqs. 8–11). A bimodal profile was defined as having significant R2 peaks at both negative and positive τ values. Bimodal timing was common, with 86 Purkinje cells (72%) exhibiting bimodal timing for at least one error parameter. For the 480 temporal profiles obtained from regressions with individual error parameters (4 parameters × 120 neurons), bimodal timing is present in 172 profiles. A small minority of the bimodal profiles exhibit more than one peak at either lead (9%) or lag (11%) times. In these instances, analysis was restricted to the largest peaks. The timing of the peaks is uniform, averaging −223 ± 150 ms and 227 ± 132 ms for the lead and lag peaks, respectively (Fig. 4H). The average R2 for each error term is larger for the feedback peak than the lead peak (Fig. 4G), although the peak R2 amplitudes are nearly equal for PDE. This pattern suggests that, on average, the firing modulation due to feedback is greater than the modulation due to the prediction. However, for XE, YE, and RE the average peak R2 at lead times is at least 60% of the amplitude of the average R2 at lag times. Therefore, bimodal timing is a common property of how simple spike firing encodes errors. Furthermore, at both leads and lags, error parameters can account for a significant fraction of the simple spike variability.
Simple spike sensitivity for the error parameters. A–C, Value of the regression coefficients (βs) versus R2 for XE, YE, and RE at the τ value corresponding to the lead or lag peak. D, DM versus R2 for PDE. E–F, Distribution of preferred directions for PDE at lead peaks (E) and lag peaks (F). G, Average lead/lag peak R2 (± SD) for each error parameter across all cells. H, Average lead/lag τ values (± SD). All regression results in Figures 4⇓–6 are based on fitting the firing residuals to the individual error parameters (Eqs. 8–11). In each panel, gray represents lead peaks and black represents lag peaks.
Examination of the βs at the time of the peak R2 values shows that firing sensitivity with XE, YE, and RE is proportional to the variability explained by the regression model (Fig. 4A–C). This relationship held whether the peak R2 occurred at a lead (gray circles) or lag (black circles). Similarly, the DM for the cosine fit to PDE scales with the peak R2 (Fig. 4D). The distributions of the preferred directions are bimodal, with maxima centered at 0 and 180° for both lead and lag timing (Fig. 4E,F). These maxima are consistent with PDE being a valid estimate of direction error, with a preferred direction at 0° reflecting increased firing with maximal error and a preferred direction at 180° reflecting increased firing when error is minimal. Therefore, the simple spike firing is highly modulated by performance errors and this error-related modulation can be of large amplitude for both predictive and feedback timing.
Sign reversal for the bimodal regression coefficients
Because the simple spike firing can encode both predictive and feedback representations for an error parameter, this raises the question of how the discharge is modulated at lead and lag timing. To investigate this question, the βs were determined as a function of τ for XE, YE, and RE (Fig. 5). For cell 23, XE, YE, and RE have bimodal R2 profiles (Fig. 3A–C) and the βs for XE, YE, and RE evolve between negative and positive values (Fig. 5A–C). For these three parameters, a negative β occurs at the τ value corresponding to the lead maximal R2, while a positive β occurs at the lag maximal R2. The confidence intervals do not include zero, confirming that the βs at each peak are significant (Zar, 1999). Cell 41 has a bimodal profile for XE (Fig. 3D) in which the β reverses sign (Fig. 5I), with a positive β at the lead (τ = −140 ms) and a negative β at the lag (τ = 400 ms). A positive β implies that the simple spike modulation increases in relation to the error and a negative β that the modulation decreases.
Examples of dual representation of error parameters in the simple spike discharge. A–D, Regression coefficients (βs) of XE, YE, RE, and the preferred direction for PDE as functions of τ for cell 23. Error bars indicate the confidence intervals for each parameter at the times of the peak R2s (Fig. 3). E–H, Firing modulation at the peak τ values for each error parameter. In E the modulation is shown at the peak lead time and in F at the peak lag time for YE. In G and H, the firing modulations at lead and lag times are superimposed in the same plot. I–L, Temporal profiles for the regression coefficients and preferred direction for cell 41. M–P, Plots of the firing modulation at the lead and lag times of the peak R2s for cell 41. The simple spike modulation depicted in the position error plots (M and N) is based on the peak τ values for XE.
A sign reversal in the βs for XE, YE, and RE denotes that the predictive and feedback signals have opposing (i.e., anticorrelated) influence on the simple spike modulation. To demonstrate this property, the simple spike firing profiles were examined at the τ values corresponding to the lead and lag peaks for each error parameter. For cell 23, the change in sign of βXE and βYE is consistent with the higher simple spike firing on the lower left side of the target at a lead of −300 ms and maximal firing on the upper right quadrant of the target at a lag of 260 ms (Fig. 5E,F, respectively). The modulation with RE (Fig. 5G) has a negative slope for the predictive firing (τ = −280 ms) and a positive slope for the feedback (τ = 40 ms). For cell 41, the sign reversal of βXE results in greater simple spike modulation on the right side of the target for the lead and greater simple spike firing on the left side of the target for the lag. The R2 temporal profiles for YE and RE for cell 41 have a single peak and each occurs at a lag (Fig. 3F,G, respectively). Accordingly, YE and RE are characterized by a single β at the times of the peak lag (Fig. 5J,K). Therefore, for RE only the simple spike firing and linear fit at the time of the peak lag is plotted (Fig. 5O).
Because the fit to PDE involves sine and cosine terms (Eq. 11), the simple spike firing modulation is not linearly related to the βs. Instead, the preferred direction can be used to describe the simple spike modulation to PDE. Therefore, we defined a sign reversal for PDE if the preferred direction changed by >90°. Cell 23 has a preferred direction of 150.8° for the lead (τ = −100 ms) and 354.2° for the lag (τ = 200 ms), a change in the preferred direction of 156.6° (Fig. 5D). As expected, the PDE for cell 23 exhibits anticorrelated firing patterns for the prediction versus the feedback representations.
For some bimodal profiles, β has a constant sign or the preferred direction does not change >90°. This is illustrated for cell 41 for the fit to PDE. The preferred direction has a bimodal profile with peaks at −140 and 360 ms (Fig. 5L). However, the sign of βs are constant. Therefore, the firing modulation with PDE is similar at the times of the predictive and feedback peaks, as assessed by the preferred direction and DM (Fig. 5P).
The reversal of sign or change in preferred direction at the peak τ values is a common and significant finding, occurring in 74% (128/172, χ2(1) = 41.0, p < 0.0001) of the bimodal profiles. The modulation strengths for XE, YE, and RE for the predictive and feedback signals are strongly correlated for both types of profiles, with a correlation coefficient of 0.85 for the ones that change signs and 0.91 for those with constant sign (Fig. 6A, filled blue circles). Similarly, the change in preferred directions shows a clear tendency to segregate into two groups (Fig. 6B). For cells that exhibit reversal of the preferred direction (red) and those that do not (blue), the change in preferred direction clusters toward 180 and 0°, respectively, with fewer values ∼90°. These findings suggest the two types of profiles may represent distinct functional groups.
Dual representation of the error parameters at the population level. A, Distribution of the standardized regression coefficients (βstd) at the lead/lag peaks for XE, YE, and RE. Blue dots depict the βstd that reversed sign (anticorrelated firing modulations) and red dots the βstd in which the sign remained constant between lead and lag peaks (correlated firing modulations). B, Distribution of preferred direction difference between the lead and lag peaks for PDE. Blue dots depict differences >90° (anticorrelated firing modulations) and red dots depict differences ≤90° (correlated firing modulations). C, Distribution of the correlation coefficients (ρ(Firing)) between simple spike firing modulation at the times of lead/lag peaks for each bimodal profile. Blue bars are bimodal profiles with sign reversal of the regression coefficients or >90° changes in preferred direction. Red bars depict bimodal profiles with constant sign. D, Average correlation (± SD) between error parameters (ρ(Parameters)) at the peak τ values. Average positive and negative correlations are computed separately. E, Distribution of the ρ(Parameters) shown in D for bimodal profiles with sign reversal (blue) and sign constant (red). F, Plot of ρ(Firing) versus ρ(Parameters) separated into bimodal profiles that exhibit sign reversal (red) or constant sign (blue).
To quantify how bimodality influences the simple spike discharge, for each bimodal profile the correlation coefficient between the firing patterns, ρ(Firing), at the lead and lag τ values was determined. A ρ(Firing) of −1.0 implies complete reversal of the firing pattern. For cell 23 (Fig. 5A–H), the firing correlations are −0.73 for XE, −0.86 for YE, −0.63 for RE, and −0.78 for PDE. This demonstrates that the predictive and feedback components have strongly opposing effects on the simple pike modulation for the three parameters. For cell 41, XE is the only error parameter that reverses sign and ρ(Firing) is −0.38. For PDE the preferred direction does not change >90° between the lead and lag peaks and ρ(Firing) is 0.36. Across the population, there is a spectrum of correlation coefficients (Fig. 6C). Bimodal profiles with sign reversal predominantly generate negatively correlated simple spike firing patterns (95/128, 74%) while bimodal profiles with a constant sign typically result in positively correlated firing (37/44, 84%). This is further evidence that bimodal profiles with sign reversal and those without may reflect two populations of Purkinje cells.
A final issue is whether bimodality is generated by an inherent structure in the error parameters. To address this question, we performed four additional analyses on the error parameters with bimodal profiles. First, we determined whether there was any correlation within an error parameter at the times of the peak R2s. The correlation coefficients (ρ(Parameter)) were computed and averaged separately for negative and positive values. For the four error parameters, there is almost no evidence for negative correlations and the positive correlations are <0.2 (Fig. 6D). The low correlations suggest that there is little covariation in an error parameter, either negative or positive, at the times of maximal modulation with simple spike firing. Second, the distributions of correlation coefficients show that the weakest correlations occur most frequently, whether a profile exhibits sign reversal or constant sign (Fig. 6E). The third analysis examined the relation between the correlations in the firing patterns (ρ(Firing)) and correlations in the parameters (ρ(Parameter)) at the times of the peak R2s (Fig. 6F). There is no relationship between the correlation in firing and the correlation in the parameters, demonstrating that there is no structure within the errors that accounts for the bimodal profiles and the associated simple spike modulation. Finally, we computed the autocorrelation function for each bimodal profile. If an error parameter has correlations at the time of the peak R2, one would expect structure in the autocorrelation at those times, such as a local minimum or maximum. The autocorrelations for 171 of the 172 bimodal profiles exhibit no such structure. Instead, the autocorrelation functions have a single maximum at 0 time lag with smooth fall-off for negative and positive τ values. Therefore, there is no evidence that correlations within an error generate the bimodal profile or the manner in which the bimodal profile modulates the simple spike discharge.
Multiple linear error model
We also regressed the simple spike firing to the ER model (Eq. 18). The analysis was based on the firing residuals (FRPVS(t + τ); Eq. 6) after removal of the variability associated with the hand kinematics. The R2 profiles (black) show that the optimal τ value occurs at a lag for the two example neurons (Fig. 7A,B). Each cell also has a prominent peak at a lead. This bimodal profile mirrors the findings from the individual regression models (Fig. 4) in which the highest R2s occur at a lag but there is also a large R2 peak at lead times.
Integration of error and kinematic signals in Purkinje cell simple spike discharge. A–B, Examples from two Purkinje cells of R2adj temporal profiles obtained by fitting the firing residuals, FRPVS(t + τ), to the ER model (black traces; Eq. 18) and sum of the individual R2 temporal profiles obtained based on fitting firing residuals to individual error parameters (gray traces; Eqs. 8–11). Similarity between the two profiles was quantified by RMSD. Cell 23 is shown in A and cell 60 in B. C, Distribution of the RMSD values across the population of Purkinje cells. D, Average maximal R2adj (± SD) based on fitting the simple spike firing (F(ER(t)); Eq. 19) and firing residuals with kinematic removed (FRPVS(ER(t)); Eq. 18) to the ER model. The third bar is the average maximal R2adj (± SD) for the control regression that fit the firing residuals with kinematic variability removed to the PVS model (FRPVS(PVS(t)); Eq. 20). E, Average maximal R2adj (± SD) calculated by fitting the simple spike firing discharge (F(PVS(t)); Eq. 22) and the firing residuals with errors removed (FRER(PVS(t)); Eq. 21) to the PVS model. The third bar is the average maximal R2 for the control regression that fit the firing residuals with errors removed to the ER model (FRER(ER(t)); Eq. 23). F, Scatter plot of the optimal R2adj values from fitting the firing residuals with kinematic variability removed to the ER model (Eq. 18) versus the optimal R2adj values from fitting the firing residuals with error variability removed to the PVS model (Eq. 21).
The next analysis tested if the individual error models (Eqs. 8–11) provide a valid estimate of the ER model (Eq. 18). This was determined by summing the individual R2 profiles for the four error parameters and comparing the sum to the profile obtained for the ER model. As shown, the sum of individual R2 profiles (Fig. 7A,B, gray traces) provides a very good approximation to the R2 profile of the ER model (black traces). The small RMSDs confirm that the individual error profiles nearly sum to the multilinear profile. The additivity is conserved across the population of Purkinje cells with a mean RMSD of 0.024 ± 0.013 (Fig. 7C). Therefore, the results obtained using the ER model can be decomposed into the sum of the single parameter models. This finding of additivity provides further evidence that the error parameters are independent. The additivity also provides additional evidence the bimodal timing is not due to interactions between the two peaks and strong justification for using the individual error models to assess the contributions of each parameter to the simple spike firing. This is particularly important for determining the timing for an individual parameter. As noted previously (Hewitt et al., 2011; see Materials and Methods), it is computationally prohibitive to determine the leads and lags for the individual error terms directly from the ER model.
The population results reveal that the ER model based on the peak R2 explains on average ∼ 31% of the simple spike firing variability whether using the firing residuals or the actual firing (Fig. 7D). The control regression demonstrates that the variability related to the error parameters was effectively removed by using the firing residuals. We also reassessed the encoding of kinematics described earlier (Hewitt et al., 2011) using firing residuals in which the variability due to the error parameters was removed. The PVS model explains ∼38% of the firing variability, whether based on the firing residuals or firing (Fig. 7E). It is worth emphasizing that for the ER model, the kinematic-related variability has been removed and for the PVS model, the error-related variability has been removed. Therefore, both classes of signals modulate the simple spike discharge independently.
The results show that both error and kinematic signals are simultaneously present in the simple spike discharge of a single Purkinje cell (Figs. 2, 7). This raises the question of whether simple spike activity typically encodes both classes of parameters or if a single cell preferentially encodes one class over the other. Comparing the R2adj of the ER and PVS models based on the firing residuals reveals no evidence of preferential encoding (Fig. 7F). Across the range of R2s, the simple spike firing is similarly modulated by both classes of signals, confirming that error and kinematic information is integrated within a single Purkinje cell.
We also examined the relation between error and kinematic encoding within cells. Determining these relationships is complicated by the fact that the simple spike firing of a single Purkinje cell is significantly correlated with multiple error and kinematic parameters. For example, 45% of Purkinje cells (54/120) significantly encode at least three of the four error parameters at predictive timing and 67% of Purkinje cells (80/120) at feedback timing. Only 18% of Purkinje cells encode a single error parameter at prediction timing and only 7% of Purkinje cells at feedback timing. Similarly, the firing of 71% of the Purkinje cells encodes four or more of the five kinematic parameters. Only one Purkinje cell encoded a single kinematic parameter. Therefore, each cell signals a complex combination of errors and kinematics that could result in numerous interactions among parameters.
To reduce this complexity, we restricted our analysis to the dominant error and kinematic parameters for each Purkinje cell (i.e., parameters with the largest R2 value). Kinematic encoding was based on the PVS model (Eq. 21) and error encoding on the individual error regressions (Eqs. 8–11). Using a contingency analysis, we determined for each kinematic parameter the probability of a significant difference between the observed and expected frequency of error parameters (χ2 test, χ2 < 0.05). The expected frequency was defined as the frequency at which an error parameter is dominant in the 120 cell population. The observed frequency was defined as the frequency that the same error parameter is dominant given that a kinematic parameter is also dominant. Four significant relationships were found among the kinematic and error parameters: (1) VXc shows an increased frequency with EX, (2) VYc with EY, (3) Sc with ER, and (4) Yc with PDE. VXc and VYc coupled with XE and YE suggests a natural pairing in which Purkinje cells that most strongly modulate with velocity respond strongly to errors along the same axis. Similarly, speed could fluctuate with the distance of the hand from the center of the target (RE) and whether the hand is moving toward or away from the target (PDE).
Decoding of error parameters
A population analysis assessed whether the error signals predicted by the simple spike firing matches the actual errors. As shown in Figure 8, the average estimates of the predicted errors obtained from the population of Purkinje cells are proportional to the observed errors for all four parameters. Determination of the correlation coefficient and slopes from the individual data points shows the slopes for XE, YE, and PDE are close to unity. Although the slope for RE is <1, there is still a strong linear relationship between the decoded and observed RE values. The SDs of the decoded estimates are quite small, showing that the estimates are relatively invariant to the selection of the training and test trials, particularly when the hand is within the target boundaries. However, the variability increases as the error magnitude increases for XE, YE, and RE. One factor in the increased variability is there are typically less data points at the larger error values. Fewer data points will decrease the accuracy of both the encoding and decoding. Another explanation for decreased decoding performance toward the boundaries of the error space is that Purkinje cells continuously monitor these parameters within a limited range. The estimates of the predicted errors are valid even when using a higher spatial resolution than used for the linear modeling showing the continuous nature of the error signals. Therefore, the simple spike firing in a population of Purkinje cells carries a highly accurate prediction of the actual errors. Finally, individual cells make predictions about the upcoming errors at different τ values ranging from −500 to −20 ms. This suggests that Purkinje cells provide prediction of the errors across a range of lead times, possibly to allow integration of motor commands having different timing requirements with the error prediction at the appropriate τ values.
Population decoding of the predicted error parameters. Plots of decoded predicted versus observed values for each error parameter. A, XE (ρ = 0.95, p < 0.0001, slope = 0.99). B, YE (ρ = 0.83, p < 0.0001, slope = 0.98). C, RE (ρ = 0.76, p < 0.0001, slope = 0.52). D, PDE (ρ = 0.98, p < 0.0001, slope = 1.05). Each plot is the average of the predicted errors (mean ± SD) versus the observed errors from 25 repetitions of the off-line decoding algorithms.
Discussion
Simple spike discharge encodes error signals
The present study demonstrates that Purkinje cell simple spike discharge in the intermediate zone encodes a wealth of information about task errors. Although previous studies are consistent with error or instructive signals in the simple spike firing (see Introduction), and functional imaging studies document cerebellum activation with movement errors (Blakemore et al., 2001; Diedrichsen et al., 2005; Miall and Jenkinson, 2005), this is the first demonstration that the simple spike discharge encodes specific performance errors. Importantly, the results establish that the error signals are independent from the kinematic signals encoded in the simple spike firing (Fortier et al., 1989; Gomi et al., 1998; Roitman et al., 2005; Medina and Lisberger, 2009; Hewitt et al., 2011).
This finding is at odds with the prevailing hypothesis that climbing fibers provide error or teaching signals in the cerebellum (see Introduction). However, the relation between complex spike discharge and errors remains problematic. Although error signals might coexist in both the complex and simple spike firing, it is not clear how the low-frequency complex spike discharge could encode, at least at the single cell level, the type of continuous error signals evaluated in this study (Ebner and Bloedel, 1987; Simpson et al., 1996). In several studies complex spike error-related signals were observed in a small fraction of trials and only after extensive averaging (Gilbert and Thach, 1977; Ojakangas and Ebner, 1994; Kitazawa et al., 1998). Several studies failed to find the expected responses to errors or perturbations (Horn et al., 1996; Ebner et al., 2002; Catz et al., 2005; Dash et al., 2010). A recent study suggests a novel role for complex spikes during learning with climbing fiber input signaling sensitivity to errors, instead of encoding errors (Marko et al., 2012). This hypothesis is based on the observations that motor learning decreases with error size as does the probability of complex spike occurrence with visual errors during saccades (Soetedjo et al., 2008; Marko et al., 2012). Finally, the complex spike discharge occurs in response to errors, that is, only as a feedback signal. We are not aware of any evidence that complex spikes predict sensory consequences, as shown in the present study for the simple spike firing.
Notably, the error measures evaluated are highly task relevant and provide continuous evaluation of tracking performance. These measures are similar to error signals involved in human motor adaptation (Berniker and Kording, 2008; Izawa and Shadmehr, 2011). Furthermore, use of these measures distinguishes the present study from most previous cerebellar electrophysiological investigations into error signaling. All four error parameters are significantly represented in simple spike firing from the majority of Purkinje cells. Position and distance errors explain more firing variability than the direction term. This difference probably reflects the quality of the sensory feedback. The target edge provides unambiguous information for the position and distance errors, while PDE is inherently a more complex measure and lacks direct sensory feedback. Nevertheless, many cells exhibit significant firing modulation with PDE and the preferred directions cluster ∼0/360° and 180°, directions at which the error magnitude is maximal or minimal. The distribution of preferred directions also suggests that directional errors are highly relevant. The decoding analyses show that, at the population level, Purkinje cells carry the information needed to predict the upcoming errors.
Purkinje cell simple spike activity integrates kinematic and error signals
Another finding established by analyzing firing residuals is that error and kinematic signals are integrated at the single cell level. The overlap of kinematic and error information occurs in functionally relevant pairs. There is no evidence for segregation by signal type among Purkinje cells. These observations may explain why motor adaptation in the presence of continuous sensory input differs from adaptation based only on endpoint or reward feedback (Magescas et al., 2009; Izawa and Shadmehr, 2011). The predictive component of the error signals could function as internal feedback used to control movements, thus avoiding the problems of delayed feedback and accounting for cerebellar involvement in compensating for motor command variability (Bhushan and Shadmehr, 1999; Wolpert and Ghahramani, 2000; Xu-Wilson et al., 2009; Shadmehr et al., 2010). The feedback component of the error signals could function to modify a forward model during motor learning (Wolpert and Ghahramani, 2000; Shadmehr et al., 2010). The presence of independent kinematic and error signals suggests the cerebellar contribution to motor control involves computations on two interacting models. The first model provides accurate kinematic information, including predictions of upcoming movements. This model is stable across several tasks (Roitman et al., 2005; Hewitt et al., 2011). The second is an error model that deals with execution variability in real time and drives motor adaptation.
Dual representations of individual error parameters
Several motor control theories require the computation of sensory prediction errors, in which the predicted and actual sensory feedback for a given motor command are compared (Wolpert and Miall, 1996; Wolpert and Ghahramani, 2000; Shadmehr et al., 2010). Both eye and limb movement studies suggest sensory prediction errors are the critical signals that drive motor adaptation, not the actual motor corrections (Wallman and Fuchs, 1998; Noto and Robinson, 2001; Mazzoni and Krakauer, 2006; Tseng et al., 2007). Importantly, it is the difference between the expected and actual sensory consequences that influences adaptation (Wong and Shelhamer, 2011). It is tempting to speculate that the dual representation of the error parameters provides the signals needed to compute sensory prediction errors. In 70% of Purkinje cells at least one R2 profile exhibits two local maxima, one with predictive and one with feedback timing. Several lines of evidence show that the bimodality cannot be explained by an inherent structure in the error parameters. For the majority of these dual representations, the regression coefficients reverse the sign for the predictive and feedback timing, resulting in opposing modulation in the simple spike firing (Figs. 5, 6). If appropriately summed, the anticorrelated predictive and feedback simple modulations could be used to compute a sensory prediction error. As required, the resulting sensory prediction error would be in the same kinematic space as a forward model (Ariff et al., 2002; Flanagan et al., 2003). Although for many bimodal profiles the correlations between spike firing at the predictive and feedback timing cluster close to −1, there is also a spectrum of weaker correlations (Fig. 6C). This suggests there are varying degrees of mismatch between the prediction and feedback about a specific error. Mismatch is expected due to the unpredictability of random tracking inducing frequent errors.
Two fundamental questions are where and how this comparison occurs. Computation of a sensory prediction error requires that the two components come together in space and time (Wolpert and Ghahramani, 2000). In this study, the predictive and feedback representations in individual Purkinje cells are separated by several hundred milliseconds. This suggests that computation of a sensory prediction error based on comparing the predictive and feedback signals would have to occur downstream of Purkinje cells, possibly in the cerebellar nuclei.
The implications of regression coefficients that have the same sign at both local maxima are more ambiguous. One possibility is that bimodal profiles with the same sign represent a particularly poor prediction that grossly mismatches the feedback. Another possibility is that these dual representations are a different class of computation resulting in the enhancement of a signal as opposed to cancellation. The latter is suggested by the sharp separation in slopes and preferred direction differences between the error parameters with constant sign versus sign reversal (Fig. 6A,B). Additionally, the distribution of correlation coefficients for firing patterns cluster toward 1 and −1, respectively (Fig. 6C). If the correlated predictive and feedback simple modulations were summed, this would act to highlight a signal. It has been shown that cerebellum-like structures have the capability to minimize sensory reafference and highlight other signals (Bell et al., 2008). Therefore, the cerebellum may be involved in the generation of sensory prediction errors and in the enhancement of salient sensory information.
Questions remain regarding how the CNS uses the error signals described and we acknowledge that correlation analysis cannot determine causality. However, having identified robust error signals allows for testing of specific hypotheses. For example, we hypothesize that the prediction and feedback components in a highly predictable task would be more closely matched and result in better cancellation of the two signals at the Purkinje cell. Another testable hypothesis is that motor learning would shift the weightings of the predictive and feedback signals. At the start of learning, the predictive and feedback signals would be poorly matched. As learning proceeds, the predictive signal would be expected to provide a better estimate of the feedback. This cancellation mechanism would be similar to the negating of sensory reafference described in cerebellar-like structures (Bell et al., 2008).
Finally, the finding of dual representations in the simple spike discharge for all error parameters suggests this type of encoding is a general principle of how Purkinje cells function. The results are consistent with cerebellar involvement in implementing a forward internal model, integrating predictions about the consequences of a motor action with the sensory feedback to fine-tune motor behavior and guide motor adaptation (Wolpert and Miall, 1996; Wolpert et al., 1998; Ebner and Pasalar, 2008; Shadmehr et al., 2010).
Footnotes
This work was supported in part by National Institutes of Health Grants R01 NS18338, R01 NS18338-25S1, F31 NS071686, and T32 GM008244. We thank Michael McPhee for generating graphics and Kris Bettin for preparation of this manuscript. Dr. Siavash Pasalar participated in the beginning phase of this study. We also acknowledge the technical assistance of Bethany Kummer, Gabriel Rodriguez, and Hao Wang. Dr. Claudia Hendrix provided statistical advice.
- Correspondence should be addressed to Dr. Timothy J. Ebner, Department of Neuroscience, University of Minnesota, Lions Research Building, Room 421, 2001 Sixth Street S.E., Minneapolis, MN 55455. ebner001{at}umn.edu
