The Journal of Neurophysiology Vol. 79 No. 3 March 1998, pp. 1409-1424
Copyright ©1998 by the American Physiological Society

Are Complex Control Signals Required for Human Arm Movement?

Paul L. Gribble¹, David J. Ostry¹, Vittorio Sanguineti², and Rafael Laboissière³

¹ Department of Psychology, McGill University, Montreal, Quebec H3A 1B1, Canada; ² Department of Informatics, Systems and Telecommunications, University of Genova, Genova, Italy 16145; and ³ Institut de la Communication Parlée, Grenoble, France 38031

	ABSTRACT

Abstract Introduction Discussion References

Gribble, Paul L., David J. Ostry, Vittorio Sanguineti, and Rafael Laboissière. Are complex control signals required for human arm movement? J. Neurophysiol. 79: 1409-1424, 1998. It has been proposed that the control signals underlying voluntary human arm movement have a "complex" nonmonotonic time-varying form, and a number of empirical findings have been offered in support of this idea. In this paper, we address three such findings using a model of two-joint arm motion based on the  version of the equilibrium-point hypothesis. The model includes six one- and two-joint muscles, reflexes, modeled control signals, muscle properties, and limb dynamics. First, we address the claim that "complex" equilibrium trajectories are required to account for nonmonotonic joint impedance patterns observed during multijoint movement. Using constant-rate shifts in the neurally specified equilibrium of the limb and constant cocontraction commands, we obtain patterns of predicted joint stiffness during simulated multijoint movements that match the nonmonotonic patterns reported empirically. We then use the algorithm proposed by Gomi and Kawato to compute a hypothetical equilibrium trajectory from simulated stiffness, viscosity, and limb kinematics. Like that reported by Gomi and Kawato, the resulting trajectory was nonmonotonic, first leading then lagging the position of the limb. Second, we address the claim that high levels of stiffness are required to generate rapid single-joint movements when simple equilibrium shifts are used. We compare empirical measurements of stiffness during rapid single-joint movements with the predicted stiffness of movements generated using constant-rate equilibrium shifts and constant cocontraction commands. Single-joint movements are simulated at a number of speeds, and the procedure used by Bennett to estimate stiffness is followed. We show that when the magnitude of the cocontraction command is scaled in proportion to movement speed, simulated joint stiffness varies with movement speed in a manner comparable with that reported by Bennett. Third, we address the related claim that nonmonotonic equilibrium shifts are required to generate rapid single-joint movements. Using constant-rate equilibrium shifts and constant cocontraction commands, rapid single-joint movements are simulated in the presence of external torques. We use the procedure reported by Latash and Gottlieb to compute hypothetical equilibrium trajectories from simulated torque and angle measurements during movement. As in Latash and Gottlieb, a nonmonotonic function is obtained even though the control signals used in the simulations are constant-rate changes in the equilibrium position of the limb. Differences between the "simple" equilibrium trajectory proposed in the present paper and those that are derived from the procedures used by Gomi and Kawato and Latash and Gottlieb arise from their use of simplified models of force generation.

	INTRODUCTION

Abstract Introduction Discussion References

According to the equilibrium point hypothesis, voluntary movements arise as a consequence of shifts in the equilibrium of the motor system. This equilibrium is dependent upon the interaction of control signals, reflexes, muscle properties, and loads, but it is under the control of central commands. In the context of this hypothesis, it has been proposed that "simple" monotonic equilibrium shifts may be used to produce smooth arm movements (Feldman 1986; Feldman et al. 1990; Flanagan et al. 1993) (Note that we use the term equilibrium shift to denote the changes to neural control signals to muscles that result in changes to the physical equilibrium of the limb.) Empirical studies in support of this view suggest that the equilibrium shift is gradual (Bizzi et al. 1984), that it is similar in form to the actual movement (Won and Hogan 1995), and that it ends substantially before the end of the movement (Feldman et al. 1995).

Alternatively, it has been claimed that the equilibrium shifts that underlie human arm movements have a nonmonotonic time-varying form and that these "complex" control signals are needed both to generate torques that are large enough to produce fast movements (Latash and Gottlieb 1991) and to compensate for the dynamics of the multijoint arm (Gomi and Kawato 1996; Hogan 1984). The recent empirical demonstration that limb impedance patterns during multijoint movement are nonmonotonic and have multiple peaks has been offered in support of the idea that control signals must be complex to compensate for dynamics.

A number of recent studies have explored the complexity of control signals for arm movements. Gomi and Kawato (1996) measured arm stiffness just before and during planar two joint movements. They report that joint stiffness is low before movement (<20 N·m/rad for the shoulder), increases during movement (to a peak value for the shoulder in the range of 40 N·m/rad), and varies over the course of movement in a nonmonotonic fashion. Using these empirically derived stiffness measures and measures for viscosity, a hypothetical equilibrium trajectory with a complex time-varying form was computed. These computations were based on the assumption that joint torques vary linearly with the difference between actual and equilibrium position and with velocity. A comparable procedure has been used to infer the form of control signals during rapid single-joint movements (Latash and Gottlieb 1991). In that study, subjects produced elbow flexion movements during ramp changes in external torque. By applying torques of different magnitudes in different directions and by assuming that joint torque varied linearly with joint angle (but was not dependent on velocity), it was possible to infer a nonmonotonic equilibrium joint angle over the course of the movement.

These claims concerning the complexity of control signals have been made in the context of highly simplified models of force generation that lack explicit representations of individual muscles and a variety of neurophysiological properties (see DISCUSSION). In the present paper, we use a model of human arm movement based on the  version of the equilibrium point hypothesis to reexamine these claims. The model is used to simulate studies by Gomi and Kawato (1996), Bennett (1993), and Latash and Gottlieb (1991). We show that the empirical patterns of kinematics and joint stiffness that have been offered as evidence for the complexity of control may be predicted using simple control signals involving constant rate changes in the limb equilibrium position. Our results underscore the need for explicit models of muscle properties, musculoskeletal geometry, and limb dynamics when testing ideas concerning the complexity of neural control (also see Gribble and Ostry 1996; Laboissière et al. 1996; Ostry et al. 1996).

	ARM MODEL

Earlier versions of the arm model have been described in Feldman et al. (1990), Flanagan et al. (1993), and Gribble and Ostry (1996). The arm model has two kinematic degrees of freedom: rotation at the elbow and at the shoulder, in the horizontal plane. Six muscle groups are modeled: single-joint flexors and extensors at the shoulder (pectoralis and deltoid) and elbow (biceps long head and triceps lateral head) and double-joint muscles spanning both joints (biceps short head and triceps long head). Muscle origins and insertions are estimated from anatomic sources (An et al. 1981, 1989; Winters and Woo 1990). Muscle moment arms for the three extensor muscles are assumed to be constant (2 cm at the elbow and 4 cm at the shoulder). Flexor moment arms vary with joint angle and are calculated on the basis of musculoskeletal geometryvalues range from 2.5 to 5 cm. The equations of motion relating joint torques to accelerations were obtained using the Lagrangian approach (see Hollerbach and Flash 1982). The inertial and geometrical constants are upper and lower arm mass: 2.1 and 1.65 kg, length: 0.34 and 0.46 m, and moment of inertia about the center of mass: 0.015 and 0.022 kgm².

For each muscle, we model the dependence of force on muscle length, on the velocity of muscle lengthening or shortening, the graded development of force over time, and the passive elastic stiffness of muscle (see Fig. 1). The model is a variant of that described by Zajac (1989), with activation and contraction dynamics and passive muscle stiffness. For each muscle we also include modeled neural inputs, length- and velocity-dependent afferent feedback, and reflex delays. The values for all model parameters except those related to the time-varying form of the neural control signal were estimated from empirical data. We used representative examples from the physiological literature and then tested the sensitivity of predicted outcomes to model parameter values (see APPENDIX C).

View larger version (13K): [in this window] [in a new window]   FIG. 1. Schematic representation of the muscle model. Model includes the dependence of force on muscle length and velocity, graded force development over time, passive elastic stiffness, and length and velocity dependent afferent input.

Control signals in the model are based on the  version of the equilibrium-point hypothesis (Feldman 1986; Feldman et al. 1990). According to the  model, neural control signals establish muscle threshold lengths (s) for  motoneuron (MN) recruitment. By changing the value of , a new threshold length is specified and force is generated in proportion to the difference between a muscle's actual length and , and in proportion to the velocity of lengthening or shortening. By setting s for all muscles, a static equilibrium configuration of the limb can be specified. By changing the values of all s in appropriate proportions (see further for details), a movement can be generated from one static equilibrium configuration to another. In this paper, it is assumed that movements are generated by "shifting," under neural control, the equilibrium position of the hand in a straight line at a constant rate from an initial equilibrium configuration to a final limb configuration.

Muscle model

Each of the six muscles is modeled separately and includes neural input, . Muscle forces are generated in the following way. Muscle activation, A (which corresponds to the effects of neural excitation on the contractile machinery of the muscle), is proportional to the difference between the current muscle length l, and the centrally specified threshold length for MN recruitment, , as well as on the rate of change of muscle length, Thus

