The model

From a computational point of view, the calculation of a motor signal, such that the executed movement equals the desired one, requires the biophysical features of the limb to be integrated into the pre-motor circuits. Not only the biological properties of the muscles (e.g. stiffness, viscosity) must be considered, but also the inertial properties of each limb, the reciprocal interaction forces between segments and their interactions with external forces such as gravity.

This could be achieved by means of internal inverse models of the biophysics of the moving limbs, embedded in neural networks [16]–[18] (Figure 1A). Such an inverse computation can be performed by a neural circuit (Figure 1B) composed of two parallel, closed, internal feedback loops: a positive one (with a gain close to, but smaller, than one to insure stability) and a negative one containing an internal forward model of the direct biophysical function of the limb (denoted H* in Figure 1B), that processes motor commands and computes predictive signals anticipating sensorimotor signals [15], [16], [19].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Main principles of the cerebellar architecture of the model.

For the sake of clarity, only the command circuit for one muscle is illustrated. (A) Structure of a command circuit accounting for the physical constraints. θ^D: desired movement; θ^P: performed movement; u: neural command. H: direct function incorporating all biophysical features of the limb. H⁻¹: internal inverse model of the direct function. (B) Control scheme used to compute an approximate inverse function. The two summing elements (positive/negative inputs) represent the cerebellar nuclei (CN) and the red nucleus (RN). H*: internal forward model of the direct function H. P and Q represent the signals originating from the cerebellar cortex (CC) and CN, respectively. (C) Direct functions Π* (representing the mechanical constraints, e.g. gravity, inertia) and μ (representing the muscle features) in the external world, and their counterparts in the CNS labeled Π* and μ*. These two internal forward models (Π* and μ*) are embedded through two internal feedback loops placed in series to calculate their approximate inverse, i.e., Π⁻¹and μ⁻¹. The direct pathways convey signals of desired position θ^D and forces F^D. In the indirect pathways, the negative output of the elements computing Π* and μ* are comparable to the inhibitory projections of the Purkinje cells of the cerebellar cortex to the neurones of the cerebellar nuclei. The P signals are comparable to simple spike activities of Purkinje cells. Dashed lines represent the climbing fibers. SC: spinal cord. Lower scheme: artificial neural networks simulating the CC. Granular, Golgi and Purkinje cells (respectively 8, 1 and 1 for each predictor) are modelled by formal neurones. s and s⁻¹ represent respectively the derivative and integration in the Laplace domain. α represents multiplicative higher orders of the position. The adaptive elements and connections are represented in grey, fixed elements in black. For the sake of the clarity only the first neural network is represented.

https://doi.org/10.1371/journal.pone.0005176.g001

This circuit, by means of its structure formed of two short and parallel feedback loops, one of which containing a predictor, computes an approximate internal inverse model of the direct biophysical function of the limb. This is nevertheless achieved without performing an explicit inverse calculation. Together, these two loops provide a feed-forward control. From an anatomical viewpoint, the internal forward model H* is thought to correspond to the cerebellar cortex, which receives sensory signals (Figure 1B, C) as well as efferent copies of motor commands (u) through the mossy fibers. The resulting Purkinje cell activity represents the simple-spike inhibitory signal that is sent to the cerebellar nuclei, and is thought to predict torque based on the motor command. This inhibitory signal, together with the efferent copy of the neural command (u), is then fed to a summing element which would represent the cerebellar nuclei (labeled CN in Figure 1B,C). The efferent copy loops through a second (downstream) summing element (labelled RN in Figure 1B–C) which represents the magno-cellular part of the red nucleus. This last summing element adds the output signal issued from the summing element representing the cerebellar nuclei to the signal coding the desired movement. It is noticeable that the motor command (u) reaches both the predictor and this summing element, similar to the messages conveyed by the excitatory mossy fibers that reach both the cerebellar cortex and the cerebellar nuclei.

Both the inverse kinematics and the inverse dynamics problem need to be solved in order to compute an adequate neural command. However, since the robotic arm used in our experiment has only one geometrical degree of freedom (DoF), the inverse kinematic problem is not of interest here. Nevertheless, the robotic arm is actuated by a pair of antagonist muscles which requires the computation of the inverse dynamics at two stages: first, the inversion of the biophysical characteristics of the muscles (including viscosity and stiffness) and second the inversion of the biomechanics of the moving segment, (including the physical constraints applied to the arm such as the gravitational and inertial torques). Notably, these inverse computations must deal with a combination of non-linearities due to the actuators and due to the different movement directions (up, down) in the gravitational field.

