J Neurophysiol 97: 331-347, 2007.
First published September 27, 2006; doi:10.1152/jn.00290.2006
0022-3077/07 $8.00
Computational Motor Control: Redundancy and Invariance
Emmanuel Guigon1,
Pierre Baraduc2 and
Michel Desmurget3
1Institut National de la Santé et de la Recherche Médicale (INSERM) U742, Action Neuroimagerie Modelisation, Université Pierre et Marie Curie, Paris; 2INSERM U534, Space and Action; and 3INSERM U371, Brain and Vision Research, Bron, France
Submitted 16 March 2006;
accepted in final form 24 September 2006
 |
ABSTRACT
|
|---|
The nervous system controls the behavior of complex kinematically redundant biomechanical systems. How it computes appropriate commands to generate movements is unknown. Here we propose a model based on the assumption that the nervous system: 1) processes static (e.g., gravitational) and dynamic (e.g., inertial) forces separately; 2) calculates appropriate dynamic controls to master the dynamic forces and progress toward the goal according to principles of optimal feedback control; 3) uses the size of the dynamic commands (effort) as an optimality criterion; and 4) can specify movement duration from a given level of effort. The model was used to control kinematic chains with 2, 4, and 7 degrees of freedom [planar shoulder/elbow, three-dimensional (3D) shoulder/elbow, 3D shoulder/elbow/wrist] actuated by pairs of antagonist muscles. The muscles were modeled as second-order nonlinear filters and received the dynamics commands as inputs. Simulations showed that the model can quantitatively reproduce characteristic features of pointing and grasping movements in 3D space, i.e., trajectory, velocity profile, and final posture. Furthermore, it accounted for amplitude/duration scaling and kinematic invariance for distance and load. These results suggest that motor control could be explained in terms of a limited set of computational principles.
|
INTRODUCTION
|
|---|
When we move our limbs to execute a motor task, we generally have many more degrees of freedom (DOF) than necessary to fulfill the requirements of the task. In a typical situation, unconstrained point-to-point arm movements involve 7 DOF for moving in a three-dimensional (3D) space. If movements are to occur in a plane, such as in a handwriting task, the physical constraint of contact still leaves infinitely many solutions to the problem of writing a letter. The coordination of kinematically redundant systems was first formulated by Bernstein (1967)
as the DOF problem. The main difficulty of Bernstein's problem is that the nervous system must conciliate two apparently conflicting abilities. On the one hand, each individual realization of a motor goal results from the choice of one among an infinite number of motor patterns. On the other hand, there is no univocal relationship between motor goals and motor patterns, a property first noted to be central to the functioning of motor systems by Lashley (1933), which he called motor equivalence. For instance, the nervous system can preserve a common kinematic pattern for executing a movement while varying the moving effector, the angular patterns of motion, or underlying muscular activations (Bernstein 1967; Burdet et al. 2001; Gribble et al. 2003; Levin et al. 2003). Solving Bernstein's problem should not only help us understand the origin of this flexibility, but also explain how the flexibility coexists with constraints inherent to the functioning of the motor system. For instance, when subjects are free to move, they automatically scale movement duration with movement amplitude and choose a trade-off between movement speed and accuracy (Fitts' law; Fitts 1954). An open question is how motor controllers in the brain solve Bernstein's problem. Here, we first present a computational approach to this problem, i.e., a solution cast in terms of principles. Our goal is to provide a solution based on a restricted number of realistic and well-supported hypotheses. Then we show how these principles apply to the problem of kinematic redundancy.
|
COMPUTATIONAL APPROACH
|
|---|
Breakthroughs into the understanding of motor functions have generally been brought about by computational studies (Bullock and Grossberg 1988; Feldman and Levin 1995; Flash and Hogan 1985; Harris and Wolpert 1998; Hoff and Arbib 1993; Todorov and Jordan 2002; Uno et al. 1989), i.e., studies that disclose functioning principles independent of brain structures or neural processes. However, since the time of Bernstein and Lashley, no adequate principled approach to kinematic redundancy and motor equivalence has been proposed and the way the nervous system tackles these problems remains mysterious (Gielen et al. 1995). In particular, the line of reasoning based on the separation between trajectory planning and trajectory execution, which attributes motor equivalence to the specification of a desired trajectory in task coordinates and which proposes a solution to redundancy using inverse kinematic transformations, has been seriously questioned (Bullock and Grossberg 1988; Cole and Abbs 1986; Sporns and Edelman 1993; Todorov 2004; Todorov and Jordan 2002). The main argument is that on-line movement corrections act to favor goal achievement rather than the following of a preplanned trajectory.
Principles
The goal of this article is to describe a principled approach that gives a unified account of motor behavior. We propose that motor control is governed by four principles: 1) separation: the nervous system processes dynamic (inertial, velocity-dependent) and static (elastic, gravitational) forces separately; 2) optimal feedback control: there exists an optimal controller for dynamic forces (dynamic controller) that is appropriate for an on-line control of movement, i.e., it is optimal for any initial state of the moving limb; 3) maximum efficiency: the nervous system attempts to reach the goal of the movement defined in the task coordinates with zero error and minimal control signals (these signals are inputs of central descending origin to motoneurons; overall control magnitude is called effort; see METHODS); 4) constant effort: all movements obeying a given set of instructions (e.g., move at preferred speed) are executed with the same level of effort. Maximum efficiency and constant effort are not incompatible. The former principle associates each movement duration to an effort level, i.e., the (minimal) effort of the optimal solution. When duration is not specified, the latter principle prescribes a level of effort used to determine the appropriate duration, i.e., the unique duration associated to this level of effort.
Below, we give a brief rationale for these principles. In the RESULTS section, we show how they are associated to properties of motor behavior.
SEPARATION PRINCIPLE.
There is substantial experimental evidence to support the idea of a separate processing of static and dynamic forces. Previous psychophysical studies showed that velocity profiles remain unchanged when moving in a known constant (Mustard and Lee 1987; Welter and Bobbert 2002) or a known elastic force field (Bock 1990; Flash and Gurevich 1992; Ghez 1979; Gottlieb 1996; Levin et al. 2003; Milner 2002; Rand et al. 2004), but that they are in general modified by time and amplitude scaling in velocity-dependent and inertial fields (Atkeson and Hollerbach 1985; Bock 1990; Gottlieb 1996; Happee 1993; Hatzitaki and McKinley 2001; Ruitenbeek 1984). This property is called kinematic invariance. These results show that, if a separation principle holds, it does not apply to any type of forces, but is specific to static forces. Even if some velocity-dependent perturbations have no influence on movement kinematics (Shadmehr and Moussavi 2000; Shadmehr and Mussa-Ivaldi 1994), our reasoning and conclusion remain valid. Electromyographic (EMG) studies revealed additive velocity-independent, tonic and velocity-dependent, phasic components that are related to the generation of antigravity torques and dynamic torques, respectively (Buneo et al. 1994; Flanders and Herrmann 1992; see also Farley and Koshland 2000; Milner 2002; Welter and Bobbert 2002). A similar additive combination between a tonic activity, related to the compensation of an external static force, and a phasic, movement-related activity was observed in neurons of primate motor cortex (Kalaska et al. 1989). The experiment of Nishikawa et al. (1999) showed that the terminal posture of 3D redundant movements is independent of movement velocity. Because the relative contribution of antigravity and dynamic torques varies with velocity, optimization of the total torque pattern would predict variations of terminal posture with velocity. This result suggests that dynamic forces are optimized independent of static forces. The study of Kurtzer et al. (2005a) provides further support to the separation principle. These authors showed that adaptation to a multiforce environment composed of a velocity-dependent force and a constant force was well described by a mechanism that processes velocity-dependent force separately from the total applied force. We emphasize that separation of static and dynamic forces is not separation of posture and movement because static forces (gravity, muscular elastic forces) are present during both posture and movement.
Optimal feedback control principle.
A solution to kinematic redundancy could likely be found in the framework of optimal control theory, which states that a unique solution to an ill-posed problem can generally be obtained as the solution that corresponds to the minimum of a cost function. Because it is well recognized that motor commands are continuously updated by internal feedback loops (review in Desmurget and Grafton 2000), it has been suggested that on-line control of movement is also optimal in the theoretical sense (optimal feedback control; Bryson and Ho 1975; Hoff and Arbib 1993; Todorov and Jordan 2002). This means that the control signals are optimal for any boundary conditions, such as following perturbation of the moving limb or changes in target position. The notion that motor control can be viewed as a continuous feedback process was first analyzed by Bullock and Grossberg (1988) and was shown to account for trajectory formation and kinematic invariance.
It is paradoxical that optimal control theory has been repeatedly and successfully used to account for many aspects of motor control (e.g., trajectory formation, muscular redundancy, postural control, locomotion; Anderson and Pandy 2001; Flash and Hogan 1985; Harris and Wolpert 1998; Kuo 1995; Uno et al. 1989), but has rarely been applied to the case of redundant manipulators (Todorov and Jordan 2002). Todorov and Jordan (2002) successfully applied optimal control to a kinematically redundant system. However, this work was limited to the linear case, i.e., a telescopic arm model rather than an articulated arm with a nonlinear dynamics. Also, no direct comparison was provided between experimental data and the predictions of the model. Does this mean that no appropriate solution can be found in this framework? In fact, the central difficulty for an optimal control approach to redundancy is to formulate the problem in such a way that it accounts for the simultaneous control of posture and movement. Most studies have not considered the case of static forces. Some studies that have actually tackled this problem have reported difficulties in solving optimal control problems in the presence of gravitational forces (Soechting and Flanders 1998; Thoroughman and Feller 2003). When a movement consists of a transition between two equilibrium postures, the boundary conditions of the optimal control problem should specify terminal equilibrium signals, e.g., muscle forces that compensate for applied static (elastic, gravitational) forces. The idea to add to the cost function a term that enforces given initial and final equilibrium postures (Dornay et al. 1996; Harris and Wolpert 1998) should lead to solutions that depend on the level and nature of the static forces. Although this issue has not been fully addressed in experimental studies, the results of Nishikawa et al. (1999) clearly suggest that dynamic and static controls are unlikely to be dealt with simultaneously. In contrast, the separation principle allows the application of optimal feedback control to kinematic redundancy problems with static forces because there is no a priori specification of the final posture of the limb.
An integral component of feedback control is the presence of a state estimator, i.e., a process that combines sensory inflow and motor outflow to provide an estimate of the current state of the system (Wolpert et al. 1995). A state estimator is necessary because 1) the state is in general not directly observable and must be inferred from sensory inputs (e.g., vision, proprioception) and 2) both inflow and outflow are noisy. It has been suggested that the nervous system acts as an optimal estimator (Baddeley et al. 2003; Wolpert et al. 1995). The optimal feedback control principle embeds an optimal estimator.
MAXIMUM EFFICIENCY PRINCIPLE.
Optimal control is associated with the choice of a cost function (which indicates a quantity to minimize) and a constraint function (which specifies constrains to satisfy) (Nelson 1983). Here, we consider a cost function based on the size of centrally generated signals that eventually generate the dynamic forces. The reason for this choice is twofold. First, it is simple and easily measurable by the CNS, compared with other cost functions encountered in motor control literature that require complex calculations and measurement processes (acceleration derivative, torque change, energy, etc.). Second, related cost functions were successfully used in recent models (Harris and Wolpert 1998; Todorov and Jordan 2002). As a constraint function, we use the initial and final boundary conditions. The rationale for this choice, which differs from the mixed error/effort cost of Todorov and Jordan (2002), is the following. The error/effort cost function contains parameters that weight the respective contribution of state errors (position, velocity, force, etc.) and effort in the cost function. Because different sets of parameters lead to different behaviors, a model based on error/effort minimization cannot provide a univocal description of motor control. This problem has little consequence for the study of motor variability (Todorov and Jordan 2002), but is critical for the study of redundancy.
CONSTANT EFFORT PRINCIPLE.
The idea of motor behavior being associated with the minimum of a cost function is appropriate when both movement amplitude and duration are specified. Otherwise, infinitely slow/fast or infinitely short/long movements could result. In cases where movement time is not given (e.g., instruction is to move at preferred speed), actual duration can be deduced by associating a desired level of effort with the instruction and using the relationship between amplitude, duration, and effort prescribed by the maximum efficiency principle. In this framework, the constant effort principle states that a given set of instructions is equivalent to a level of effort. For these instructions, movements of different amplitudes, directions, or against different loads are executed with the same effort.
A simple example
To illustrate the model, we consider the control of an inverted pendulum in the gravity field (Fig. 1A). The dynamics of the pendulum is given by
 | (1) |