(1)

where

(2)

The muscle lengths and velocities associated with positive values of A thus define an "activation area"a region in state-space in which muscle activity is observed (Feldman et al. 1990; Levin and Feldman 1994).

The parameter µ specifies the dependence of the muscle's threshold length on velocity and provides damping due to proprioceptive feedback. In the simulations presented here, we have assumed that µ is 0.06 s and is constant over time. The parameter µ likewise reflects the relative magnitudes of spindle primary to spindle secondary afferent effects on motoneuron recruitment (see DISCUSSION). Because primaries and secondaries may be affected differentially by fusimotor inputs to intrafusal muscle fibers, then in principle, µ may be modulated through central commands (Feldman et al. 1990). In the present paper, the magnitude of µ was estimated by matching the predicted viscosity of the model in statics to corresponding empirical estimates of viscosity (see APPENDIX A for details and APPENDIX C for sensitivity analysis).

A reflex delay, d, of 25 ms has been used for all muscles. This value was estimated from observed delays in the unloading response of human arm muscles (Houk and Rymer 1981). Thus taking into account time-varying central neural commands (t) and a reflex delay, d, muscle activation A(t) is

(3)

Changes to  and thus to muscle activation are associated with MN recruitment and changes in MN firing rates. The resulting muscle force, (representing force due to the contractile apparatus), is approximated using an exponential function of the form

(4)

where c is a form parameter and is the same for all muscles, and  is a magnitude parameter related to force-generating capability and is specific to each muscle. Note that this relationship is consistent with the size principle, where c is associated with the MN recruitment gradient. As the difference between actual and threshold muscle length  increases, progressively larger motor units are recruited and larger increments in force are observed. The value of the c parameter was estimated from empirical force-length data for cat gastrocnemius muscle (Feldman and Orlovsky 1972). (These data were recorded in a preparation with intact dorsal and ventral roots and with activation produced in a physiological manner by means of stimulation delivered at the level of the brain stem.) To estimate c, passive stiffness, which was assumed to be linear (see DISCUSSION), first was subtracted from total force. Then, a regression technique was used to approximate the set of force-length relations reported by Feldman and Orlovsky (1972) with a single exponential function, subject to the constraint that the magnitude and form parameters,  and c, remained constant, and a single parameter, , was free to vary. With this method, we obtained a value of c = 0.112 mm¹, which is similar to that obtained for human elbow muscles (Feldman 1966).

The  parameter is associated with the force-generating abilities of individual muscles. We have assumed that values of vary in proportion to the modeled muscle's physiological cross-sectional area (PCSA). The following estimates of PCSA, obtained from Winters and Woo (1990), were used: biceps short head, 2.1 cm²; biceps long head, 11 cm²; deltoid, 14.9 cm²; pectoralis, 14.9 cm²; triceps lateral head, 12.1 cm²; and triceps long head, 6.7 cm². PCSA estimates reflect relative force-generating abilities. To convert values of PCSA into estimates of muscle force-generating ability that match empirical data, values of  were computed by scaling PCSA values by 1 N/cm². This was done to match the predicted stiffness of the model in statics to corresponding empirical measures of static joint stiffness (see APPENDIXES A and C for procedure and further details). Note that the procedures that establish values for both µ and  are based on achieving a correspondence between model behavior and empirical data in statics only. The overall time-varying form of stiffness and viscosity during movement is not dependent upon this procedure.

The graded development of muscle force due to calcium kinetics is modeled using a second-order, low-pass filtering of the muscle force, , where M represents the instantaneous value of muscle force

(5)

The filter is critically damped with a single parameter, . A value of  = 15 ms was chosen; this leads to an asymptotic response to a step input in ~90 ms. This corresponds to empirically observed times from onset of stimulus to maximum force in human adductor pollicis muscle (Hainaut et al. 1981).

The dependence of muscle force on the velocity of muscle lengthening or shortening was estimated from data for cat soleus muscle at different rates of stimulation of the motor nerve (Joyce and Rack 1969). Force levels for different stimulation rates were normalized before fitting (velocities were not normalized) (see Zajac 1989 for discussion). The data were fit to a sigmoidal function of the form

(6)

where F is the resulting force, and coefficients f1-f4 have values of 0.82, 0.50, 0.43, and 58 s/m, respectively. The final force, M, in Eq. 6 is given as the sum of force F (generated by MN recruitment), and passive force (muscle force in the absence of neural input). As a simplification, we have assumed that passive force is linearly dependent on the difference between current muscle length, l, and muscle resting length, r (see DISCUSSION). Values of r were computed as the muscle lengths associated with a shoulder angle of 45° and an elbow angle of 90°. The passive stiffness of muscles in the arm model was assumed to vary linearly with physiological cross-sectional area and was scaled to match the passive component of the force-length relation shown in Feldman and Orlovsky (1972). The following values of passive stiffness are used: biceps short head, 36.5 N/m; biceps long head, 190.9 N/m; deltoid, 258.5 N/m; pectoralis, 258.5 N/m; triceps lateral head, 209.9 N/m; and triceps long head, 116.3 N/m.

Organization of control signals

Following earlier work with the model, we describe how control signals are organized to provide two types of commandsone that shifts muscle threshold lengths of agonist and antagonist muscles in opposite directions, generating movement between equilibrium positions and another that shifts threshold lengths in the same direction, independently modulating the level of muscle coactivation.

In the simulations presented in this paper, we have generated movements using a series of equally spaced equilibrium positions or via-points to produce straight-line equilibrium shifts in hand space. Because the number of muscles (6) exceeds the kinematic degrees of freedom of the model (2), there are an infinite number of sets of s associated with a given via-point, each set corresponding to different levels of muscle force. For each via-point, there exists a set of s that minimizes total muscle force. In the present model, movements are produced by shifting s at a constant rate between these points of minimum force. This is analogous to the "R" command in previous versions of the model (Feldman et al. 1990; Flanagan et al. 1990).

In addition, we define a cocontraction command that shifts all s in the same direction. In statics, this increases muscle forces but produces no net change in joint torques (and hence no movement). In dynamics, the application of the cocontraction command increases the size of the activation area for agonist and antagonist muscles (see previous text) and, as a result, muscles reach their threshold lengths earlier during movement. This is analogous to the "C" command described in previous versions of the model (Feldman et al. 1990; Flanagan et al. 1990). For a given equilibrium position, there are an infinite number of possible cocontraction commandswe have chosen one that increases muscle forces in the most equal proportions without changing net joint torque. The  shifts associated with the cocontraction command can be scaled in magnitude and applied in combination with the  shifts, which yield movement between equilibrium positions. In all simulations presented here, the magnitude of the cocontraction command was constant throughout movement. It should be noted that for purposes of these simulations, the cocontraction command initially was defined in force space, and hence the units of the cocontraction command are N. For a cocontraction command of a given force level in N, a vector in force space (a set of muscle forces) associated with that average increase in muscle force and zero net joint torque was found. The vector in  space associated with this change in muscle force in statics then was computed and served as the cocontraction command.

Parameters of the model have been obtained either directly from the physiological literature, or as in the case of  and µ, by matching the model's performance to empirical data in statics (see APPENDIX A). The free variables in the model are the simulated neural control signals for hand position and the cocontraction command. In each case, the start time, the rate of equilibrium shift, and the duration of equilibrium shift may vary.

Several additional points about the organization of control should be raised. In the  model, the neurally specified equilibrium (as determined by the set of muscle s) is similar to the mechanical equilibrium of the physical system in the absence of external loads. However, with the same neurally specified equilibrium and different loads, the physical equilibrium will vary (Feldman 1986). To obtain a correspondence between physical and neurally specified equilibria, we assume that the system can set individual muscle s to achieve a desired physical equilibrium and cocontraction level in statics. This implies that the system takes account of muscle geometry and mechanics (in statics) for planning movements (Gribble et al. 1997). Although this requires knowledge of muscle forces and geometry to set neural commands, it should be emphasized that the present approach is different from proposals in which commands are established by solving the inverse dynamics of the systema procedure that, in the present model, would require neural representations of the equations of motion of the limb, reflex delays, the dependence of force on velocity, and the gradual development of muscle force over time.

	SIMULATIONS

As an argument against the equilibrium-point hypothesis, it has been claimed that complex equilibrium shifts are necessary to account for empirically observed patterns of multijoint movement. Nonmonotonic patterns of joint stiffness have been measured during multijoint movement and have been offered in support of these claims (Gomi and Kawato 1996). In this section, we use the arm model to address this issue. We assess the extent to which empirically observed time-varying patterns of stiffness associated with multijoint movements can be predicted using constant-rate equilibrium shifts and constant cocontraction commands. We also assess a related claimthat simple control signals require high levels of stiffness to generate suitable torques for rapid movements. We use the model to demonstrate that using simple equilibrium shifts, simulated rapid single-joint movements have stiffness levels that are comparable with those measured empirically (Bennett 1993). The sensitivity of the findings to changes in model parameters is reported in APPENDIX C.