The motor command (denoted u_i) computed in this model (Figure 1A–C) can be compared to the activity of a pool of motoneurons allowing muscular contractions of the i^th muscle, which produces a force f_i, by means of a biophysical process described by the direct function denoted μ_i. It must be noted that here i ∈{1,2} since we consider two muscles (for the sake of clarity, in Figure 1 we used simplified notation since only one circuit command allowing to compute the motor command for one muscle is illustrated). Then, these forces act on the joint and produce the resulting torque to accelerate the movable segment. Here, the torque T_i is related to the force f_i with respect to the rotation centre O, according to the following equation:(1)where “×” denotes the cross product (for simplicity arrows above vectors have been omitted) of two vectors r_i, whose origin is the insertion point of the muscle and whose extremity is the point of application of the force f_i on the moving limb. According to Newton's law, the equation expressing the total torque exerted at the joint by the actuators and the mechanical forces applied to the one DoF arm moving in the sagittal plane is as follow:(2)(3)(4)where T_f1 and T_f2 are the respective torques developed by the muscle 1 and 2. J is the moment of inertia and , and , are respectively the angular position, velocity and acceleration of the limb. B and K denote respectively the coefficients of viscosity and stiffness of the joint. Finally, m, g and ρ denote respectively the mass of the limb, the gravitational acceleration and the radius drawn from the axis of rotation to the gravity vector. The direct transfer function of the dynamic mechanical constraints (e.g. gravity, inertia) applied to the limb are denoted by Π (Figure 1C). It must be noted that the direct functions μ and Π act in series, since the first one provides the force exerted by the muscles (given the neural command), whereas the second provides the movement caused by the torque resulting from the various forces applied to the arm (see the right side of Figure 1C, labeled “World”). Elaboration of the motor commands requires the successive computations, in the reverse order, first of the approximate inverse function of Π (labelled Π⁻¹), which provides the desired force from the desired movement, and second of the approximate inverse function of μ, (labelled μ⁻¹), which provides the neural command from the desired force (see the left part, labeled “CNS”, of Figure 1C).

The inverse computation of these two functions μ⁻¹ and Π⁻¹ is performed by means of the general computational scheme depicted in Figure 1B. Therefore, the neural control circuit shown in Figure 1C includes two internal forward models (denoted μ* and Π* mimicking respectively the two direct functions μ and Π in order to predict the desired movement and force) put in series and embedded within two distinct cerebellar neural network modules. The predictors μ* and Π* were modelled using artificial neural networks whose architecture was designed by replicating the well-known connectivity of the cerebellar cortex [15], [16] (for reviews see [20], [21]).

Motor learning.

The synapses between parallel fibers and Purkinje cells are well-known neural sites of plasticity in the cerebellum, modulated by teaching signals arising from the Inferior Olive and conveyed via climbing fibers. A supervised learning procedure schematizing this circuitry was incorporated into our model. Based on previously results, no other sites of cerebellar plasticity were introduced [15], [18]. Learning was modeled by modifying progressively the input weights of the parallel fiber signals, according to the correlation of these signals with a teaching signal calculated from the quadratic error between the desired and performed positions and forces (Figure 1C).(5)(6)

A teaching signal was sent whenever the difference in position was superior to 0.5° (simulated) or to 1° (for robot movements, Figure 2B) and the difference in force >20 g according to the following learning rules:(7)(8)where and .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Learning curves for robotic experiment.

(A) Three learning curves obtained for the three mass conditions M0, M1 and M2. (B) Typical desired and performed arm displacements from a horizontal position (i.e. 0°) and with an amplitude of 25°. DaD: Desired angular Displacement. PaD: Performed angular Displacement.

https://doi.org/10.1371/journal.pone.0005176.g002

Details and demonstrations of the learning rules are given in [16], [18]. Typical learning curves with an asymptotic profile are shown in Figure 2A.

Experimental design

Motor task and gravitational constraints applied to the robot arm.