where 
defines the position of the pendulum; I0, m0, and L0 are the inertia, mass, and length of the pendulum, respectively; and u is a control (here a torque). The control problem is to find u(t) (t in [t0; tf]) such that (t0) = 0 and (tf) = f (Fig. 1A). According to the separation principle, we can write
where ustat(t) compensates for the effect of static forces (see Eq. 3 below) and udyn(t) is the solution to the control problem in the absence of static forces. To obtain udyn(t), we apply the optimal feedback control principle and the maximum efficiency principle. At each time tc, a control udyn(t) is computed for the remaining duration (t in [tc; tf]) based on the planned displacement from the currently estimated position of the pendulum [^(tc), which is generally different from (tc) arising from noise or perturbation] to the target position f (Fig. 1B). The control udyn(tc) is applied and the process starts again at the next time step. The mathematical solution to this problem is obtained using the Euler–Lagrange equation (e.g., Hogan 1984; Kirk 1970)
 | (2) |
with boundary conditions (tc) = ^(tc) and (tf) = f. For this problem, we have
 | (3) |
which reflects the fact that the state of the pendulum is not known, but can only be estimated.
View larger version (16K):
[in this window]
[in a new window]
|
FIG. 1. A: inverted pendulum in the gravity field. B: in the presence of noise, the estimated trajectory of the pendulum (gray) deviates from its actual trajectory (black). At a given time tc, the control u(tc) is calculated based on the estimated error (gray arrow) and is applied to the current state of the pendulum (black arrow). C: application of the model to a displacement between 0 and 90° in 1 s. Left: angle of the pendulum, velocity in inset. Right: dynamic (plain) and static (dashed) control. For simplicity, we used I0 = 1, m0 = 1, L0 = 1. D: same as C in the presence of additive Gaussian noise ( = 0.4). Left: amplitude/duration scaling. Dotted lines: 90° in 1 s. A nonzero intercept (here 0.05 s) was necessary for appropriate functioning of the controller. Right: angle of pendulum. Dotted lines: 1 s, 90°.
|
|
In the absence of noise and perturbations, we see that the control calculated at t0 for the whole movement duration [t0; tf] (open-loop control) is the same as the control recalculated at each time step (closed-loop control). An example of solution to Eq. 2 in the absence of noise and perturbations is shown in Fig. 1C. An open question is whether the controller is stable in the presence of noise. We simulated the same example in the case where
where n(t) is a zero-mean Gaussian noise with variance 2. Because of noise, the pendulum will in general never reach its target position. Thus movement duration could not be fixed in advance and was determined by an amplitude/duration law—i.e., at each time, the remaining duration was a function of the residual distance to the target (Fig. 1D, left). The pendulum reached the vicinity of the target and then smoothly oscillated around it (Fig. 1D, right). This result was robustly observed over a broad range of parameters (variance of noise, slope of the amplitude/duration law). Thus stable control of a naturally unstable system can be obtained with the principles of the model in the linear case. Stability arises from the action of the dynamic controller, which attempts to reduce the distance between the estimated state of the pendulum and the target, but not of the static controller, which simply compensates for gravity.
General sketch of the model
The model is summarized in Fig. 2A. It has three components: 1) a dynamic controller, which calculates, for given target (goal) and estimated states, the appropriate control to master the dynamic forces and progress toward the goal (dashed arrow on the left); 2) a static controller, which calculates for each estimated state the appropriate control that maintains equilibrium against the static forces (dashed arrow on the right); 3) a state estimator, which calculates a state estimate from sensory inflow and motor outflow (gray arrows).
Scope of the article
The four principles have been proposed to address Bernstein's problem. However, the present article concerns the part of the problem related only to kinematic redundancy. The goal of this article is to show that the dynamic controller can solve the problem of kinematic redundancy in a way that is consistent with experimental observations. Thus we focused on the dynamic controller and we used the following simplifications (Fig. 2B). First, because neither perturbations nor noise was introduced, no state estimator was needed and optimal feedback control was strictly equivalent to open-loop optimal control (see Fig. 1B). Thus simulations described below were obtained using a method to solve open-loop control problems. Second, the static controller was not modeled. We assumed that there were no static forces or, equivalently, that the static forces were exactly cancelled at each time step during the movement. The complete model including static and dynamic forces was not simulated.
|
METHODS
|
|---|
Modeling approach
The principles were cast in computational terms so they needed to be translated into simulations for comparisons with experimental observations. A central issue was the representation of the neuromuscular system. In fact, it is a complex machinery that contains both intricate neural circuits and noisy nonlinear elements for action and sensation, so an open question is the degree of biological realism that can guarantee that a simulation is a useful reflection of biological operations. Although there is no general answer to this question, we addressed it using two principles: 1) start with a simple model; if it does not work, then try a more complex one; and 2) of two models that provide similar results, the simplest is the best.
According to the first principle, we built on simple hypotheses. At a basic level of details, we focused on four main properties of the muscle: 1) it generates forces when it is stimulated; 2) it behaves like a linear spring; 3) it performs low-pass filtering; and 4) it is inserted around a joint and generates a torque. We also assumed that 1) a muscle is made of a single type of fibers that are all innervated by the same motoneuron; 2) a motoneuron receives a unique and specific control signal (see Controlled object). Because linear springlike forces are static forces, they were removed according to the separation principle. The results obtained with these hypotheses were in adequate correspondence with experimental data. For comparisons, we built more complex models. A muscle is in fact a nonlinear spring so we included a force/velocity relationship. A muscle can act on several joints so we considered biarticular muscles. A muscle can be controlled through muscular synergies. We first describe the results obtained with the simplest model that best highlight the power of the proposed principles. Then we show results obtained with biarticular muscles (see Muscular redundancy), nonlinear muscles (see Influence of muscle properties), and muscular synergies (see Muscular synergies).
Controlled object
The controlled object was modeled as a rigid, multilink, articulated system with N DOF. It was actuated by one pair of antagonist muscles per DOF. Each muscle i (1 
i 2N) was controlled by a motoneuron and the ensemble motoneuron + muscle (neuromuscular system) was modeled as a second-order low-pass filter (van der Helm and Rozendaal 2000), which transforms a neural control signal (ui) into a muscular force (Fi) according to
 | (4) |