Multijoint movements

We simulated the procedure used by Gomi and Kawato (1996) to estimate joint stiffness and viscosity matrices during the course of a multijoint movement. A movement similar to the one performed by subjects in the Gomi and Kawato (1996) study was tested. The hand moved from an initial position 50 cm forward from the shoulder and 20 cm to the left of the shoulder, to a final position 50 cm forward and 20 cm to the right of the shoulder, in the horizontal plane. The total movement duration was 1 s. The control signal used to generate the movement shifted in a straight line in hand space, at a constant rate, from the initial position to the final position. The duration of the equilibrium shift was 0.7 s. The cocontraction command was 5 N until movement onset, constant throughout the movement and 5 N after movement end. The magnitude of the cocontraction command during movement was set such that the same  shifts applied in statics result in an average peak muscle force of 30 N.

In each trial, at one of nine points in time during the movement (corresponding to the times used by Gomi and Kawato 1996), small force perturbations were introduced at the hand, in one of eight different directions. The magnitude of the perturbations was 0.1 N and resulted in 5- to 7-mm displacements of the hand. A total of 72 perturbed movements were simulated (9 time points × 8 directions). Using the simulated perturbation trials and the regression technique described by Gomi and Kawato (1996) (see APPENDIX B), we calculated joint stiffness and viscosity matrices for each of the nine points in time at which perturbations were applied.

Hand stiffness matrices were computed from the estimated joint stiffness matrices R using the Jacobian transformation (see Gomi and Kawato 1995 for details), and hand-stiffness ellipses were used to visualize limb stiffness at the hand. Figure 2, top, shows hand-stiffness ellipses estimated during the simulated movement. The size and orientation of the ellipses are comparable with those reported by Gomi and Kawato (1996), and likewise are larger than the corresponding ellipses during statics (see Fig. 9).

View larger version (26K): [in this window] [in a new window]   FIG. 2. Simulated hand-stiffness ellipses and joint-stiffness matrices for the arm model during multijoint movement. Constant-rate equilibrium shifts and constant cocontraction commands were used to produce the simulated movements.

Figure 2, bottom, shows the elements of the estimated joint-stiffness matrices for the arm model during movement. The terms of the joint-stiffness matrix, R, relate joint torques at the shoulder due to shoulder motion (R_ss), torques at the shoulder due to elbow motion (R_es), and so on. The basic form of the matrices is similar to those reported by Gomi and Kawato (1996), even though the equilibrium trajectory we used to generate the simulated movement was simple in shape. At the beginning of movement onset the shoulder term, R_ss, increases sharply from ~18 to ~40 N·m/rad, then decreases in the middle of movement to ~20 N·m/rad, increases again around movement end to 40 N·m/rad, and finally decreases after the end of movement to ~15 N·m/rad. The other three terms in the stiffness matrix follow roughly the same form but show a less pronounced decrease in the middle of the movement. The elbow term, R_ee increases from ~5 N·m/rad at movement start to 20-25 N·m/rad during movement, and the two double-joint terms, R_se and R_es, increase from ~2 N·m/rad at movement start to ~7-10 N·m/rad during movement. R_ee, R_es, and R_se all decrease at movement end to levels comparable with those at the onset of movement.

Using the empirically derived joint-stiffness and viscosity matrices, Gomi and Kawato (1996) compute a hypothetical equilibrium trajectory (see APPENDIX B). Their calculations are based on the assumption that joint torques can be represented with the following linear equation

(7)

where R and D are stiffness and viscosity matrices derived from the perturbation procedure, _in are the calculated joint torques (see APPENDIX B), q_eq is the equilibrium trajectory, and q and are the unperturbed movement position and velocity, respectively.

To show that the Gomi and Kawato (1996) results can be predicted using simple control signals, we used their procedure to compute a hypothetical equilibrium trajectory using the stiffness and viscosity estimates from our simulations. The trajectory that results from this calculation is shown in Fig. 3. The top panel shows the equilibrium trajectory used to generate the movement based on the  model (···), the simulated movement trajectory (- - -), and the hypothetical equilibrium trajectory derived using Gomi and Kawato's equations (), plotted in hand space. Figure 3, middle, shows the horizontal components of these trajectories plotted against time, and Fig. 3, bottom, shows the tangential velocities of the hand trajectories plotted against time.

View larger version (13K): [in this window] [in a new window]   FIG. 3. Simulated multijoint movement (- - -) generated using a constant-rate equilibrium shift (···) and a constant cocontraction command. Using the stiffness patterns shown in Fig. 2 and the algorithm proposed by Gomi and Kawato (1996), a nonmonotonic trajectory is obtained ().

The hypothetical equilibrium trajectory computed using Gomi and Kawato's procedure is "complex" in shape and does not resemble the simulated movement, which is smooth, relatively straight and looks like the movements made by subjects in the Gomi and Kawato (1996) study. Nor does it resemble the equilibrium trajectory that was used to generate the movementthe equilibrium trajectory used in the simulations is a simple constant-rate monotonic shift from one position to another. Gomi and Kawato's hypothetical equilibrium trajectory first leads then lags the simulated movement. The tangential velocity of the hypothetical equilibrium trajectory has multiple peaks and does not resemble the velocity profile of the simulated movement, which is smooth and bell-shaped. We suggest that the discrepancy between the equilibrium trajectory based on the model and the trajectory computed using Gomi and Kawato's equations arises from their use of a simplified model of force-generation (see DISCUSSION).

A number of additional points should be noted. Direct estimates of joint viscosity are not provided by Gomi and Kawato (1996). However, the present estimates correspond to values reported elsewhere. Specifically, the simulated estimates of joint viscosity have maximum values of ~2.5-3.0 Nms/rad, which is in the range of 5-7% of corresponding maximum joint stiffness. This is comparable with the relation between joint viscosity and stiffness during cyclical one-joint movements (Bennett et al. 1992) and with values for multijoint stiffness and viscosity in statics (Gomi and Osu 1996; Tsuji et al. 1995). It also should be noted that the simulations reported above have been based on constant-rate shifts in the hand equilibrium position. We also have carried out these simulations using constant-rate shifts in  space. The time-varying form and the magnitudes of joint-stiffness and viscosity matrices obtained under these conditions were comparable with those reported in the preceding text.

In summary, the kinematics and time-varying patterns of stiffness for the multijoint limb movements reported by Gomi and Kawato (1996) can be predicted using simple, constant-rate equilibrium shifts. Hypothetical equilibrium trajectories computed using Gomi and Kawato's force-generating mechanism are complex in shape even though the equilibrium trajectories that were used to generate the simulated movements are simple in form.

Single-joint movements

It has been claimed, in the context of the equilibrium point hypothesis, that high levels of stiffness are required to generate rapid movements using simple equilibrium shifts (Flash 1987; Gomi and Kawato 1996). Empirical estimates of stiffness during rapid multijoint movements are unavailable, however, stiffness estimates during rapid single-joint movements have been reported and may be compared with simulations using the present model. In this section, we show that using constant-rate equilibrium shifts, rapid single-joint movements can be generated that have stiffness levels comparable with those reported empirically (Bennett 1993). We use the arm model to simulate Bennett's procedure for estimating stiffness, and we explore the extent to which simulated stiffness during single-joint movements of various speeds matches empirical values.

Movements comparable with those performed by subjects in the Bennett (1993) study were simulated. The model was constrained to produce only single-joint elbow movements by fixing the orientation of the shoulder. Constant-rate equilibrium shifts and constant cocontraction commands were used to generate 1 rad elbow flexion and extension movements at various speeds. The duration of the 1 rad equilibrium shift was varied to produce the different movement rates (see Fig. 5, bottom right). The peak velocities of the simulated movements were comparable with those reported by Bennett (1993). Movements were generated both in the absence (unperturbed) and presence (perturbed) of external torque. In the perturbed movements, a positional perturbation (8-10°) was introduced shortly after movement start and remained on until just before movement end. The perturbation torque was controlled using a simulated servo that used proportional, integrative, and derivative control to maintain the magnitude of the positional perturbation at a constant level during the movement and to yield a velocity profile for the perturbed movements that matched the velocity profile of the unperturbed movements. Bennett's rationale for delivering perturbations in this manner was that by keeping the velocity profile of the unperturbed and perturbed movements virtually identical, any change in joint torque could be associated with the positional perturbation, and velocity-dependent changes in torque could be eliminated. Stiffness was estimated as the difference between the average torques during the unperturbed and perturbed movements, divided by the average magnitude of the positional perturbation.

View larger version (14K): [in this window] [in a new window]   FIG. 5. Simulated joint stiffness during single-joint movements. Empirical values (Bennett 1993) appear as a dashed line. Simulated movements were generated using constant-rate equilibrium shifts and constant cocontraction commands. When the cocontraction command is 3 N, predicted joint stiffness does not vary with movement velocity (, marked by cocontraction level 3 N). When the cocontraction command increases in proportion to the rate of equilibrium shift, predicted joint stiffness(, with cocontraction values indicated below) matches empirically observed values. Right: to reproduce empirically observed patterns, both the cocontraction command and the rate of equilibrium shift increase with movement speed.