Specific combinations of movement direction, initial position and amplitude of limb motion (i.e. joint angle) in the sagittal plane can lead to important non-linear variations on the gravitational effects exerted on the limb [5]. Therefore, the dynamics of the limb cannot be computed with straightforward classical command methods or by a simple inverse calculus. As a consequence, in order to test the learning capacities of this cerebellar model to compute the inverse dynamics of gravity torques, we manipulated the gravitational effects applied to the arm through three experimental sessions during which various initial positions, movement amplitudes and directions were combined during vertical arm pointing movements. For both the simulated and the real robotic limb, vertical arm movements were performed to specific targets (Figure 3): up to 40°, 35°, 30°, 25°, 20°, 15°, 10°, 5°, 0°, and down to −5°, −10°, −15°, −20°, −25°, −30°, −35°, −40°. The origin (0°) of this frame of reference was located at the axis of rotation of the robotic arm, which corresponds to the human shoulder joint level. Figure 3 illustrates the three behavioral sessions; each session is distinguished by a particular variation of the gravitational torques.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Schematic explanation of the setup and the three simulated tasks.

Representation of the simulated arm and its two antagonist muscles. F1: force developed by muscle 1. F2: force developed by muscle 2. R: radius of the sprocket. Black arrow: gravitational torque (GT) exerted on the segment. First row: Session I. Initial position (green circle) at 0°; Required movements: upward (red arrow), downward (blue arrow). Traces to the right: gravitational torque over time as a function of target movement amplitude and direction (upward: red, downward: blue). For all movements of session I, the gravitational torque varies monotonically and is independent of movement direction. Second row: Session II. Initial positions at 20° and −20°. The gravitational torque varies monotonically but depends on movement direction. Third row: Session III. Initial positions at 20° and −20°. For each movement (amplitudes: 40° and 60°) the gravitational torque varies non-monotonically.

https://doi.org/10.1371/journal.pone.0005176.g003

In the first session, all movements started from the same initial position, namely at 0°. Eight upward (i.e. 5°, 10°, … 40°) and eight downward (i.e. −5°, −10°, … −40°) pointing movements towards targets at multiples of 5° were performed. Under these conditions, the gravitational torque is a decreasing monotonous function of movement amplitude, but is independent of movement direction (Figure 3A, compare red and blue gravitational torques traces). Thus, the gravitational torque is identical for up and down movements of similar amplitude, although the arm moves against gravity (upward movements) and with gravity (downward movements).

In the second session (Figure 3B), arm movements were performed from the initial positions of 20° and −20°. From these initial positions, the arm had to perform 20° rotations, i.e. from the initial position of 20° down to 0° or up to 40° and from the initial position of −20° up to 0° or down to −40°. Under these conditions, the gravitational torque is a function of movement amplitude and direction, as well as of the initial position, and varies for a given direction in a monotonous fashion. For example, from the initial position of 20°, a downward displacement implies an increase in gravitational torque, whereas the same movement started from the −20° position involves a decrease in gravitational torque. In addition, the amplitude of these variations is not only opposite but also differs in magnitude. For instance, from the 20° initial position, when the arm moves upward and downward, the gravitational torques decrease and increase, respectively (i.e. opposite changes); in addition, the amounts of these changes were different (due to the non linearity of the gravitational torque). Since the initial arm positions of 20° and −20° are symmetrical with respect to the horizontal axis, gravitational torque values are identical at movement onset for movement in session II.

In the third session (Figure 3C), arm pointing movements were performed from the same initial positions as in session II, however, movement amplitudes were at least twice as large (i.e. 40°, 60°). Under these conditions, the gravitational torque varies again as a function of movement amplitude, direction and initial position, but now in a non-monotonous fashion since the arm crosses the horizontal plane. The gravitational torque increases from the initial position until the horizontal axis is reached and then decreases towards the target. Furthermore, for all sessions, the inertia of the arm was systematically varied by applying additional weights to its extremity. This imposed further changes in the magnitudes of the gravitational torques. We tested to which extent the simulated and real robot arm were able to perform accurate pointing movements under these experimental conditions. Performance of the simulated or real robot arm was assessed after learning.

Simulation and robotic experiments