where 
is a time constant and 
(z) = [z]+ ([z]+ = z if z > 0; otherwise, [z]+ = 0). This latter function was used to express the fact that a muscle exerts only a pulling force. The terms for these quantities ei and ai correspond to excitation and activation, respectively. Electromyographic (EMG) activity was defined as [ei]+. Torques were calculated at each DOF as the difference between the forces generated by antagonist muscles scaled by a coefficient (
; see following text).
The dynamics of the controlled object and the neuromuscular system was described by
 | (5) |
where x describes the state of the object (angular position and velocity) and the muscles (activation, excitation). Equation 5 contained Eq. 4 for each muscle and the equations for movement dynamics (see Kinematic chains).
Optimal controller
Point-to-point trajectories of the controlled object were obtained as solutions of an optimal control problem: find deterministic controls u(t) = {ui(t)} (1 i 2N) in [t0; tf] such that x(t) is a solution to Eq. 5 with boundary conditions
 | (6) |
and the quantity
 | (7) |
is minimum (we use the generic term effort for the quantity E; Todorov and Jordan 2002). Function 
expresses constraints on the final state of the system (see following text). Any kind of constraint can be handled. For instance, was x(tf) – xf for nonredundant objects (xf is the desired final boundary value).
Constant effort principle
The optimal control model can be used to calculate control signals that drive the arm from an initial position to a target position (movement of amplitude A) in a given time (duration D = tf – t0). By construction, the control signals are unique and their associated effort is E (Eq. 7). We can repeat this operation for different amplitudes and different durations and build a surface E = E(A, D). This relationship is monotonic because it increases with A and decreases with D (larger or faster movements require larger controls). Conversely, if A and E = Ec are given, it defines a unique duration D. The relationship between movement amplitude and duration corresponding to a given effort (Ec) was obtained by building the surface E = E(A, D) and searching the intersection between this surface and the plane E = Ec. The surface was obtained by interpolation between calculated points on a grid.
Kinematic chains
We considered several kinematic chains that illustrate different aspects of redundancy (Fig. 3). A kinematically nonredundant system was used to study kinematic invariance and muscular redundancy (Fig. 3C). A chain is specified by a set of joint angles (qi, i = 1...N) that define a kinematic transformation (from angular to Cartesian coordinates). The equations for movement dynamics were derived using the Euler-Lagrange method (Spong and Vidyasagar 1989). Briefly, four steps were necessary: 1) the rotation vector of each segment was calculated; 2) translational and rotational kinetic energies were calculated; 3) the Lagrangian L was obtained as the sum of kinetic energies; and 4) torques Ti at each DOF were calculated using
which leads to
 | (8) |