Figure 4 shows an example of the simulation results for one movement speed (peak velocity of this movement is 3.00 rad/s). Figure 4, top, shows elbow angle plotted against time for an unperturbed (top trace) and perturbed (bottom trace) movement. The two vertical dashed lines indicate the boundaries of the segment of data used to estimate stiffness. The boundaries were chosen to contain the portion of the trial for which the positional perturbation was constant, and the velocity profiles of the unperturbed and perturbed movements were matched. A comparable procedure for choosing the boundaries was used in Bennett (1993). Figure 4, middle, shows the difference between the perturbed and unperturbed elbow angles plotted against time, and Fig. 4, bottom, shows the velocity of the unperturbed and perturbed movements plotted against time.

View larger version (15K): [in this window] [in a new window]   FIG. 4. Simulation of the procedure used by Bennett (1993) to estimate joint stiffness during single-joint movement. A servo-controlled positional perturbation was applied throughout single-joint movements of different speeds. The servo matched velocities during perturbed movements to those during unperturbed movements.

In an initial set of simulations, the magnitude of the cocontraction command was low (3 N) for all five movement speeds. This was done to investigate the extent to which measured stiffness changes may be observed as the velocity of the equilibrium shift increases in the absence of changes in the magnitude of the cocontraction command. (A value of 3 N allowed us to match predicted stiffness at the slowest movement speed.) The resulting stiffness estimates are shown in Fig. 5 (labeled 3). Stiffness is plotted against peak movement velocity, and the magnitude of the cocontraction command used for each simulation (in this case, 3 N), is shown. With the cocontraction command held constant at 3 N for all five movement speeds, estimated stiffness in the simulations does not increase as it does in the empirical data reported by Bennett (1993), which is shown in Fig. 5 as a dashed line. Estimated stiffness remains relatively constant at 2-4 N·m/rad for all movement speeds, suggesting that in the context of the present model, the changes in stiffness observed by Bennett (1993) that accompany increases in movement speed may be associated with an increase in the magnitude of the cocontraction command.

We tested this possibility in a second set of simulations, in which the cocontraction command was scaled in proportion to the rate of equilibrium shift. The magnitude of the cocontraction command was raised to levels for which stiffness during simulated movements matched the data reported by Bennett (1993). Note that in all simulations, the cocontraction command was held at a constant level throughout the duration of the movementonly the overall magnitude of the command was manipulated. Also note that the durations of equilibrium shifts used to generate movements with different levels of cocontraction were the same as those used for movements with low cocontraction. The resulting stiffness estimates appear in Fig. 5. By increasing the magnitude of the cocontraction command in proportion to the rate of equilibrium shift, the simulated stiffness of the arm model matches the Bennett (1993) data. In Fig. 5, top right, the magnitude of the cocontraction command is plotted against peak movement velocity. In Fig. 5, bottom right, the relationship between the rate of equilibrium shift and peak movement speed is shown. In both cases, the relationships are close to linear. Note, as well, that when cocontraction remains constant at 3 N, the simulated movement velocities are less than the velocities of movements produced by the same equilibrium shifts and higher values of cocontraction.

In summary, single-joint movement simulations can be generated at a variety of speeds using constant-rate equilibrium shifts and constant cocontraction commands. Simulated joint stiffness matches empirical values reported by Bennett (1993) when the magnitude of the cocontraction command is increased in proportion to the rate of equilibrium shift.

A further claim about the equilibrium point hypothesis is that nonmonotonic equilibrium shifts are necessary to generate torques large enough to account for rapid movements. Latash and Gottlieb (1991) report an empirical study in which perturbations were used to derive the form of hypothetical equilibrium trajectories underlying rapid single-joint movements. Using the arm model, we have simulated Latash and Gottlieb's experiment and the mathematical procedures used to infer the hypothetical equilibrium shifts. We show that even when constant rate equilibrium shifts are used to generate the simulated movements, the hypothetical equilibrium trajectories computed using Latash and Gottlieb's procedure are complex in shape. As in Gomi and Kawato's procedure, we suggest that this discrepancy is due the use of a simplified account of force generation.

We used a constant-rate equilibrium shift and a constant cocontraction command to simulate 50° elbow flexion movements. The peak movement velocity was 5 rad/s. At the onset of movement, a ramp increase in external torque was applied. The torque ramp continued throughout the movement and reached a maximum at movement end. Eight movements were simulated, each with a different magnitude of final external torque: 10, 6, 4, 1, 1, 4, 6, and 10 N·m. One trial also was simulated in the absence of external torque.

Every 5 ms during the simulated movement, values of elbow torque and joint angle for each of the eight perturbed and one nonperturbed trial were computed. These values were fit to the following linear equation, which Latash and Gottlieb (1991) use to describe the static dependence of joint torque on joint angle

(8)

T is a vector of nine joint torques,  is a vector of nine joint angles, and k₁ and k₂ are the intercept and slope, respectively, of the torque-angle relationship. By repeating this calculation every 5 ms over the course of the simulated movements, time-varying estimates of k₁ and k₂ are obtained. In Latash and Gottlieb's formulation, because torque is assumed to be a linear function of joint angle, the intercept of the torque-angle relation (the joint angle at which torque is 0), provides a measure of the hypothetical equilibrium joint angle. The ratio k₁/k₂ computed over time thus provides a measure of the hypothetical equilibrium trajectory.

Figure 6 shows the elbow angle plotted against time for a simulated unperturbed movement (), the hypothetical equilibrium trajectory computed using Latash and Gottlieb's procedure (- - -), and the equilibrium trajectory used to generate the simulated movement (···). The simulated movement is smooth with kinematics that resemble those of movements performed by subjects in the Latash and Gottlieb (1991) study. The nonmonotonic function derived from Latash and Gottlieb's equations first leads then lags the movement.

View larger version (17K): [in this window] [in a new window]   FIG. 6. Simulated kinematics () of single-joint movement, using a constant-rate equilibrium trajectory (···) and a constant cocontraction command. The algorithm proposed by Latash and Gottlieb (1991) to compute the equilibrium trajectory underlying the movement results in a nonmonotonic function (- - -).

	DISCUSSION

Abstract Introduction Discussion References

We have used a model of arm movement based on the  version of the equilibrium point hypothesis to examine a number of claims about the complexity of control signals underlying single- and multijoint movement. We have shown that multijoint movements with time-varying joint stiffness comparable with that measured empirically (Gomi and Kawato 1996) can be predicted using constant rate equilibrium shifts and constant cocontraction commands. We have shown that rapid single-joint elbow movements also can be predicted using constant-rate equilibrium shifts and constant cocontraction commands. Moreover, the magnitude of elbow joint stiffness during simulated movements of various speeds matches corresponding empirical measures (Bennett 1993).

The stiffness and viscosity patterns that were derived from our simulations were used in conjunction with the models proposed by Gomi and Kawato (1996) and Latash andGottlieb (1991) to reconstruct postulated equilibrium shifts underlying movement. The resulting trajectories were found to be nonmonotonic, as reported by Gomi and Kawato (1996) and Latash and Gottlieb (1991), even though the control signals that underlie the simulated movements were simple in form.

The hypothetical control signals derived from measures of limb impedance are dependent on the nature of the model of force generation. Neither the Gomi and Kawato (1996) nor the Latash and Gottlieb (1991) formulations include explicit muscle modelsinstead, a single "motor" generates torques at the joints. These joint torques are linearly dependent on joint stiffness, on the difference between the equilibrium and the actual trajectory, and on joint velocity and viscosity (Gomi and Kawato 1996). Hence to generate torques that first accelerate and then decelerate the limb, these hypothetical equilibrium trajectories must first lead and then lag the actual limb position.

Differences between the results of using the force generation mechanisms in the Gomi and Kawato (1996) and Latash and Gottlieb (1991) formulations and that proposed in the present paper are shown in Fig. 7. In Fig. 7A, we show the movement that results from the use of a constant rate equilibrium shift and the muscle model presented here. Figure 7B shows the same equilibrium shift used in conjunction with the force generating mechanism used by Gomi and Kawato (see Eq. 7). Figure 7C shows the postulated equilibrium trajectory required to generate the movement shown in Fig. 7A, when the Gomi and Kawato force generating mechanism is used.

View larger version (13K): [in this window] [in a new window]   FIG. 7. Effects of the shape of the equilibrium shift and the model of force generation on simulated multijoint movements. A: simulated movements that result when a constant-rate equilibrium shift is used with the model proposed in this paper. B: a constant-rate equilibrium shift was used with the force-generation model proposed by Gomi and Kawato (1996). Net deceleration torques are generated only when the limb position crosses the equilibrium (in this case, after the end of the equilibrium shift). C: to produce smooth movements comparable with that shown in A using the force-generation mechanism proposed by Gomi and Kawato, a nonmonotonic equilibrium shift is required.