Note that the learning phase was carried out only during the session I and for the 3 (i.e. 10°, 20° and 30°) of the 8 amplitudes. In addition, learning took place under different mass conditions. In the vertical plane the increase of the inertia of the moveable segment implies changes in both inertial and gravity torques and constitutes a good example for testing the capacity of our model to integrate both intrinsic (inertia) and extrinsic (gravity) parameters. Therefore, the capacity of generalization of our cerebellar neural network was tested in two ways: 1/In the first session, the performance in the other amplitudes than those used for training were assessed, i.e. interpolations to 15° and 25° and extrapolations to 5°, 35° and 40° for up and down directions. 2/In the second and third session, by quantifying (the difference between desired and performed) movements which differed significantly from those imposed during the learning procedure in the first session. Note that while adaptations to changes in gravitational force can be tested empirically, the initial learning of internal representations of gravity can (currently) only be studied by means of simulations and robotic experiments using anthropomorphic limbs.

Simulated task

Experimental procedure with the real robotic arm

The robot experiments were strictly identical to those performed during the “Simulation 1”. However, since the mass at the tip of the segment could not be changed after each trial for technical reasons, robot experiments similar to those presented during “Simulation 2” could not be performed.

Data analysis

Motor performance was assessed by taking into account both the dynamic and static components of the movement. The static component was defined as the period when the displacement reached and remained within ±5% of its final value. The dynamic component was defined as the transient displacement from the movement onset up to the beginning of the static component. For the dynamic and static components, the performance was quantitatively assessed by the RMSE (Root Mean Square Error) between the desired and performed displacement according to the following formula:where MD is the movement duration, θ^D and θ^P represent the desired and performed angular displacement, respectively. The RMSE corresponding to the dynamic component of the displacement provided by the model was labeled RMSE_D, the static component was labeled RMSE_S. Whereas the RMSE_D can be considered as transient and related to the dynamics of the system (e.g. delays, stiffness), the RMSE_S is observed once the system reached its steady state and is similar to the constant errors analyzed in studies investigating arm pointing movements in humans (for a review see [23]).