Coefficients Aij and Ci were obtained with a tool for symbolic calculus. Equation 8 was then inverted to obtain relationships between angular accelerations and torques, which were inserted into Eq. 5. For illustration, the resulting equations are given below for N = 2. The other cases involve lengthy equations that cannot be reproduced here (e.g., for N = 7, Eq. 8 requires more than 200,000 elementary operations [+, –, x, cos, sin]).
View larger version (10K):
[in this window]
[in a new window]
|
FIG. 3. Kinematic chains. A: 4 degrees of freedom (4-DOF), 2-link chain in 3-dimensional (3D) space. B: 7-DOF, 3-link chain in 3D space. C: 2-DOF, 2-link chain in 2D space.
|
|
4-DOF, 2-LINK CHAIN IN 3D SPACE (FIG. 3A).
In this case, N = 4. The effector had two segments (upper arm, forearm) and two joints (shoulder, elbow) with 3 DOF at the shoulder and 1 DOF at the elbow. Arm kinematics was represented by
where (r1, r2, r3) are the endpoint coordinates; (q1, q2, q3, q4) are the shoulder azimuth angle, the shoulder elevation angle, the shoulder intrinsic (humeral) rotation angle, and the elbow rotation angle, respectively; and L1, L2 are the upper arm and forearm lengths. As a result of redundancy, final arm configuration was specified indirectly by the desired endpoint coordinates in function (Eq. 6). Initial arm configuration was (q1, q2, q3, q4) = (120, 140, 30, 90°).
7-DOF, 3-LINK CHAIN IN 3D SPACE (FIG. 3B).
In this case, N = 7. The effector had three segments (upper arm, forearm, hand) and three joints (shoulder, elbow, wrist) with 3 DOF at the shoulder, 2 DOF at the elbow, and 2 DOF at the wrist. Arm kinematics was represented by
where (r1, r2, r3) are the endpoint coordinates; (q1, q2, q3, q4, q5, q6, q7) are the shoulder azimuth angle, the shoulder elevation angle, the shoulder intrinsic (humeral) rotation angle, the elbow rotation angle, the elbow intrinsic rotation angle, the wrist flexion/extension angle, and the wrist pronation/supination angle, respectively; and L1, L2, L3 are the upper-arm, forearm, and wrist lengths. Initial arm configuration was (q1, q2, q3, q4, q5, q6, q7) = (120, 120, 90, 90, 60, 90, 10°). As a result of redundancy, final arm configuration was specified indirectly by the desired endpoint coordinates and desired hand orientation in function (Eq. 6). In fact, the target was an elongated object and the constraints were that 1) the endpoint of the arm matches the center of gravity of the object and 2) the hand (i.e., the vector in the plane of the hand and perpendicular to its long axis) is oriented parallel to the long axis of the object.
2-DOF, 2-LINK CHAIN IN 2D SPACE (FIG. 3C).
In this case, N = 2. The controlled system had two segments (upper arm, forearm) and two joints (shoulder, elbow) with 1 DOF per segment. Kinematics was represented by
where (r1, r2) are the endpoint coordinates; (q1, q2) represent the shoulder and elbow angles, respectively; and L1, L2 are the upper arm and forearm lengths. In general, initial arm configuration was (q10, q20) = (45, 90°). Desired final configuration was (q1f, q2f). The state vector was expressed by
where v1(t) = dq1/dt and v2(t) = dq2/dt (angular velocities). The dynamics of the arm + muscles was expressed by the following series of equations
where
 | (9) |