In Fig. 7A, the equilibrium shift is accompanied first by torque at the shoulder, and then at about the time of peak hand velocity, the net torque changes sign and decelerates the limb. These torques arise from activity in agonist, then antagonist shoulder muscles as a result of the shifting equilibrium associated with changes to muscle threshold lengths. However, when a monotonic equilibrium shift is used with the Gomi and Kawato force generating mechanism (Fig. 7B), net deceleration torques occur only when the limb position passes the equilibrium. The resulting overshoot and oscillation arises because in this formulation the only way in which deceleration torques may be produced is when the equilibrium lags the current position of the limb. In Fig. 7C, it may be seen that to produce simulated movements comparable with those in Fig. 7A, a nonmonotonic equilibrium shift is needed that first leads and then lags the limb position. Only in this way can acceleration and deceleration torques comparable with those in Fig. 7A be produced.

A number of additional differences may be noted between the force-generating mechanism used by Gomi and Kawato (1996) and the model proposed here. In their formulation, torque varies linearly with the difference between actual and equilibrium joint angles and with joint velocity. Similarly, in Latash and Gottlieb's formulation, torque is linearly proportional to joint angle, but there is no velocity-dependent contribution to torque. A number of empirical observations suggest that in the intact preparation, the relation between muscle force and muscle length is nonlinear (Asatryan and Feldman 1965; Feldman and Orlovsky 1972). Increases in muscle activation are associated with the recruitment of progressively larger motor units, and hence larger force increments are observed (Henneman et al. 1965). In the model presented here, the gradient associated with the recruitment of motoneurons is nonlinear and is approximated using an exponential function (see Eq. 4).

Another difference concerns the dependence of force on the velocity of muscle lengthening or shortening. In Gomi and Kawato's formulation, the relation between force and velocity is assumed to be linear, and in Latash and Gottlieb's model, there is no effect of velocity on force. The nonlinear dependence of force on the velocity of muscle lengthening and shortening has been well documented (see Partridge and Benton 1981; Winters 1990; Zajac 1989 for reviews). In the muscle model used in this paper, both velocity-dependent afferent contributions to motoneurone activity and the dependence of force on the velocity of muscle shortening or lengthening are included.

A number of aspects of the  model and the organization of control signals to muscles facilitate the use of constant-rate equilibrium shifts. In Fig. 8A, a simulated movement comparable with that performed by subjects in the Gomi and Kawato (1996) study is shown. In Fig. 8, B-D, we show the effects on predicted movement kinematics of failing to model a number of phenomena, each of which is present in the  model and absent in the Gomi and Kawato (1996) and Latash and Gottlieb (1991) formulations. Each of these properties is removed from the model, and movements in their absence are examined. Figure 8B shows the effects of failing to account for the velocity-dependence of muscle threshold length (µ was set to 0). In Fig. 8C, the ability to move at different cocontraction levels is excluded (the cocontraction command was 0). In Fig. 8D, the dependence of muscle force on velocity is removed. It can be seen that these elements of the  model contribute differentially to damping during movement and to the ability of the model to produce smooth movements using simple equilibrium shifts.

View larger version (12K): [in this window] [in a new window]   FIG. 8. A: simulated movement comparable with that performed by subjects in the Gomi and Kawato (1996) study. B-D: effects of failing to model a number of physiological properties on the predicted movements. B: effect of removing the velocity dependence of muscle threshold length. In C, the ability to move at difference cocontraction levels is eliminated. In D, the dependence of active muscle force on velocity is absent. , predicted horizontal hand position; - - -, equilibrium hand position over time.

The time-varying patterns of joint stiffness and viscosity observed by Gomi and Kawato (1996) are predicted readily by the  model. Muscle stiffness is proportional to muscle force. During a movement, there are two peaks in muscle forces that are due to the activation of agonist and antagonist muscle groups and that correspond to limb acceleration and deceleration, respectively. Joint stiffness varies with joint torque and is therefore high when acceleration or deceleration is high and low during the middle of movement when acceleration decreases. The same is true for viscosity. Damping in the  model is provided by both intrinsic muscle properties (the dependence of force on velocity) and reflexes (the dependence of muscle threshold length on velocity). By differentiating Eq. 4, it can be shown that the reflex contribution to viscosity, D, varies with stiffness R (symbols as in APPENDIX B)

(9)

The intrinsic component of damping arises from the sigmoidal form of the force-velocity relationship and changes with speed. Thus at low speeds, the system has both high stiffness and viscosity, which facilitates acceleration and deceleration. In the middle of movement, both stiffness and viscosity are small, allowing the system to move at higher speeds.

In this paper, we have shown that in the context of simple point to point movements, control signals need not be complex. However, it should be emphasized that the goal of this paper has not been to demonstrate that control signals are simple rather than complex but to show that inferences about the form of the control signals are dependent on the nature of the neuromuscular plant. Indeed, more complex control signals presumably are required to produce movements with more complicated kinematics (Gribble and Ostry 1996) or to control movements in the presence of external loads (Shadmehr and Mussa-Ivaldi 1994). It should be noted that the simulations reported in the present paper have been limited to cases where direct empirical data exists against which simulation results can be compared. For this reason, movements at different speeds, different directions, or movements with more complicated kinematics have not been tested.

Even in the context of the present model, the system must use information about the consequences of the specification of s to achieve desired equilibrium positions. For example, to produce a movement in a desired time or at a specific velocity, the system must specify the rate of  shift and the level of cocontraction that will result, in conjunction with muscle properties and dynamics, in an appropriate movement. This means that the system has information about the relationship between control signals to individual muscles and resulting equilibrium trajectories, as well as information about limb dynamicsboth of which must be used to produce movement.

There have been a number of recent demonstrations that the nervous system has knowledge of its own dynamics and of the dynamics of external loads (Eliasson et al. 1995; Flanagan and Wing 1993, 1997). In one such demonstration, Flanagan and Wing (1993) showed that when transporting an object held in a precision grip, in both point-to-point and cyclical movements, the grip force exerted by subjects varied directly in anticipation of the load force as determined by the mass and the acceleration of the object. It also has been suggested that the nervous system can make adjustments to central commands to compensate for the presence of external loads acting on the limb. Shadmehr and Mussa-Ivaldi (1994) have shown that subjects can learn to make reaching movements in the presence of externally imposed force fields. With practice, hand trajectories in the force fields converge to a path similar to that observed with no external forces. This is not inconsistent with the present approachmodels such as the one presented here may be used to explore the nature of changes to central control signals that accompany motor learning.

In previous reports, we have demonstrated the importance of including models of muscle mechanics and limb dynamics when exploring the nature of control signals underlying movements. In Ostry et al. (1996), we used a model of human jaw and hyoid movement to demonstrate that context-sensitivity in jaw movements during speech need not be represented in control signals but may arise from dynamics and muscle properties. In Gribble and Ostry (1996), we used a previous version of the arm model to show that the power law relation between movement curvature and velocity observed during drawing movements may arise from mechanical and dynamical properties of the arm, and similarly, need not be explicitly planned in control signals.

The  model proposes that muscle activation is dependent upon the difference between the actual and threshold muscle length and the rate of muscle length change. Spindle primary and secondary receptors play a major role in determining this activation but presumably other sensory afferents contribute as well (for example, Ostry et al. 1997). In the model, activation is represented in terms of the combined effects of position- and velocity-dependent inputs on motoneuron recruitment rather than muscle spindle firing rates per se.

The relative gains of the velocity and position dependent inputs to motoneurons are reflected in the model in the parameter µ. µ is the ratio of velocity dependent to position dependent feedback gains in motoneuron recruitment. One possibility for relating µ to physiological parameters might be to examine the ratio of spindle primary to secondary discharge rates with changes in muscle length and velocity (for example, Houk et al. 1981; Matthews 1981). However, as noted above, in the present formulation of the  model, the effects of afferents are given in terms of their role in motoneuron recruitment and muscle activation and not in terms of spindle firing rates. Thus a straightforward comparison with these data is not possible without consideration of factors such as synaptic density and synaptic efficiency, which presumably affect the relationship between spindle firing rates and motoneuron activation.

It may be noted that whereas the model assumes that the effects of position and velocity dependent inputs to motoneurons summate, Houk et al. (1981) propose that the effects of changes in muscle length and velocity on spindle firing rates are multiplicative. In their formulation, primary and secondary discharge rates are dependent on the product of muscle length and velocity to the exponent 0.3. Although it may be difficult to reconcile these results with the additivity proposed in the  model, the present Houk et al. (1981) formulation needs to be extended to account for spindle discharge during muscle shortening and for firing rates in statics before an attempt is made to relate information about spindle afferent discharge rates to µ.