In order to compare the effect of the experimental conditions (mass, amplitude, direction and initial position) on RMSE_D and RMSE_S, a statistical analysis has been conducted. According to the normality of the distribution of the data (tested using a Lillefort's test) we used either parametric (e.g. t-test, ANOVA) or non-parametric (e.g. Wilcoxon, Kruskall-Wallis tests) statistical methods.

Simulation 1

The goal of the first simulation was to establish the performance of the proposed model after learning one specific weight matrix for each of the three mass conditions. Figure 4 (first column) shows performed (dashed trace) and desired (continuous trace) movement amplitudes over time for all movements and sessions under the mass condition M0. Figure 4A depicts the upward (red) and downward (blue) movements for the session I (first row) for the session II (second row) and for the session III (third row). Clearly, the network successfully learned to control movements against and with gravity.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Performance of the model (M0 condition) for simulation 1 and the robotic experiment.

Left column: arm movements performed during simulation 1 in the sagittal plane with no additional mass (M0) during the first (first row), second (second row) and third session (third row). Performed (dashed trace) and desired (continued trace) angular displacement. Red and blue traces: up and down movements respectively. Second column: Arm movements under zero gravity conditions after learning under normal gravity condition, Session I, II and III, no additional mass (M0). Third column: Robot arm movements in the sagittal plane with no additional mass (M0) during the first, second and third session.

https://doi.org/10.1371/journal.pone.0005176.g004

Motor performance after learning (under training conditions).

Table 1 (column ‘Training’) shows that after learning the overall performance of our model was good; the average RMSE_D and RMSE_S between the desired and performed movement ranged between 0.04° and 0.08°. Statistical analysis did not reveal any significant effect (t-tests p>0.05). Notably, there was no significant difference between RMSE_D and RMSE_S errors. Interestingly, values of RMSE_D and RMSE_S were similar for the three weight conditions (matrixes). Furthermore, movement direction had no significant impact on the RMSE_D (up: 0.07±0.02°; down: 0.06±0.02°) and RMSE_S (up: 0.05±0.02°; down: 0.07±0.02°) values.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Pointing errors for simulation 1.

https://doi.org/10.1371/journal.pone.0005176.t001

Simulation 2

The goal of Simulation 2 was to examine motor performance after having learned a single weight matrix, which had, in addition, to cover a wider range of mass conditions (6 different masses up to a 50% increase). The generalization over these different mass conditions was tested.

Motor performance after learning (under training conditions).

The grand average RMSE_S (for various masses and movement amplitudes) was 0.48±0.32°, that is, larger than for simulation 1 (t-test, p<0.001). Similarly, the grand average of the RMSE_D increased to 0.44±0.30° compared to simulation 1 (t-test, p<0.001) but was not significantly different from the RMSE_S (Wilcoxon test, p>0.60). However, and in contrast to simulation 1, the performance depended of the mass. Table 2 (column ‘Training’) shows the average RMSE_S and RMSE_D for each mass condition. Clearly, the best performance (i.e. lowest RMSEs and RMSE_D) was found for intermediate masses (M2_T, M3_T), while the performance for small (M0_T, M1_T) or great masses (M4_T, M5_T) significantly increased. These differences were statistically significant (Kruskall-Wallis test, p<0.001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Pointing errors for simulation 2.

https://doi.org/10.1371/journal.pone.0005176.t002

Motor performance for novel movement (intra- and extrapolation).

Table 2 shows the average RMSE_D and RMSE_S for each mass condition and experimental session. Similar to the training session, the average performance was better for intermediate weights and poorer for light and heavy weights. Figure 5 shows the RMSE_D (left column) and RMSE_S (right column) for each of the three sessions. Each graph depicts the RMSE as a function of mass condition (i.e. for intra- and extrapolated masses) and as a function of movement amplitude and initial position.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Performance of the model for simulation 2.

Pointing error as a function of weight (mass condition M0_Iep, M1_Iep, M2_Iep, M3_Iep, M4_Iep and M5_Iep) and movement amplitude (and initial position). First column: dynamic error (RMSE_D). Second column: static error (RMSE_S). Each black dot represents an error measure (dynamic or static) obtained for a given experimental condition. IP: initial position.

https://doi.org/10.1371/journal.pone.0005176.g005

Session I: with values of 0.44±0.30° and 0.45±0.32°, the grand averages (i.e. including all experimental conditions) of RMSE_D and RMSE_S were similar to the training condition (Table 2). As in the training condition, both RMSE_D and RMSE_S varied as a function of mass. This is illustrated in Figure 5 top row, which shows a V-profile of the RMSE, indicating that the RMSE varied as a function of mass, but not as a function of movement amplitude or movement direction. This was statistically confirmed using a Kruskall-Wallis test coupled with a multiple comparisons test for each factor (for both type of error; mass: p<0.001, amplitude and direction: p>0. 95).

Session II: the grand averages of RMSE_D (0.06±0.01°) and RMSE_S (0.15±0.04°) were smaller than those in session I (t-test, p<0.001) and smaller than those in the training session (Table 2). Figure 5 (middle row, note reduced Z-axis scale) shows that the RMSE_D varied little (p>0.10 for all comparisons) whereas the RMSE_S clearly varied as a function of mass and of movement direction (p<0.0014) but not of initial position and direction (p>0.70).

Session III: the grand averages of RMSE_D (0.15±0.05°) and RMSE_S (0.39±0.24°, Table 2) tended to be smaller than those in session I (t-test, p<0.01, only for the RMSE_D) but larger than those in session II (t-test, p<0.01, for both type of RMSE). Figure 5 (bottom row) shows changes of RMSE_D and the RMSE_S as a function of mass, of movement amplitude and initial position. In both cases, the largest RMSE was generally found for movements of large amplitudes and either heavy or light masses. However, no significant effect of the mass, movement amplitude and initial position on both type of error was found (p>0.23 for all comparisons, except for the case where the RMSE_D was higher for a U movement of 60° amplitude than for a D movement of 40° amplitude (p = 0.009).

Robot experiments

The aim of the robot experiments was to verify whether the simplified calculations of the gravity torque acting on the segment would be adequate for controlling an anthropomorphic robot arm. Only the equivalent of simulation 1 was performed with the robot. Figure 4 (right column) shows the movement of the robot arm for all movements and sessions under mass condition M0.

Motor performance after learning (under training conditions).

Table 3 (column ‘Training’) shows that the overall robot performance was not as accurate as for the simulated movements. Over all three mass conditions the average RMSE_D (1.51±0.71°) was higher (t-test, p<0.05) than the RMSE_S (1.08±1.01°). Table 3 shows that the RMSE_D and RMSE_S did not vary as a function of mass, of amplitude or movement direction (Kruskall-Wallis test coupled with a multiple comparison, p>0.14 for all comparisons).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Pointing errors for robot experiments (conditions equivalent to simulation 1).

https://doi.org/10.1371/journal.pone.0005176.t003

Motor performance for novel movements (intra- and extrapolation).

Table 3 (column ‘Average Iep’) shows that performance for novel movements were less precise than those included in the sample used for training. However, the decrease was not statistically significant (t-test, p>0.10). As in the training condition, there was no significant difference between the three mass conditions (i.e. between the three matrixes). Furthermore, contrary to the simulation results, the RMSE_D and RMSE_S did not vary as function of the session. RMSE_D (averaged over the three weight conditions) for session I, II and III were respectively 2.01°±0.91, 1.79°±0.89 and 3.99°±1.34° while for RMSE_S, they were respectively 1.25°±0.93°, 2.49°±1.82 and 3.12°±1.94°. Lastly, motor performance of the robot arm did not vary as a function of movement direction.

Figure 6 compares the performance of simulated and robot movements. Figure 6A shows RMSE_D and RMSE_S as a function of mass and movement amplitude conditions in session I (i.e. Inter and extrapolated position). For session II and III, RMSE_D and RMSE_S are illustrated as a function of initial position, movement amplitude and mass condition. Figure 6B shows the data for the corresponding robot movements. As already shown in Figure 4, the robot experiment produced errors larger than those in the corresponding simulation 1. Moreover, Figure 6 shows that errors on the three movement variables (amplitude, initial positions and directions of the movements) tended to be qualitatively different between the robot experiment and the simulation. For the simulation in session 1, maximal dynamic and static errors were found for the largest movement amplitudes (e.g. 40° up and down), independently of the mass condition. Although the robot experiment replicated the independence on the mass condition, the dependence on amplitude changed: largest and smallest movement amplitudes provoked maximal errors (e.g. 40° as well as 5° and 10°). For session II, except for the absolute values, no striking qualitative difference was found between simulation and robot experiment. In session III, the simulation showed maximal errors preferentially for large movement amplitudes and specific mass conditions. In contrast, the robot experiment seemed to be less sensitive to the mass condition and produced less symmetrical error surfaces.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Comparison between simulated and robot movements.

Distribution of the RMSE_D (left column) and RMSE_S (right column) for the three sessions for the simulation 1 (A) and the robotic experiment (B). Both type of error are represented as a function of movement amplitudes during the session I (i.e. intra- and extrapolated positions) and of initial positions and movement amplitudes during respectively the session II and III.

https://doi.org/10.1371/journal.pone.0005176.g006

The cerebellum computes both sensorimotor predictions and dynamic inverse computations

In our model, sensory and motor signals are used to predict the gravitational torque exerted on the arm and then to compute inverse dynamics of the arm. This is in agreement with previous studies which have suggested that the cerebellum implements an inverse model to perform the dynamic and/or kinematic inverse computation [9], [10], [24], [25]. Moreover, the proposed model also reinforces the idea proposed by prior investigations that the cerebellum computes sensorimotor predictions [11], [12], [26]–[28]. For instance, it has been suggested [29] that the cerebellum incorporates a forward model that predicts the sensory consequences of movements. The authors concluded that the cerebellum was involved in sensorimotor prediction and was therefore a plausible neural site to implement a forward model. Our model also supports the idea that the cerebellum includes both forward and inverse models [30]–[32], and it provides an anatomical plausibility and sensorimotor learning scheme for both internal models. Indeed, by taking into account the functional and the anatomical features of the cerebellar pathways, our model provides a coherent theoretical framework that implements the coexistence of these two internal models replicating both the functional role of the cerebellar cortex and its corresponding anatomical structure. Specifically, it shows that the cerebellar cortex could possibly implement a predictive neural network, which is equivalent to a forward model in order to predict the effects of the gravitational and inertial torques applied to the segment. In our model, this prediction is propagated through the cerebellar connectivity via three biologically plausible internal feedback loops: a first one in the cerebellar cortex and the two others in the cerebellar pathways. This model suggests that the cerebellar cortex is primarily responsible of sensorimotor predictions, whereas the entire cerebellum is in charge of the inverse dynamics computation. Hence, it is compatible with previous studies that have proposed that the cerebellum integrates an inverse model, as well as with studies that have suggested a predominant role of the cerebellum in sensorimotor prediction.

Computing an inverse model from a forward model is a technique used for industrial plant control [33] and has already been applied to motor control [12], [34], however, the proposed circuits were hardly biological plausible. On the contrary, the significance of the present model is that its structure is deduced from physical constraints, and it is consistent with the cerebellar connectivity. It controls arm pointing movements in a vertical plane based on the acquisition of a direct predictive model of the limb biomechanics, which is sufficient to approximate the inverse dynamics via multiple looped structures. It must also be noted that the principle of using such multiple recurrent loops, particularly through the forward model, makes our cerebellar model somewhat comparable to recurrent neural networks, reinforcing thus, the idea that this type of neural network structure is able to correctly simulate dynamic behavior [35]–[38]. Indeed, although these studies developed neural models without specifically addressing any particular anatomical structures; they revealed that recurrent neural networks efficiently learned temporal patterns by feeding back the copies of the current sensorimotor prediction outputs to the next sensorimotor inputs. Finally, one important theoretical and computational consequence of the architecture of this model is that it offers a biologically plausible solution that avoids the well-known, but artificial two-step sensorimotor learning scheme: most published models require, first, the learning of the forward model [34], or the assumption that the forward model is a priori known, and subsequently the learning of the inverse model [39], [40].

The importance of sensorimotor prediction of gravitational torques in neural control

Our model approximates inverse dynamics by means of predictions, i.e. by a forward model based on sensory and pre-motor signals. This internal forward model is acquired through motor learning and provides a control for vertical arm movements.

Previous studies, using simulation [34] or experimentation [7], [8], have highlighted the importance of internal forward models for sensorimotor predictions. For instance, using an object manipulation task, it has been shown that subjects learned sensorimotor predictions prior to the control of objects; the learning of a forward model precedes that of an inverse model [41]. Similarly, by means of simulations, it has been shown that even an approximate internal forward model was able to train an inverse model [34]. Furthermore, other investigations [8], [42], using a motor imagery paradigm, asked subjects to execute or to imagine horizontal and vertical single-joint arm pointing movements, and showed that movement durations were very similar for both conditions. Isochrony between executed and imagined movements is achieved by well trained internal forward models which precisely predict the gravito-inertial forces acting on the arm. More recently it has been shown that subjects improved 3D arm pointing movements during a speed/accuracy trade-off during mental practice [43]. The authors suggested that an efferent copy of the neural drive during motor imagination would be available to the forward model which, in turn, was thought to train the inverse model. These studies are in accordance with the present results indicating a major role of the forward model in motor control and learning.

Moreover, models of physical laws that are able to predict object motion in the gravitational field may also be represented in the brain [3] (for a review see [2]). For instance, during ball catching in the gravitational field, it has been found that the vestibular network and the medial cerebellum were activated [44] (for the implication of the vestibular system see [45]–[47]). It must be noted that the predictive mechanisms embedded in our cerebellar model is not in contradiction with the idea that the posterior parietal cortex could also be involved in motor predictions. Since both parietal cortex and cerebellum are anatomically connected through a functional loop [48], [49], they could share a neural substrate for forward models and may have different roles in motor prediction [50], [51]. Thus, prediction related to a high cognitive level could be implemented in the parietal cortex, whereas sensorimotor predictions (e.g. self-generated movements) may be performed by the cerebellum [12], [50].

Beyond the cerebellum and the parietal posterior cortex that seem to be the main candidates for forward models [31], [50], others studies, without always mentioning explicitly such a computational concept, revealed that others neural regions would also incorporate predictive mechanisms. For instance, the basal ganglia that projects to the cerebral cortex through multiple parallel channels [52], [53] would also perform predictions based on the environmental states to select the appropriate action in a given context [54], [55]. Also, highly neural recurrent structures such as the hippocampus [56], [57], the prefrontal [58], [59] and the orbitofrontal [60], [61] cortex seem to perform associative prediction that would reflect outcome expectancies, providing thus, an internal representation of the consequences likely to follow a particular act.

However, two main differences appeared between these predictive mechanisms and those embedded in the cerebellum and by extension in our model. First, our cerebellar model used a learning rule that simulated supervised learning that takes place in the cerebellum and driven by the climbing fibers messages issued from the inferior olive. Here, pointing errors were detected and used to compute the teaching signals, as it seems to occur in the inferior olive. This type of learning contrasts with the reinforcement learning used in the basal ganglia [54], [55] or hippocampus [57] and the unsupervised learning that takes place in the cortex [54]. However, more generally and independently of the learning procedure used, a key feature of our model is its capability to compute an inverse mapping by means of a direct mapping embedded into recurrent loops. Therefore, we can wonder if it might be possible that similar mechanisms, albeit a different learning method, might take place in others predictive structures as those mentioned to compute inverse mappings. Second, our cerebellar model predicts the sensorimotor output of the system contrary to the predictive mechanisms above described that perform predictions at a higher and more general behavioral level. Such lower and higher levels of prediction have been investigated in a series of computational studies in order to understand their interactions using multiples environments and tool manipulations [30], [32], [38], [62], [63]. The lower level would correspond to the sensorimotor processes of detailed environment interactions to provide an accurate control of limb motion. The higher cognitive level would embed abstractions of those lower sensorimotor processes level to infer behaviors or plan goal-oriented movements.

Our cerebellar model implemented predictive mechanisms related to the sensorimotor (short-timescale) predictions of the internal state of the arm, addressing thus, mainly the lower level. However, it is of interest to note that our model predicts the internal state of the arm (inertia) but also its mechanical interaction with the external (gravity) context in which actions takes place. In humans, the estimation of action context in terms of mechanical interaction between the body limb and the environment or a tool is essential for the performance of skilful movements. For instance, when we make a reaching movement while rotating our torso, we compensate for the velocity-dependent Coriolis forces that arise from the rotation and act on our arm [64]. Likewise, the development of a new forward model in the context of microgravity, allows astronauts to adapt their actions during space flights. Therefore, from a general point of view, the cerebellar model presented here takes in to accounts these dynamical interactions between a given environment and the body limbs but is not able to deal with multiple environments or tools. Thus, this model could be naturally embedded into a more general neural structure including a higher level. For instance, this could be explicitly done by means of a modular structure incorporating switching mechanisms such as gate-selection [30], [32], [38], or by means of an emerging functional hierarchy using neural networks with neurons having different timescale [63].

Performance and limits of our model of the cerebellar pathways

The movement accuracy achieved by our model during the robot experiment (RMSE_S between 0.33° and 5.81°) was comparable to previously reported data for human arm pointing movements with one [5] or two [65] geometrical DoF, with constant pointing errors less than 5°. The model was able to learn a set of upward and downward movements and to generalize by interpolation and extrapolation to other types of movements including those with non-monotonous profiles of gravitational torque. However, the pointing errors varied as a function of the type of generalization. First, the error increased when movement amplitudes were tested outside those of the training range. Second, the error also varied as a function of the time-varying profile of gravitational torque exerted during the movement: monotonic variations of gravitational torque (session I and II) provoked smaller errors than non-monotonic variations (session III). Furthermore, although learning was faster and produced smaller errors when based on a specific weight matrix per inertial condition, a global weight matrix across different inertial conditions provided much better generalization.

As previously mentioned, the two types of errors presented here refer to two different aspects of the performance of the model. Therefore, a reduced dynamic error reflects the capacity of our model to take into account the dynamic features of the system (e.g. non-linearities, delays, stiffness) whereas a small static or constant error would suggest the presence of a small and constant bias in the sensorimotor transformations [23]. Concerning the simulations 1 and 2, our model of cerebellar pathways was able to capture adequately the dynamics of the moving limb. However, the further the model had to extrapolate from the training conditions, the stronger was the constant bias in the sensorimotor transformations (even if both types of error remain small, <1.2°). Furthermore, whereas for the two simulations the dynamic error was always inferior or equivalent to the static error, during the robotic experiment session this trend was inverted (dynamic higher than static error). This is due to the fact that contrary to the simulation, when the learning is performed during the experimental session, some mechanical features (e.g. stiffness) were either not accurately taken into account or neglected by the model of cerebellar pathways.

Conclusions

This study presents a command circuit comparable to the cerebellar pathways that learns the inverse dynamics of an anthropomorphic robot arm, including the effects of the gravitational forces. Learning was achieved through an internal forward model allowing the computation of an approximation of the inverse dynamics. After learning, this circuit was able to drive arm movements in the vertical plane, with an accuracy comparable to that of human movements. The model suggests that the cerebellar cortex is a plausible neural site for learning internal predictive forward models of the gravitational forces, and that the whole cerebellum is likely able to perform approximate inverse computations.