are muscular shoulder and elbow torques, respectively, and
with Lc1 = c1L1 and Lc2 = c2L2 (distance from rotation axis to the center of gravity of the segments), and where m1, m2 represent the mass of the segments, and I1, I2 are the moments of inertia of the segments. The boundary conditions were
and
(Velocities and forces were zero at the beginning and end of the movement.) More generally, any boundary conditions can be used (nonzero velocities or forces).
Numerical methods
The problem defined by Eqs. 5, 6, and 7 is an optimization problem that can be transformed into a differential equation problem using the calculus of variations (Kirk 1970). The differential equation problem is a two-point boundary-value problem, i.e., a differential system with boundary conditions at initial and final times. Solutions to this problem were obtained numerically using a gradient method (Bryson 1999). Function (z) = [z]+ is not differentiable and was replaced by the differentiable function z 
log [1.0 + exp(
z)]/.
Parameters
General parameters were = 0.04 s (time constant of muscle filtering) and = 10 (characteristic of muscle force generation). Other parameters were specific to each kinematic chain. Each segment i (1: upper arm; 2: forearm; 3: hand) is defined by a mass (mi in kilograms), a length (Li in meters), a center of gravity (ci in percentage of the length), and moments of inertia along and perpendicular to its main axis (Ji and Ii in kilograms per square meter). For each DOF, there is a coefficient (in meters) that translates force into torque.
For N = 2, parameters were m1 = 2.52, L1 = 0.33, c1 = 0.5, I1 = 0.023, m2 = 1.3, L2 = 0.4, c2 = 0.5, I2 = 0.011, and sh = el = 0.04. The choice of sh and el is somewhat arbitrary. In fact, these parameters result from the interplay between: 1) the innervation ratio of the muscles acting at shoulder and elbow; 2) the moment arm of these muscles; and 3) the modulation of force production by firing rate and recruitment in pools of motoneurons. Rather than doing a complex estimation based on the contribution of these factors, we explored the influence of sh and el on the behavior of the model. The model was mostly insensitive to the values of sh and el as long as the ratio el/sh was neither too small (>0.05) nor too large (<1.2).
For N = 4, parameters were as defined for N = 2 and J1 = 0.0013, 1 = 10, 2 = 10, 3 = 10, and 4 = 5.
For N = 7, parameters were as defined for N = 4, but m2 = 1.3, L2 = 0.25, I2 = 0.0074 (shorter forearm), and m3 = 0.49, L3 = 0.15, c3 = 0.25, I3 = 0.001, J3 = 0.0005, 5 = 1, 6 = 1, and 7 = 1.
When a mass ma was added to a segment (length L, inertia I, mass m, center of mass c) at position La, the new inertia of the segment was
 | (10) |