A number of limitations of the model presented here should be noted. As a simplification, we have assumed that passive muscle force varies linearly with the difference between current muscle length and muscle rest length. However, passive stiffness increases with muscle length, particularly for large amplitude stretches (see Winters 1990; Zajac 1989 for reviews). Nevertheless, the results of a number of studies indicate that the passive force-length relationship is well approximated by a linear function for extensions of up to ~75% of a muscle's maximum departure from resting length (Feldman and Orlovsky 1972; Matthews 1959; Nichols 1973). All muscle length changes in the simulations presented here fall within this range.

Note also that no attempt has been made in the present model to incorporate the contribution to force of tendon compliance. However, at least in the case of upper arm muscles, the effect of the tendon on the composite force-length dependence of the muscle plus tendon is small (Zajac 1989). In addition, although the model includes velocity- and position-dependent afferent inputs, no attempt has been made to implement force feedback due to tendon organ afferents, reciprocal or recurrent inhibition (see Bullock et al. 1996; Feldman et al. 1990). In the context of the model, reciprocal inhibition may provide additional joint stiffness and also may affect velocity-dependent reflex damping (see Feldman et al. 1990).

The present version of the arm model is a two-dimensional, two df planar model with six musclesit is possible that more detailed predictions may be obtained by extending the model to three dimensions, by modeling more mechanical df (supination/pronation at the elbow or adduction/abduction and rotation at the shoulder, for example) or by including more muscles. An additional limitation concerns the balance of forces produced by the cocontraction command. At present we know of no empirical reports documenting the relative balance of muscle cocontractive forces during movement. In the present version of the model, we have assumed that the cocontraction command increases all muscle forces in the most equal proportions, while maintaining the constraint that there is no motion. Other cocontraction commands associated with different distributions of muscle force also may be used.

	ACKNOWLEDGEMENTS

  The authors thank A. Feldman and E. Saltzman for comments.

  This research was supported by National Institute of Deafness and Other Communication Disorders Grant DC-00594, by National Sciences and Engineering Research Council (Canada), and Fonds pour la Formation de Chercheurs et l'Aide à la Recherche (Québec).

	APPENDIX A: EMPIRICALLY BASED ESTIMATES FOR  AND µ

In previous versions of the model, values for  have been set in proportion to the physiological cross-sectional area (PCSA) of muscles and thus provide relative estimates of muscle force-generating ability. Values of µ have been set so that damping in simulated movements is comparable with that observed empirically. In the present paper, we have established parameter estimates for  and µ on the basis of both PCSA values and empirical measures of joint stiffness and viscosity in statics. In statics, predicted stiffness and viscosity vary directly with values of  and µ, respectively, and thus it is possible to establish parameter estimates for and µ by matching simulated stiffness and viscosity values to those observed empirically.

Empirical estimates of stiffness and viscosity in statics have been reported by a number of researchers (Gomi and Kawato 1995; Mussa-Ivaldi et al. 1985; Tsuji et al. 1995). In these experiments, subjects grasped a 2 df manipulandum while maintaining a limb posture in the horizontal plane. Perturbations displaced the hand from the rest position in different directions, and the resulting restoring forces were measured. Gomi and Kawato (1995) and Tsuji et al. (1995) fit the time-varying response of the limb to a second-order equation containing constant inertia, stiffness, and viscosity terms. Mussa-Ivaldi et al. (1985) estimated stiffness only by measuring restoring forces associated with positional displacements. The results of these experiments provide estimates of joint stiffness and viscosity of the following form

(A1)

(A2)

where shoulder and elbow torques T_s and T_e opposing the perturbation are related to joint displacements _s and _e by the stiffness matrix R and to joint velocities _s and _e by the viscosity matrix D. The terms in the stiffness matrix R relate torques at the shoulder due to shoulder displacements (R_ss), torques at the shoulder due to elbow displacements (R_es), and so forth. The four corresponding terms in the viscosity matrix D relate restoring torques at each joint to corresponding changes in joint velocity.

View larger version (18K): [in this window] [in a new window]   FIG. A1. Simulated hand-stiffness ellipses and joint-stiffness and viscosity matrices for the arm model for 5 static postures.

Using these techniques, Gomi and Kawato (1995) and Tsuji et al. (1995) report estimates of static joint stiffness that range from 5 to 20 N·m/rad for R_ss, 4 to 13 N·m/rad for R_ee, and 1 to 7 N·m/rad for R_es and R_se, which tended to be about equal in value. Values reported by Mussa-Ivaldi et al. (1985) are somewhat higher, ranging from 16 to 40 N·m/rad for R_ss, 13 to 41 N·m/rad for R_ee, and 2 to 25 N·m/rad for R_es and R_se. Joint viscosity matrices reported by Tsuji et al. (1995) had values in the range of 7-12% of the corresponding values of joint stiffness. Values for different subjects ranged from 0.4 to 1.0 Nms/rad for D_ss, 0.2 to 0.6 Nms/rad for D_ee, and 0.1 to 0.5 Nms/rad for D_es and D_se.

In the context of the model, we used a procedure comparable with that used by Tsuji et al. (1995) to estimate stiffness and viscosity matrices in statics. The aim is to obtain empirically based parameter estimates for  and µ. Values for  were obtained by scaling estimates of physiological cross-sectional area by a single constant value such that the predicted joint stiffness of the model matches empirical estimates. Similarly, the magnitude of µ was established such that the joint viscosity of the model corresponds to empirically measured viscosity (Tsuji et al. 1995). Note that this procedure was carried out after the values of all other muscle model parameters were fixed.

To estimate joint stiffness and viscosity, a number of different postures corresponding to those used by Gomi and Kawato (1995) were tested (see Fig. A1, top). The limb was in equilibrium at each posture (velocity was 0), the cocontraction command was 5 N, and the simulated limb was displaced by 5 mm in eight different directions in hand space. The time-varying shoulder and elbow joint torques T_s(t) and T_e(t) opposing the perturbation were computed and related to the joint displacements _s(t) and _e(t) and joint velocities _s(t) and _e(t) to compute the joint stiffness matrix R and the joint viscosity matrix D (see Eqs. A1 and A2).

PCSA estimates for each muscle were scaled by 1 N/cm² to yield joint-stiffness matrices comparable with those reported by Gomi and Kawato (1995) and Tsuji et al. (1995). Over five different workspace positions, values ranged from 11 to 13 Nms/rad for R_ss, 6 to 7 Nms/rad for R_ee, and 2 to 3 Nms/rad for R_es and R_se (see Fig. A1, middle).

With µ set to 0.06 s, predicted joint viscosities in the model correspond to measured joint viscosities (Tsuji et al. 1995). Values of the viscosity matrices over the five postures tested ranged from 0.7 to 0.9 Nms/rad for D_ss, 0.3 to 0.4 Nms/rad for D_ee, and 0.1 to 0.3 Nms/rad for D_es and D_se (see Fig. A1, bottom).

It should be emphasized that the procedures described above establish values for µ and , which match empirical and model data in statics only, and in no way constrain the time-varying form of stiffness and viscosity matrices during movement.

	APPENDIX B: PROCEDURE FOR ESTIMATING LIMB IMPEDANCE AND HYPOTHETICAL EQUILIBRIUM TRAJECTORIES DURING MOVEMENT

The dynamics of a two-joint limb are given by

(B1)

where q is position, I(q) is the matrix of limb inertia, H(q, ) is the coriolis-centrifugal force vector, _m (q, , u) is torque due to muscle activation, _ext is external torque, and u is the descending motor command (see Gomi and Kawato 1996).

To estimate stiffness (R), viscosity (D), and inertia (I) matrices from the perturbation data, Eq. B1 is differentiated to give the following variational equation (Gomi and Kawato 1995, 1996)

(B2)

where q, , and represent time-varying differences in position, velocity, and acceleration imposed by the external torque perturbations, and _ext denotes the time-varying external torque imposed by the manipulandum. The coriolis-centrifugal force vector H and the inertia matrix I are calculated from limb geometry.

The form of Eq. B2 is straightforwardthe left-hand side represents torques due to the dynamics of the limb (inertial, coriolis, and centrifugal forces), and the right-hand side represents change in joint torques due to muscles (D  Rq), and external torques _ext. Note that when this variational equation is used, torques are linearly dependent on position and velocity.

Although Gomi and Kawato (1996) estimated the terms representing limb dynamics (Eq. B2, left) from the perturbation data, we calculated them directly using model parameters. Our values obtained with the model by direct calculation are comparable with those used by Gomi and Kawato (1995).

Following Gomi and Kawato (1996), Eq. B2 was linearized to solve for stiffness and viscosity matrices R and D by subtracting the left-hand side of Eq. B2, (represented in following sections by vector Z), from the external torques _ext

(B3)

To obtain numerical estimates for the terms in the matrices, data from all eight perturbation directions were used at each of the nine time points at which the values in the matrices were to be estimated. As in Gomi and Kawato (1996), we used 0.28 s of data after the onset of the perturbations for the regressions.

View larger version (11K): [in this window] [in a new window]   FIG. C1. Range of parameter variation in sensitivity analyses for c (force-length relation),  (graded force development), and f4 (force velocity). Vertical axes indicate force (N) for the force-length relation and normalized force for graded force development and the force-velocity relation. Minimum (- - -), center (), and maximum (···) levels are indicated for each parameter. Values of c from left to right are 10, 112, and 220 mm1. Values of  from left to right are 1, 15, and 30 ms. Values of f4 from shallow to steep are 10, 58, and 120 s/m. (c = 112 mm1,  = 15 ms, and f4 = 58 s/m are used for the simulations reported in the main body of the paper.) For force velocity, negative velocities correspond to muscle lengthening and positive values represent muscle shortening.

Gomi and Kawato (1996) use their empirically derived joint stiffness and viscosity measures to compute estimates of the equilibrium trajectories underlying the movements performed by subjects. To do this, they assume that joint torques can be represented with a linear equation

(B4)

where _in is joint torque, R and D are stiffness and viscosity, respectively, q_eq is the postulated equilibrium trajectory, and q and are the unperturbed movement position and velocity, respectively. Given estimates of stiffness and viscosity, Gomi and Kawato (1996) decompose torque _in into inertial, coriolis, and centrifugal components, and solve Eq. 7 for the postulated equilibrium trajectory q_eq

(B5)

	APPENDIX C: ANALYSIS OF THE SENSITIVITY OF PREDICTED LIMB IMPEDANCE TO MODEL PARAMETER VARIATION

Simulations were carried out to assess the sensitivity of the results to changes in model parameters. We varied the values of all muscle model parameters and examined the effect of these variations on simulated joint stiffness and viscosity in statics and during movement. The effects of the following parameters were assessed: , the scale parameter for muscle force-generating ability; the magnitude of the cocontraction command; µ, the dependence of muscle threshold length on velocity; d, reflex delay; , the time constant for graded force development; c, the form parameter for the force-length relationship; parameter f₄, which determines the form of the force-velocity relationship; and k, passive muscle stiffness.

Sensitivity in statics

The sensitivity of stiffness and viscosity to parameter change was determined quantitatively by varying each parameter, one at a time, and computing the associated change in simulated joint stiffness and viscosity. For each simulation, the Tsuji et al. (1995) procedure described in APPENDIX A was used. The parameter values used in the main body of the paper were selected as center values for the sensitivity analyses. Each parameter was varied in 20 equal steps over a range of approximately ±100% of the center values. (However, note that nonzero values generally were used for the lower limits.) When varying values of a given parameter, the values of the other muscle model parameters were set to their center values (see Tables C1 and C2 for center values, in parentheses, and ranges for each parameter). As an aid to visualization of the parameter ranges that were tested, Fig. C1 shows the effect of extreme parameter values on the force-length and force-velocity relationships as well as on the time course of graded force development.

View larger version (19K): [in this window] [in a new window]   FIG. C2. Estimates of predicted stiffness (Rss) and viscosity (Dss) resulting from model parameter variation. The figure shows predicted values of stiffness and viscosity resulting when each of the muscle model parameters are varied one at a time. Ranges of parameter variation are indicated. , force-length scale parameter; coc, cocontraction level; delay, reflex delay; gfd, graded force development; c, force-length form parameter; fv, force-velocity; k, passive stiffness; dashed lines indicate ranges of empirical estimates of stiffness and viscosity reported by Tsuji et al. (1995). Note that , cocontraction, and c have the greatest effect on simulated stiffness, while all parameters but , c, and k affect simulated viscosity. Numerical estimates of the slopes of these relations and those for the other 3 terms in the stiffness and viscosity matrices are given in Tables C1 and C2.

View larger version (19K): [in this window] [in a new window]   FIG. C3. Estimates of simulated stiffness (Rss) and viscosity (Dss) during movement resulting from changes in individual model parameters. Means (bold) ±1 SD are shown. Abbreviations as in Fig. C2. Numerical estimates of the sensitivity of all 4 terms in the stiffness and viscosity matrices to model parameter change are given in Fig. C4.

The estimates of stiffness and viscosity resulting from parameter variation are given in Fig. C2 for the R_ss and D_ss terms of the stiffness and viscosity matrices. For each parameter, the abscissa corresponds to the parameter range that is tested, and the ordinate gives the predicted stiffness or viscosity. In Fig. C2, dashed lines give the range of empirical estimates reported by Tsuji et al. (1995). Simulated stiffness and viscosity are seen to vary linearly with changes in the value of each parameter (except for stiffness at low levels of the cocontraction command). Linear relationships were obtained for the other terms of the stiffness and viscosity matrices as well. The figure also shows that , cocontraction, and c are the primary determinants of stiffness, whereas all variables except , c, and k affect simulated viscosity.

Tables C1 and C2 give numerical estimates of the sensitivity of each of the terms in the stiffness and viscosity matrices to changes in the muscle model parameters. The percent change in stiffness and viscosity resulting from a 1% change in each parameter is shown. For example, a 1% increase in  results in a 0.75% increase in stiffness for R_ss and a decrease in viscosity of 0.05% for D_ss. As indicated in Fig. C2, the main parameters affecting stiffness are c, , and cocontraction, whereas all variables except for , c, and k affect simulated viscosity.

View larger version (21K): [in this window] [in a new window]   FIG. C4. Sensitivity of predicted stiffness and viscosity during movement to changes in individual model parameters. Horizontal axis shows time, and the vertical axes show the change in stiffness (or viscosity) resulting from a 1% change in each model parameter. , Rss and Dss; - · -, Ree and Dee; - - - and ···, es and se terms, respectively. Abbreviations as in Fig. C2.

Sensitivity during movement

A procedure comparable with that described above for statics was used to assess the sensitivity, during movement, of simulated stiffness and viscosity to parameter change. For each set of parameters, the Gomi and Kawato (1996) procedure was used to estimate stiffness and viscosity during movement.

As in the statics analyses, each of the parameters was varied one at a time, and the effects on stiffness and viscosity during movement were examined (note that during a given simulated movement, all parameters were held constant throughout the movement.) The parameters were varied in 20 equal steps over the same ranges reported in Tables C1 and C2 (except for the magnitude of the cocontraction command, which varied between 10 and 60 N). When varying a given parameter, the values of the other parameters were set to the middle values in their range (the center value of the cocontraction command was 30 N during movement.)

The resulting stiffness and viscosity estimates over the course of simulated movements are given in Fig. C3 for R_ss and D_ss, respectively. The mean SD is shown. It can be seen that the basic double-peaked form of the stiffness and viscosity functions remains unchanged when parameter values are varied. Changes to , cocontraction, and c have the greatest effects on stiffness and tend to raise or lower the overall level of the function. In contrast, the magnitude of the graded force development and force-velocity parameters differentially affect simulated stiffness over timethe effects of graded force development are most prominent at the peaks of the function, whereas the magnitude of the force-velocity parameter affects simulated stiffness during the middle of movement. Lower values of f₄ (which are associated with more shallow force-velocity gradients) result in higher levels of stiffness in the middle. For the same range of parameters, simulated viscosity during movement is affected most by the cocontraction command, graded force development, c, and force velocity. The effects of cocontraction, graded force development, and force velocity are greatest at the peaks of the function, whereas the effects of c are greatest in the middle.

At each of the nine time points during the movement at which stiffness and viscosity were estimated, numerical estimates were obtained of the sensitivity of stiffness and viscosity to parameter change. For each of the four terms in the stiffness and viscosity matrices, the change in stiffness (or viscosity) resulting from a 1% change in each model parameter is shown for the nine time points during movement at which simulated limb impedance was computed. Figure C4 indicates that the sensitivity of simulated stiffness and viscosity varies during movement and tends to be greatest at the beginning and end when muscles are active. The sensitivity of stiffness to changes in the force-velocity parameter, and the sensitivity of viscosity to changes in c show different patternsin each case sensitivity is greatest during the middle of movement.

	FOOTNOTES

  Address for reprint requests: D. J. Ostry, Dept. of Psychology, McGill University, 1205 Dr. Penfield Ave., Montreal, Quebec H3A 1B1, Canada.

  Received 24 April 1997; accepted in final form 14 November 1997.

	REFERENCES

Abstract Introduction Discussion References