Comparison with experimental data
The model is a parameter-free model, i.e., all the parameters are well defined and constant for a given individual at a given time (e.g., forearm inertia). As a consequence, we did not try to provide the best fits to the experimental data because 1) fitting is general associated with arbitrary parameter adjustments and 2) fitting quality may not be sufficient to estimate the validity of a model (e.g., see the debate between minimum-jerk and minimum torque change models of motor control; Flash 1990). The results reported here are meant to demonstrate the power of the concepts underlying the model, rather than account for a few selected aspects of the behavior. When possible, we indicate in the text a reference to one or more published figures that can be used for comparisons with the model.
|
RESULTS
|
|---|
We describe implications of the proposed principles for point-to-point movements. Except for cases dealing with kinematic redundancy, the 2-DOF arm was used. Except for cases illustrating the constant effort principle, the simulated movements were specified by their amplitude and duration.
Kinematic redundancy
3D POINTING MOVEMENTS.
We simulated optimal point-to-point trajectories of a 4-DOF arm in 3D space (Fig. 3A). An example is illustrated in Fig. 4A for a forward/upward/leftward movement. The trajectory was curved, planar, independent of speed, and the tangential velocity had a bell-shaped profile (Fig. 4A). Further examples of trajectories and velocity profiles are shown for movements in a frontal (Fig. 4, B and C) and sagittal (Fig. 4, D and E) plane (for comparison, see Fig. 3 in Flanders et al. 1996; Fig. 1 in Gottlieb et al. 1997). Hand path curvature varied with movement direction and went through one cycle (Fig. 4F). We note that quantitative data on movement curvature are scarce. The most complete study is for vertical movements in a sagittal plane (Flanders et al. 1996; Pellegrini and Flanders 1996). Curvature could be described by: curvature = cos (direction + phase), which also holds for the model (for comparison, see Fig. 2 in Pellegrini and Flanders 1996). However, we cannot expect to obtain an exact fit because the phase should depend on arm posture and the hypothesized mechanical actions of the muscles (in the data of Pellegrini and Flanders 1996, the phase varied across subjects). Soechting and Flanders (1998) showed that neither minimum-torque change nor minimum-muscle force change models can explain the pattern of curvature for vertical movements. The terminal posture was independent of movement velocity (Fig. 4, G and H). These observations were qualitatively similar for all tested movements and were consistent with experimental observations (Alexandrov et al. 1998; Atkeson and Hollerbach 1985; Flanders et al. 1999; Hermens and Gielen 2004; Klein Breteler et al. 1998; Nishikawa et al. 1999; Pellegrini and Flanders 1996; Soechting and Lacquaniti 1981; Torres and Zipser 2004; Zhang and Chaffin 1999). The terminal posture was also independent of forearm inertia (<2 deg of variations for movements in Fig. 4, B and D) over the tested range (using a weight of 0.6–1.6 kg attached to the forearm at 3L2/4 from the elbow; Hermens and Gielen 2004). We verified that the terminal posture depended on initial posture (Gielen et al. 1997; Hermens and Gielen 2004; Soechting et al. 1995). However, a quantitative study was not pursued on this point because these studies did not report actual initial postures of the arm that were used. Initial hand positions are not sufficient in this case. A quantitative analysis would require establishing the relationship between initial and final arm angles across a large range of movements (different initial postures, different directions). This issue was addressed in the context of grasping movements (see following text; Fig. 7).
View larger version (27K):
[in this window]
[in a new window]
|
FIG. 4. Optimal control for a redundant arm. A: 2 views of a simulated trajectory for a forward/leftward/upward movement (20 cm, 600 ms). Inset: tangential velocity profile. B: movements in 8 directions (15 cm, 600 ms) in the frontal plane projected on this plane (left), on a sagittal plane (right), and on a dorsal plane (bottom). U, up; F, forward; D, down; B, backward. C: tangential velocity profiles for the movements in B. D: movements in 8 directions in a sagittal plane projected on this plane (left), on a frontal plane (right), and on a dorsal plane (bottom). R, right; L, left. E: tangential velocity profiles for the movements in D. F: hand path curvature (in millimeters) measured as the mean deviation from straight line (counterclockwise deviation >0) for data in B (left) and D (right). G: final angular positions (circle: q1; box: q2; down triangle: q3; up triangle: q4) for the movements in B at 3 velocities (600, 800, 1,000 ms). Offsets have been used to reveal superimposed points. H: same as G for data in D.
|
|
View larger version (8K):
[in this window]
[in a new window]
|
FIG. 5. Upward/downward movements. Simulation of the 4-DOF arm (same model as in Fig. 4) for 400- to 1,200-ms, 30-cm upward and downward movements. Initial posture for upward movement was as in Fig. 4. Initial posture for downward movements was the final posture of upward movements. A: trajectories (black: upward; gray: downward) projected on a frontal plane (as viewed when facing the subject). B: normalized velocity profiles for upward (black) and downward (gray) movements (solid: 600 ms; dashed: 400 ms).
|
|
View larger version (9K):
[in this window]
[in a new window]
|
FIG. 6. Grasping movements with a 7-DOF arm. A: sample trajectory and velocity profile for a 20-cm, 600-ms movement (270 deg = toward the body). Target is indicated by a vertical cylinder. B: normalized time course of orientation error vs. time course of target distance error along the movement (coarticulation) for the movement in A. C: variations of hand azimuth with hand-centered movement direction (20-cm, 600-ms movements in the horizontal plane toward a vertical cylinder; movement direction is 0 deg rightward).
|
|
View larger version (11K):
[in this window]
[in a new window]
|
FIG. 7. Variations of final posture as the function of movement starting point (high, open circle; down, closed circle) and object position (Sa, sagittal; Cl, close; La, lateral) for a 7-DOF arm. Initial posture was (100, 105, 115, 80, 60, 90, 10) for high and (80, 160, 60, 80, 60, 90, 10) for down. Movement amplitude (in centimeters) and final hand orientation (two angles, in deg) were 53.5, 105.4, 74.2 (downSa), 52, –29.3, 82.2 (downCl), 55.8, 41.1, 67.4 (downLa), 31.3, 176.9, –11.1 (highSa), 26.7, –142.4, –13.0 (highCl), 12.3, 171.3, –29.3 (highLa). Movement duration was 600 ms.
|
|
According to the separation principle, trajectories should not be influenced by gravity, which is at variance with experimental observations (Atkeson and Hollerbach 1985; Papaxanthis et al. 2003). However, results on up–down and down–up movements are ambiguous because movement trajectories may depend not only on gravity, but also on initial and terminal positions. Thus direction-dependent hand paths can be expected in an optimal control model in the absence of gravity. This was found in the model (Fig. 5). In particular, we observed differences between upward and downward trajectories (Fig. 5A) and between the corresponding velocity profiles (Fig. 5B) (for a comparison, see Fig. 2 in Papaxanthis et al. 2003). These results indicate that the observed differences between upward and downward movements are not incompatible with the separation principle. More behavioral data, including systematic variations of movement kinematics with initial and final position, would be necessary to distinguish between gravity and position effects.
Grasping movements
An additional source of redundancy arises for the control of a distal segment, such as a hand that grasps an object. Additional degrees of freedom are related to hand orientation and aperture and additional constraints are provided by the shape of the object. We addressed the control of hand orientation with a 7-DOF arm (Fig. 3B). We simulated movements toward an elongated object at different locations (see METHODS).
The trajectories were similar to those observed with a 4-DOF arm (Fig. 6A). The model exhibited coarticulation of hand transport and rotation along the path (Fig. 6B). There was a linear relationship between movement direction and hand azimuth (Fig. 6C). These results are consistent with experimental observations (Bennis and Roby-Brami 2002; Cuijpers et al. 2004; Desmurget et al. 1996; Roby-Brami et al. 2000; Torres and Zipser 2004).
To assess how final posture depends on initial posture, we reproduced the experiment of Desmurget et al. (1998). Movements from two initial postures (high and low) toward three targets (sagittal, close, lateral) were simulated. Final shoulder and elbow angles were generally different for movements toward the same target, but starting from different postures (Fig. 7). The variations were quantitatively similar to those observed experimentally (see Fig. 2 in Desmurget et al. 1998). A difference was found with respect to the variations of forearm rotation angle that could be explained by differences in initial hand orientation. This result was not self-evident a priori. Indeed, not all optimal control models are able to explain the effects of initial posture on final posture (Admiraal et al. 2004).
Kinematic invariance
The proposed principles also have consequences that are not directly related to redundancy. These consequences were explored for a 2-DOF arm (Fig. 3C). When subjects perform movements of different amplitudes or against different inertial loads, they tend to use a single velocity profile that is scaled in time and amplitude (Atkeson and Hollerbach 1985; Bock 1990; Gordon et al. 1994; Kaminski and Gentile 1989). Invariant velocity profiles were observed as a consequence of the constant effort principle, i.e., for movements of different amplitudes when movement time was chosen to obtain a given effort level (Fig. 8A; for comparison see Fig. 3 in Gordon et al. 1994). In fact, the assumption that all movements are executed with the same effort leads to an implicit relationship between amplitude and duration, illustrated in Fig. 8B. For the eight directions tested, movement duration and peak velocity were linearly related to amplitude (Fig. 8, C and D; for comparison, see Fig. 4 in Gordon et al. 1994). The separation principle predicts that amplitude/duration scaling should be similar for upward and downward movements, previously observed experimentally (Virji-Babul et al. 1994). Moreover, movement duration and peak velocity varied with movement direction. As observed experimentally (Gordon et al. 1994; Turner et al. 1995), the duration was longer and peak velocity smaller for directions parallel to the forearm (Fig. 8, E and F) in which the inertial load to be moved was larger (for comparison, see Fig. 6 in Gordon et al. 1994; Fig. 9 in Sober and Sabes 2003; see also Pellegrini and Flanders 1996, Fig. 3 and for vertical movements, Fig. 5 in Murata and Iwase 2001; see also Fernandez and Bootsma 2004; Jagacinski and Monk 1985; Smyrnis et al. 2000). The shape of velocity profiles was described by the ratio c between peak and average velocity (Ostry et al. 1987). This factor varied with the level of effort (Fig. 8G, inset), but was almost constant across amplitudes (Fig. 8G). Exceptions were observed in directions of larger curvature (Fig. 8, E and G). We tested the hypothesis that kinematic invariance is actually related to movement curvature. We simulated straight movements in the 135° direction by forcing the curvature to be as small as possible. We determined the amplitude/duration relationship and found that it was still linear and the velocity profile was strictly invariant. It should be noted that, in fact, kinematic invariance is an emergent property, given that no desired trajectory or velocity profile was specified for the simulated movements.
View larger version (25K):
[in this window]
[in a new window]
|
FIG. 8. A: kinematic invariance under constant effort (E2 = 10) for a planar arm. Actual (left) and normalized (right) velocity profiles for 10, 15, and 20 cm rightward movements. B: movement duration vs. amplitude for rightward movements (circles correspond to profiles in A). Inset: correspondence between colors and directions; thick (thin) black line: minimum (maximum) inertia. 0° is rightward. C: movement duration vs. amplitude for 8 directions. D: peak velocity vs. amplitude. E: plain line: duration vs. direction for a given amplitude (15 cm). Dotted line: curvature vs. direction. Curvature was measured as the mean deviation (in millimeters) from straight line (positive for counterclockwise deviation). F: peak velocity vs. direction for a given amplitude (15 cm). G: c = vpeak/vavg vs. amplitude. Inset: distribution of c for 100 movements (5–20 cm, 200–500 ms, direction 0°). See METHODS for construction of this figure.
|
|
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 9. A: kinematic invariance for added mass under constant effort (E2 = 10) for a planar arm. Actual (left) and normalized (right) velocity profiles for a 15-cm rightward movement with 0, 1, and 2 kg added mass (m). Mass was a punctual mass at the arm endpoint. B: scaling factor in time: k(m) = MT(m)/MT(m = 0), where MT is movement time. C: peak velocity vs. added mass (black line). Relation predicted by vpeak(m) = vpeak(m = 0)/k(m) (gray line). D, top: normalized shoulder EMG for the 3 movements in A (flexor: black; extensor: gray). Bottom: normalized elbow EMG for the 3 movements in A.
|
|
Constant effort also predicts invariance for inertial loads (Fig. 9A). These results were obtained by 1) adding a mass at the arm endpoint and 2) calculating optimal control solution for the system arm + mass (new inertia was calculated with Eq. 10). In this case, larger loads lead to slower movements. Interestingly, the timescaling factor (Fig. 9B) was very close to the amplitude-scaling factor (Fig. 9C) as observed experimentally (Figs. 2, 3, and 4 in Bock 1990; see also Hatzitaki and McKinley 2001). We note that scaling in time and amplitude occurred in the model in the absence of scaling at the level of joint torques, EMGs (Fig. 9D), and control signals (Gentner 1987; Zelaznik et al. 1986).
Muscular redundancy
Motor control systems master not only kinematic redundancy, but also muscular redundancy, e.g., because of the presence of biarticular muscles. An open question is how muscular redundancy would affect the preceding results obtained with a minimal agonist/antagonist organization. We first note that optimal control can actually tackle muscular redundancy (e.g., Dornay et al. 1996). Two biarticular muscles were introduced in the 2-DOF model: a shoulder/elbow flexor and a shoulder/elbow extensor. Joint torques were (see Eq. 9)
where a5 and a6 are excitations for the biarticular flexor and extensor, respectively, and the 
parameter specifies the contribution of the biarticular muscles relative to the monoarticular muscles. We chose = 1. On the one hand, the presence of biarticular muscles had little influence on movement kinematics (Fig. 10, A and B). On the other hand, the contribution of monarticular shoulder and elbow muscles was dramatically modified (Fig. 10, C and D). Part of their original contribution was now subserved by the biarticular muscles (Fig. 10E). These observations were graded with respect to . These results indicate that the covariations reported in preceding sections should be qualitatively immune to the presence of muscular redundancy. This does not mean that biarticular muscles are useless. In fact, from the perspective of a controller, mono- and biarticular muscles are not different—that is, they are exploited at best to satisfy the required constraints.
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 10. Influence of biarticular muscles for movements of a planar 2-link arm. A: simulation of a 30-cm, leftward, 500-ms movement (inset). Normalized velocity profiles are shown. Plain line: model with biarticular muscles. Dashed line: model without biarticular muscles. B: angular trajectories (thick: shoulder angle; thin: elbow angle). C: EMG in shoulder flexor (thick) and shoulder extensor (thin). D: EMG in elbow flexor (thick) and elbow extensor (thin). Elbow flexor is agonist in the absence of biarticular muscles and antagonist in the presence of biarticular muscles. E: EMG in biarticular flexor (thick) and extensor (thin) muscles.
|
|
Influence of muscle properties
The preceding results were obtained using muscles modeled as force-generating elements. We assessed the influence of the force/length and force/velocity (called viscosity for simplicity) relationship in muscles (Hill model). Briefly, the force F generated by a muscle (Eq. 4) was used to calculate the maximal isometric force
where 
is muscle length and K is a parameter. The actual velocity-modulated muscle force was
where d/dt is muscle velocity, a' = –0.4Fiso, b' = –b(a' + Fiso)/(a + Fiso), a = 250, b = 0.5, K = 0.1 (Fig. 11D, inset). For simplicity, we assumed that muscle length was proportional to joint angle. As shown in Fig. 11, viscosity had a dramatic effect on EMGs (Fig. 11, C and D), but little influence on movement kinematics (Fig. 11, A and B). This result was confirmed over a broad range of movements. We did not find any remarkable properties related specifically to the presence of viscosity, i.e., properties that would be absent in the absence of viscosity. In fact, optimal control builds an internal model of the neuromuscular system. The muscular properties are exploited by the controller to reach its goal (for related ideas, see Todorov 2000). For instance, if a muscle generates less force because it is shortening, the controller could increase the control (change in EMG) or, if it is too costly, use another DOF (change in kinematics).
View larger version (13K):
[in this window]
[in a new window]
|
FIG. 11. Influence of force/velocity relationship in muscles for movements of a planar 2-link arm. A: simulation of a 30-cm, leftward, 500-ms movement (inset). Normalized velocity profiles are shown. Plain line: model with nonlinear muscles. Dashed line: model with linear muscles. B: angular trajectories (thick: shoulder angle; thin: elbow angle). C: EMG in shoulder flexor (thick) and shoulder extensor (thin). D: EMG in elbow flexor (thick) and elbow extensor (thin). Inset: velocity-modulated muscle force (Fvel) as a function of the rate of change of muscle length (dl/dt) for different levels of maximal isometric for (Fiso).
|
|
Muscular synergies
We assumed that 1) the number n of control signals was larger than the number of muscles; 2) each control was defined by a fixed synergy of muscles; and 3) the s synergies were uniformly distributed in muscular space. Formally, the problem was similar to the problem defined by Eqs. 5, 6, and 7 with the following change. The goal was to find minimum controls U(t) = {Uj(t)} (1 j s) such that the controls {ui(t)} (1 i 2N), defined by
where 
ij are random coefficients drawn from a uniform distribution in [–1; 1] are appropriate to displace the articulated segments between given initial and final positions. Simulations did not reveal any qualitative differences between synergistic (s = 500) and direct (s = 2N) control.
|
DISCUSSION
|
|---|
A proper solution to the degrees of freedom problem should be able to explain 1) how a unique solution is chosen for each realization of a motor act and 2) why this solution is different each time. The central idea implicit in the proposal of Bernstein and explicitly formulated by Todorov and Jordan (2002) is that the nervous system continuously and efficiently tracks a goal rather than a desired trajectory. This idea was expressed here by two principles (optimal feedback control and maximum efficiency). We have shown that these two principles combined with a separation principle for static and dynamic forces lead to a realistic solution to redundancy. A fourth principle was added (constant effort), which is not directly related to redundancy, to obt
TOP
ABSTRACT
INTRODUCTION