• AN, K., HUI, F., MORREY, B., LINSCHEID, R., CHAO, E. Muscles across the elbow joint: a biomechanical analysis. J. Biomech. 14: 659-669, 1981.[Medline]
• AN, K., KAUFMAN, K., CHAO, E. Physiological considerations of muscle force through the elbow joint. J. Biomech. 22: 1249-1256, 1989.[Medline]
• ASATRYAN, D. G., FELDMAN, A. G. Functional tuning of the nervous system with control of movements or maintenance of a steady posture. I. Mechanographic analysis of the work of the limb on execution of a postural task. Biophysics 10: 925-935, 1965.
• BENNETT, D. J. Torques generated at the human elbow joint in response to constant position errors imposed during voluntary movements. Exp. Brain Res. 95: 488-498, 1993.[Medline]
• BENNETT, D. J., HOLLERBACH, J., XU, Y., HUNTER, I. Time-varying stiffness of human elbow joint during cyclic voluntary movement. Exp. Brain Res. 88: 433-442, 1992.[Medline]
• BIZZI, E., ACCORNERO, N., CHAPPLE, W., HOGAN, N. Posture control and trajectory formation during arm movement. J. Neurosci. 4: 2738-2744, 1984.[Abstract]
• BULLOCK, D., GROSSBERG, S., AND GUENTHER, F. Neural network modeling of sensory-motor control in animals. In: Advances in Motor Learning, edited by H. Zelaznik. Champaign, IL: Human Kinetics Publishers, 1996, p. 261-292.
• ELIASSON, A., FORSSBERG, H., IBUTA, K., APEL, I., WESTLING, G., JOHANSSON, R. Development of human precision grip. V. anticipatory and triggered grip actions during sudden loading. Exp. Brain Res. 106: 425-433, 1995.[Medline]
• FELDMAN, A. G. Functional tuning of the nervous system with control of movement or maintenance of a steady posture. II. Controllable parameters of the muscle. Biophysics 11: 565-578, 1966.
• FELDMAN, A. G. Once more on the equilibrium-point hypothesis ( model) for motor control. J. Mot. Behav. 18: 17-54, 1986. [Medline]
• FELDMAN, A. G., ADAMOVICH, S. V., LEVIN, M. F. The relationship between control, kinematic and electromyographic variables in fast single-joint movements in humans. Exp. Brain Res. 103: 440-450, 1995.[Medline]
• FELDMAN, A. G., ADAMOVICH, S. V., OSTRY, D. J., FLANAGAN, J. R. The origin of electromyogramsexplanations based on the equilibrium point hypothesis. In: Multiple Muscle Systems: Biomechanics and Movement Organization, edited by J. Winters, and S. Woo . New York: Springer-Verlag, 1990, p. 195-213
• FELDMAN, A. G., ORLOVSKY, G. N. The influence of different descending systems on the tonic reflex in the cat. Exp. Neurol. 37: 481-494, 1972.[Medline]
• FLANAGAN, J. R., OSTRY, D. J., FELDMAN, A. G. Control of human jaw and multi-joint arm movements. In: Cerebral Control of Speech and Limb Movements, edited by and G. Hammond . New York: Springer-Verlag, 1990, p. 29-58
• FLANAGAN, J. R., OSTRY, D. J., FELDMAN, A. G. Control of trajectory modifications in target-directed reaching. J. Mot. Behav. 25: 140-152, 1993. [Medline]
• FLANAGAN, J. R., WING, A. M. Modulation of grip force with load force during point-to-point arm movements. Exp. Brain Res. 95: 131-143, 1993.[Medline]
• FLANAGAN, J. R., WING, A. M. The role of internal models in motion planning and controlevidence from grip force adjustments during movements of hand-held loads. J. Neurosci. 17: 1519-1528, 1997.[Abstract/Free Full Text]
• FLASH, T. The control of hand equilibrium trajectories in multi-joint arm movements. Biol. Cybern. 57: 57-74, 1987.[Medline]
• GOMI, H. AND KAWATO, M. Task dependent stiffness of human multi-joint arm during point-to-point movement. Nippon Telegraph and Telephone ISRL-95-4, 1995.
• GOMI, H., KAWATO, M. Equilibrium-point control hypothesis examined by measured arm stiffness during multijoint movement. Science 272: 117-120, 1996.[Abstract]
• GOMI, H. AND OSU, R. Human arm stiffness and viscosity in interaction with environments on a horizontal plane. Nippon Telegraph and Telephone ISRL-96-3, 1996.
• GRIBBLE, P. L., LABOISSIÈRE, R., OSTRY, D. J. Control of human arm and jaw motion: issues related to musculo-skeletal geometry. In: Self-Organization, Computational Maps and Motor Control, edited by P. Morasso, and V. Sanguineti . North Holland: Elselvier, 1997, vol. 118, p. 483-506
• GRIBBLE, P. L., OSTRY, D. J. Origins of the power law relation between movement velocity and curvature: modeling the effects of muscle mechanics and limb dynamics. J. Neurophysiol. 76: 2853-2860, 1996.[Abstract/Free Full Text]
• HAINAUT, K., DUCHATEAU, J., DESMEDT, J. Differential effects on slow and fast motor units of different programs of brief daily muscle training in man. In: Motor Unit Types, Recruitment and Plasticity in Health and Disease, edited by and J. Desmet . Basel: Karger, 1981, vol. 9, p. 241-249
• HENNEMAN, E., SOMJEN, G., CARPENTER, D. O. Functional signifigance of cell site in spinal motoneurons. J. Neurophysiol. 28: 560-580, 1965.[Free Full Text]
• HOGAN, N. An organizing principle for a class of voluntary movements. J. Neurosci. 4: 2745-2754, 1984.[Abstract]
• HOLLERBACH, J., FLASH, T. Dynamic interactions between limb segments during planar arm movement. Biol. Cybern. 44: 67-77, 1982.[Medline]
• , sect. 1, vol. II, part 1, chapt. 8, p. 257-323. HOUK, J., RYMER, W. Neural control of muscle length and tension. In: Handbook of Physiology. The Nervous System. Motor Control., . Bethesda, MA: Am. Physiol. Soc., 1981
• HOUK, J., RYMER, W., CRAGO, P. Dependence of dynamic response of spindle receptors on muscle length and velocity. J. Neurophysiol. 46: 143-166, 1981.[Free Full Text]
• JOYCE, G. C., RACK, P.M.H. Isotonic lengthening and shortening movements of cat soleus muscle. J. Physiol. (Lond.) 204: 475-491, 1969.[Abstract/Free Full Text]
• LABOISSIÈRE, R., OSTRY, D. J., FELDMAN, A. G. Control of multi-muscle systems: human jaw and hyoid movements. Biol. Cybern. 74: 373-384, 1996.[Medline]
• LATASH, M. L., GOTTLIEB, G. L. Reconstruction of joint compliant characteristics during fast and slow movements. Neuroscience 43: 697-712, 1991.[Medline]
• LEVIN, M. F., FELDMAN, A. G. The role of stretch reflex threshold regulation in normal and impaired motor control. Brain Res. 657: 23-30, 1994.[Medline]
• MATTHEWS, P. The dependence of tension upon extension in the stretch reflex of the soleus muscle of the decerebrate cat. J. Physiol. (Lond.) 147: 521-546, 1959.
• , sect. 1, vol. II, part 1, chapt. 6, p. 189-228. MATTHEWS, P. Muscle spindles: their messages and their fusimotor supply. In: Handbook of Physiology. The Nervous System. Motor Control., . Bethesda, MA: Am. Physiol. Soc., 1981
• MUSSA-IVALDI, F. A., HOGAN, N., BIZZI, E. Neural, mechanical, and geometric factors subserving arm posture in humans. J. Neurosci. 5: 2732-2743, 1985.[Abstract]
• NICHOLS, T. Reflex and non-reflex stiffness of soleus muscle in the cat. In: Control of Posture and Locomotion, edited by R. Stein, K. Pearson, R. Smith, and J. Redford . New York: Plenum, 1973, p. 407-410
• OSTRY, D. J., GRIBBLE, P. L., GRACCO, V. L. Is context-sensitivity in speech kinematics centrally planned? J. Neurosci. 16: 1570-1579, 1996.[Abstract/Free Full Text]
• OSTRY, D. J., GRIBBLE, P. L., LEVIN, M. F., FELDMAN, A. G. Phasic and tonic stretch reflexes in muscles with few muscle spindles: human jaw opener muscles. Exp. Brain Res. 116: 299-308, 1997.[Medline]
• , sect. 1, vol. II, part 1, chapt. 3, p. 43-106. PARTRIDGE, L., BENTON, L. Muscle, the motor. In: Handbook of Physiology. The Nervous System. Motor Control., . Bethesda, MA: Am. Physiol. Soc., 1981
• SHADMEHR, R., MUSSA-IVALDI, F. Adaptive representation of dynamics during learning of a motor task. J. Neurosci. 14: 3208-3224, 1994.[Abstract]
• TSUJI, T., MORASSO, P. G., GOTO, K., ITO, K. Human hand impedance characteristics during maintained posture. Biol. Cybern. 72: 475-485, 1995.[Medline]
• WINTERS, J. Hill-based muscle models: a systems engineering perspective. In: Multiple Muscle Systems: Biomechanics and Movement Organization, edited by J. Winters, and S. Woo . New York: Springer-Verlag, 1990, p. 69-93
• WINTERS, J. AND WOO, S.-Y. (Eds.) Multiple Muscle Systems: Biomechanics and Movement Organization. New York: Springer-Verlag, 1990.
• WON, J., HOGAN, N. Stability properties of human reaching movements. Exp. Brain Res. 107: 125-136, 1995.[Medline]
• ZAJAC, F. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. CRC Crit. Rev. Biomed. Eng. 17: 359-415, 1989.

0022-3077/98 $5.00 Copyright ©1998 The American Physiological Society