Inverse dynamic optimization including muscular dynamics, a new simulation method applied to goal directed movements

Pergamon

J. Biomechonics.

Vol. 27. No. 1, pp. 953-960, Elsevier

Snence

1994 Ltd

Printedin Great Britain cm-9290/94$7.00+0.00

CO21-9290(93)EO004-9

TECHNICAL

NOTE

INVERSE DYNAMIC OPTIMIZATION INCLUDING MUSCULAR DYNAMICS, A NEW SIMULATION METHOD APPLIED TO GOAL DIRECTED MOVEMENTS R. HAPPEE* Man-Machine-Systems Group, Delft University of Technology, Faculty of Mechanical Engineering and Marine Technology, Mekelweg 2, 2628 CD Delft, The Netherlands Abstract-This paper presents a new method for estimating muscular force and activation from experimental kinematic data. The method combines conventional inverse dynamics with optimization utilizing a dynamic muscle model. The method uses only very limited computational power, which makes it a useful tool especially for complex systems like the shoulder or the locomotor system. The net torques/forces are calculated by using conventional inverse dynamics. A solution of the load sharing problem is determined by minimization of the weighted sum of squared muscle forces. The load sharing problem is solved with a dynamic constraint reflecting physiological muscle properties. This constraint takes into account the nonlinear dynamics of the contractile element (CE) and the series elastic element (SE), active state dynamics and neural excitation dynamics. This physiological constraint is determined with an inverse muscle model. With this model, muscular states and neural inputs are also estimated. The method of inverse dynamics requires position, velocity and acceleration signals as input. A method to prepare such signals from noisy measured data is presented.

INTEODUCTION Simulation studies of human movement can be of great importance in the fields of preventive health care, rehabilitation, basic motion research and sports. Simulations can provide insight into matters which cannot be measured with a reasonable effort. Among these are the activity of deeper lying muscles and the forces occurring in joints and ligaments. Simulations require a model and input signals. Models and parameters of several parts of the movement system can be found in the literature (Audu and Davy, 1985; Hatze, 1981; Winters and Stark, 1985, 1987,1988; Yamaguchi et al., 1990; Zajac, 1989). For voluntary movements it is very logical to define the efferent muscular signals (alpha activation) as inputs which yields a forward simulation (Fig. 1). However, the muscular inputs cannot be measured and must therefore be estimated somehow. This can be done with dynamic optimization (or optimal control): input patterns are computed such that a certain mathematical criterion is optimized. Such a criterion has to reflect general aspects such as the metabolic cost of movement and task specific aspects such as maximal energy or impulse or minimal movement duration. Due to the complex nonlinear properties of the movement system, dynamic optimization generally cannot be performed analytically (Seif Naraghi and Winters, 1990), so numerical methods are used. These numerical methods yield extremely

Received in final form 6 October 1993. TN0 Crash-Safety Research Centre, Biomechanics section, Schoemakerstraat 97,2628 VK Delft, The Netherlands, Tel + 31-15-697024. Fax + 31-15-624321. E-mail: [email protected].

large computational costs for the analysis of large-scale systems (Pandy et al., 1992). Alternatively, simulations can be performed with the method of inverse dynamics. In such inverse simulations the displacements and the external forces are used as inputs (Fig. 2). These displacements and forces must be recorded in relevant experiments. Limb displacements can be recorded with video systems. Displacements and forces can be recorded with mechanical devices. The method of inverse dynamics has the advantage that the solution for each position (sample) can be computed independently of the solution for other samples. This greatly reduces computational complexity. Due to these practical advantages, the method of inverse dynamics can be used to analyze complex systems like the locomotor system (Crowninshield and Brand, 1981; Patriarco et al., 1981; Pedersen et al., 1987; Pedotti et a[., 1978) and the shoulder mechanism (Karlsson and Peterson, 1992; Van der Helm, 1993a,b). Until now, muscular dynamics have not been considered in the method of inverse dynamics. It is not taken into account that muscular dynamics constrain the rate of force increase and even more strongly constrain the rate of force decrease. Thus, muscle forces may be calculated that cannot be realized with any muscular input or that can only be generated with large metabolic costs. In this paper, a method is presented that is based on inverse dynamics and that takes into account muscular dynamics at limited computational costs. This method is used to predict the coordinated control of an agonist/antagonist muscle pair in goal directed movements. METHODS AND RESULTS The method of inverse dynamics requires the trajectory for all degrees of freedom as input data including its first and second derivatives towards time. For each degree of freedom

953

954

Technical Note

Fig. 1. Forward simulation: u=neural input, I, =muscle length, F=muscle force, T=net torque, r,=extemal load, @=joint angle, x=position, L(0) describes muscle lengths as a function of joint angles, thereby its derivative L’(0) describes the moment arms.

the desired net torque (or force) can be. calculated analytically from the trajectory and the external forces (Fig. 2). As a following step it has to be determined how these desired net torques are obtained by muscular forces: the load sharing problem has to be solved. When the number of muscles exceeds the number of degrees of freedom the load sharing problem has an infinite number of solutions. So, to obtain a unique solution extra information has to be added. This is usually done by formulating a criterion that is to be optimized (minimized or maximized). This yields an inverse dynamic optimization (also called static optimization). The criterion optimized can represent the metabolic cost of muscle contraction but can also be task dependent. Several criteria have been proposed in the literature (Crowninshield and Brand, 1981; Dul et al., 1984a, b, Happee and Van der Helm, 1993a; Patriarco et al., 1981; Pedersen et al., 1987). An important class of criteria (J) is formed by summed muscle forces (F) with power p and weighting factors w: J= i

wy;.

(1)

i-1

The weighting factors w can be chosen as a function of muscle size parameters (Happee and Van der Helm, 1993a). With linear criteria (p=l), a minimal number of muscles having the largest moment arms and/or physiological crosssectional areas will be recruited (Pedersen er al., 1987). With more than proportional criteria (pz 1) ‘less favorable’ muscles will also be recruited. This is realistic as generally, in experiments, all potential synergist muscles are found to be active. Reasonably realistic results are generally obtained with n = 2 (Hannee and Van der Helm. 1993~ Pedotti et aI.. 1978jor with p-2 3 (Crowninshield and Brand, 1981; Peder: sen et al., 1987). Such a solution of the load sharing problem can be determined independently per position or per sample. The solution for one sample does not affect the solution for other samples. However, due to muscular dynamics, muscle force depends on the preceding input and movement history. The fact that neural inputs are bounded constrains the rate of

X’ X”

__D

inverse $ geometry

w inverse b inertia

T _F)

force increase and decrease. We will develop a method to include this physiological constraint in the solution of the load sharing problem. Muscle-joint model

In Happee (1992b) muscular activity during goal directed movements was studied experimentally. In Happee (1992a) dynamic optimixations were performed with a model matching these movements. In this model the shoulder has been described as a joint with only one degree of freedom: an anteflexion/retroflexion in the sagittal plane. Two muscles have been defined in the model: an anteflexor and a retroflexor. This yields the simplest model that can be used to study agonist/antagonist coordination. [For advanced shoulder models see Veeger et al. (1991), Van der Helm et al. (1992) and Happee and Van der Helm (1993b).] The anteflexor reflects the lumped effects of mainly the pectoralis major and the deltoideus pars anterior. The retroflexor reflects the lumped effects of mainly the latissimus dorsi and the deltoideus pars posterior. A lumped-parameter Hill-type. muscle model was adopted from Winters and Stark (1985, 1987,1988). The model equations are given in the left column of Table 1 and the parameters are given in Table 2. Neural excitation and active state dynamics are described by two first-order systems in series. In the literature these dynamics are often described together by one first-order system. Such a reduction of the model order considerably affects the predicted neural inputs (Happee, 1992a). The muscle model contains three states: neural excitation (E), active state (A) and the length (l,,) of the contractile element (CE). A series elastic element (SE) is included. An active force-length relation has not yet been implemented in the CE but this can be done easily by taking F,.. to be dependent on the muscle length (which is time dependent). This is, however, not important as long as only relatively small changes of muscle length are considered. A passive force-length relation of a parallel elastic element (PE) can be used to define F,,. In the current model F, is chosen to be zero. This is acceptable as no extreme joint angles are considered. Force-velocity

load Inverse + ‘.., sharing i-_-D muscular 4 j I> dynamics ,,, .,.I’ “._._ ...._.‘............_..__.........-..

D w

-kl

Fig. 2. Inverse simulation, symbols as in Fig. 1.

NEW .“.

Technical Note

955

Table 1. Model equations: Y = neural input (0 Q u 6 lk E = neural excitation; A = active state; A, = active state resealed such that it is greater than A,i,; CE = contractile element; SE = series elastic element; PE = parallel elastic element; Z, = total muscle length; I,. = length CE; 1.. = length SE; 0 = joint angle; LB0= muscle length for 0 = 0; F = active muscle force; F, = active force relative to the isometric force that would result from the current active state; F, = passive muscle force; E’ = dE/dt; i = sample number, DT = time step between two samples; parameter definitions and values are given in Table 2 Element

Continuous form

Discrete analogon

Neural excitation

u-E E’ = T ne

Ei=Ei_r+(ui_r-Ei_r)[l-exp(-F)]

Active state

if E > A then T. = T,,, else T. = Td.

ifEi-r r Ai_I then T. = TaE,else T. = Td

E-A A’ = T.

A,=A,_,+(Ei_,-A,-,)[I-exp(

Contractile element (force-velocity)

-F)]

Simulation of relations in the left column, assuming a constant active state over DT and with linear interpolation of [, over DT. A variable time step Kutta-Merson algorithm is used for this simulation.

A, = &in + (1 - &,,.)A F, = A rn.X r G.=,,[,,($-P,)-P,]

Series elastic (not dynamic)

5, = IIn- J,, F=

Parallel elastic (not dynamic)

F,, = 0

F,, = 0

Joint

T,,, = C(F + F,,)ma

Not discretired

0” =

- B&” + T,,,

1, = LB0 -2a0 Table 2. Model parameters: the two muscles differ only in the sign of their moment arms

Description

Parameter

Value

Maximum isometric force Minimal or resting active state Moment arm Limb inertia defined towards joint Joint damping Series elastic stretch with maximal isometric force Series elastic shape parameter Hill, maximal shortening velocity Hill, shape parameter Time constant, increasing active state Time constant, decreasing active state Time constant, neural excitation

Fmm Amin ma

2500 N 0.005 *0.04 m 0.25 kg mz 0.2Nmsrad-r 0.0201 m 3 1 ms-’ 0.25 10 ms 50 ms 40ms

J,

B* SE., S&I, MV,,’ MV,,* T*, Tcl. T“L

*Derived parameters: PI =0.54, P r= 1.38; Ps=tan(l/Pr relations described in the literature generally consist of separate equations for lengthening and shortening muscle. This leads to discontinuous derivatives of the force-velocity relation. A discontinuous first derivative proved to give major problems in the method presented. Therefore, the force-velocity relation was described by a continuous function (Table 1; Fig. 3) which approximates the Hill curve from Winters and Stark (1985). The model parameters given in Table 2 were chosen as follows. The maximal isometric joint

-PJ;

P*=MV,,/(Ps+tan

Equivalent to (for ma = f 0.04) 1OONm

30 25 rads-’

Pl).

torques were measured in one subject. Limb inertia was computed from cadaver data (Veeger et al., 1991). The minimal active state A,,, was taken as the resting active state from Hatze (1981, p. 34). A,,, must be taken larger than zero to assure that F, is finite. For the other parameters the values were taken in the range described by Winters and Stark (1988). This model is the basis for the current paper. An inverse model will be developed with approximately equal dynamic relations and with identical parameters.

956

Technical Note

-1

1

0 lengthening

velocity

Fig. 3. Force-velocity relation (Tables 1 and 2).

Test data

A fast goal directed movement was forward simulated using the continuous model in Table 1. The reference input used here was derived from dynamic optimization. The parameters optimized were the input amplitudes at certain time points (nodal points in Pandy et al., 1992). Intermediate inputs were linearly interpolated. Time points were defined with time steps of 0.01 s from t =0.21 until t =0.45 s. Outside this interval a zero input was assumed. The criterion optimized was a sum of squared muscular inputs and of the squared position error starting from t=0.40 s. These reference inputs and the resulting movement are shown in Fig. 4 and will be used as a reference to develop and test the method of inverse dynamics. Inverse muscle model The exact inverse of the model used contains differentiating elements. Such elements do not exist in real life and cannot be simulated exactly. Therefore, the continuous model was approximated by a discrete model (Table 1) that can be inverted. The model consists of three first-order systems in series. For each of these three systems a discrete analogon was determined, assuming a constant input over one time step. Due to this assumption, the discrete model will be exactly equivalent to the continuous model only for infinitely small time steps. Assuming constant inputs, for a linear first-order system the output can be analytically shown to be an exponential function of time. This analytical solution was used to define the discrete analogon for neural excitation (E) and active state (A) dynamics. The third firstorder system consists of the force-velocity relation (CE) and the series elastic element (SE). The equations of CE and SE are strongly nonlinear. The discrete analogon of this system is therefore chosen as a forward simulation which has been performed with a variable time step Kutta-Merson algorithm. This system is only discrete in terms of defined inputs and outputs. The three systems have no direct terms; a change of the input does not instantaneously influence the output. We now have a model in which the input u at sample i-3 first affects the force at sample i. The assumed constant inputs in the three discrete systems all cause a delay of approximately one half time step. This effect will be compensated by shifting (delaying) the estimated input with one and a half time steps. The discrete relations in Table 1 have been inverted. The inverse of the neural excitation and active state dynamics was determined analytically. The CE/SE dynamics are inverted iteratively: the active state during the time interval before sample i is iterated until the desired force at sample i is reached. With this iterative linear extrapolation, the active state is estimated to an accuracy of lO-‘O in 3-7 iterations. In the next paragraph a dynamic physiological constraint on

muscle force will be introduced. Often, the forces computed will equal their lower bound and sometimes they will equal their upper bound. In these cases, it is obvious that the active state will be at its lower or upper bound and the iterative procedure is not necessary. With this inverse model, in three steps, the active state, neural excitation and neural input can be computed. The inverse model requires muscle length and force as input data (length was also an input in the forward muscle model, Figs 1 and 2). The inverse is computed in double precision with a time step of 10 ms. With larger time steps a less precise description of inputs was found. With time steps below 0.5 ms, numerical inaccuracy caused some slight input fluctuations. The muscle model was inverted in three stages for the three first-order submodels. Alternatively, the inversion can be performed in one stage using the iterative procedure now used for the last stage. This yields a procedure similar to that of Pierrynowski and Morrisson (1985). However, we found that the one-stage inversion yielded numerical instability for time steps below 10 ms. Such an instability was also encountered by Caldwel and Chapman (1991) and results from the fact that force responds slowly to the neural input. Thereby, the force at the end of a small time interval depends only weakly on the neural input in this interval. Thus, to obtain the desired forces, large input variations will be necessary, which may cause instability. For this reason the threestage inversion was preferred. Dynamic physiological constraints

Applying the inverse model, sometimes neural inputs u are estimated to be negative or greater than one. This indicates that the desired combination of muscle force/length cannot be met due to muscular dynamics. From the fact that neural inputs range between zero and one, dynamic constraints are derived. These constraints are based on the discrete forward dynamic relations in the right-hand column of Table 1. First, a constraint on the neural excitation E is derived. The lower bound Elb equals E for u = 0 and the upper bound Eub equals E for I(= 1. This constraint is used to determine Alb and &, between which the active state A is constrained. From Alb and Aubr the constraints on muscle forces Irib and Fub are computed. For each sample, these constraints are determined from the states E, A and I,,. These states are determined by force and displacement in the previous samples. So, F,, and F,, describe a time-dependent constraint on muscle force. This physiological constraint reflects muscular dynamics and will be used solving the load sharing problem. Load sharing problem First the desired net torques and the moment arms and musculo-tendon lengths are determined as a function of time (Fig. 2) and then the load sharing problem is solved. A quadratic criterion J is used with weighting factors w on muscle forces F, with constraint C, and satisfying the movement equations: J=?

w,Ff, *=1 C=(FlbjiFj
(2) j=l,2,.

, N).

(3)

Because the two mucles in our model have identical properties equal factors w, are chosen. Since the minimum of a function is not affected by multiplication of the function by any positive value, factors w,= 1 are chosen. Muscular force is constrained between the lower bound Fu, and the upper bound Irut,. For the first sample, Fib and Fub are chosen as the minimal and maximal isometric forces implying a stationary activity U=E=A=O or U=E=A=l. For following samples, the constraints Fib and FYbare derived from the inverse muscle model as described above. With the inverse muscle model also the muscular states and the neural input are computed (see the two previous paragraphs). Because we

951

Technical Note

~

_Jin,,oYe!l

.l

,,,,,,,,,,,,,,

.2

.3

.4 time

.5

.7

.6

(s)

Fig. 4. Inverse solution of the load sharing problem, forces and inputs are determined from the net torque and the joint angle. The test data (net torque, joint angle) and the reference forces have been generated by forward simulation of the reference inputs.

used a quadratic criterion, the load sharing problem could be solved with a quadratic programming method. Figure 4 shows the result of this procedure. The desired net torque is delivered exactly. The muscular forces are slightly lower than the reference forces, which indicates that the inverse solution is more optimal than the reference solution. This, however, does not disqualify the reference solution as this has been optimized with respect to a different criterion. The inputs are quite similar to the reference inputs used to generate the test data. Processing noisy measured data As described above, the method of inverse dynamics was applied to ideal simulated data. Now a scheme will be. presented to evaluate measured trajectories of goal directed movements. Since the method of inverse dynamics requires multiple differentiations it is very important that the trajectory is sufficiently smooth. In measured trajectories, owing to measurement noise, this will generally not be the case. Noise reduction by low-pass filtering would distort signals. Averaging of signals from several movements (ensemble aver-

aging) will also result in a distortion of signals. The idea now is to approximate the trajectory by a smooth function which should describe only the relevant information. In Happee (1993) a linear model, with which the movement trajectory of a goal directed movement can be described, was presented. This model will be used here to generate a smooth approximation of a measured trajectory. The original model contained a third-order filter. To match the order of the (inverse) muscle model, the order was increased to five without adding extra parameters: H(s)=

1 (T~s+l)(T,s+l)(tTis+l)

1 7. s

(4)

This filter is applied to an input signal of three blocks (Fig. 5). The parameters of this input (height and width of the blocks) and the parameter I’, of the filter are estimated from the measured trajectory. The parameters are estimated by a minimization of (the summed square of) the difference between the actual movement and the model response. The iterative minimization is executed with the Gauss-Newton algorithm with a Levenberg-Marquardt term. The filter contains two

Technical Note

958 2

i

1.5

\s

’

x

.5

.C 30”

measured ..._...__.._ mode,

0 -.5

‘;

50

.-E 5

0

:

0

.6

.7

.a

.9

1

ll,lllll,,llll.//,llll~l.l,.,l,,,,,,,,,,,,,,, .2 .3 .4 .6 .I .5 time (s)

.7

.a

.9

1

.l

.2

.3

.4

.5

I

-50 0

n

Fig. 5. Description of an experimental trajectory by the response of a linear model (upper traces) to an input of three blocks (lower trace).

:

EMG (mv)

0

.1

.2

lsj5

.3 tirt-2

.6

.7

.a

.9

1

agonist .4 3 Q .2 .c 0 antagonist .4 zi 0. .2

.s

0 1 20 ? 2

net

4

I

.3

.4

I

1

/

I

.6

.7

.a

.9

torque

0

I-20

0

.l

.2

s time

1

(s)

Fig. 6. Estimated inputs for the experimental trajectory in Fig. 5 and the EMG (upper traces) of the agonist deltoideus pars anterior and the antagonist latissimus dorsi.

Technical Note poles with time constant T1 and one with T,/2. Several other combinations of time constants were tried. The chosen set was found to give the best approximation of data. This choice appeared not to be very critical because the approximation of the data depended only weakly on the set of time constants. The response of the filter was again determined for the estimated parameters; now the filter was used with the outputs position, velocity and acceleration which are desired for the inverse dynamics. As can be seen in Fig. 6, for this smooth trajectory realistic inputs could be estimated with the method of inverse muscular dynamics. DISCUSSION A method to estimate muscular activity for a given movement trajectory has been presented. First, from the trajectory and the external forces, the desired net torques are calculated by using the well-known method of inverse dynamics. Then, a solution of the load sharing problem is determined by minimization of the sum of squared muscle forces and satisfying a dynamic physiological constraint on muscular force. This constraint takes into account the dynamics of the contractile element (CE) and the series elastic element (SE), active state dynamics and neural excitation dynamics. Hannaford et al. (1986) also solved this inverse problem but in a very different way. In fact, they performed a dynamic optimization in which the difference between the experimental trajectory and the model response was minimized. Figures 1 and 2 describe the difference between such forward simulations and our inverse method. Computational

costs

The main advantage of our method lies in its limited computational costs. The computational costs of our inverse method and of a dynamic optimization will now be compared. Our method requires only one simulation of the inverse limb dynamics to compute muscle lengths and desired joint torques as a function of time. The computational cost of this simulation is termed CPUli,,+. For each time step, per muscle, on average six simulations are required to determine bounds on muscle force and to compute iteratively the active state required to obtain the desired force. The cost of one such simulation is CPU,,, (summed over all time steps). The load sharing problem is solved at a cost CPU,,,, for each time step. When we now consider the number of muscles (N) and the number of time steps (NT), then the computational cost for the inverse simulation will be: CPUinve,oc= CPU,i,,,b+6NCPU,,,+NTCPU,oad.

(5)

For a dynamic optimization as performed in Happee (1992a) and to generate the test data, parameters describe input amplitudes at certain time points. Thus, the number of parameters to be estimated equals the number of time points (NT) multiplied by the number of muscles (N). With NT = 25 and N = 2, 50 parameters were estimated. A Gauss-Newton algorithm was used with a Levenberg-Marquardt term. Gradients were determined by applying parameter perturbations. This required as many simulations as the number of parameters. One more simulation,was required with the updated parameter set. In these simulations, all muscles, plus the limb, were simulated together. The computation time for such a simulation will be about CPUlimb+N CPU,,,. Reasonable results were obtained in about 40 iterations (Iter). The computational cost of this dynamic optimization will be about: CPU,, c Iter(N NT + 1)(CPUlimb+ N CPU,,,).

(6)

With N =2 we obtain a cost of CPUlimb+ 12 CPU,., + NT CPU,.,, for the inverse dynamics. For the dynamic optimization, with Iter=40 and NT=25 the cost is

959

2040CPU1i,a + 4080 CPU,,, . So, for dynamic optimization, far more evaluations of limb and muscle dynamics are required. On the other hand, dynamic optimization does not require solving the load sharing problem. Because we used a quadratic criterion, the load sharing problem could be solved with a quadratic programming method. Thus, the CPU,..* was about 50% of the total CPU time. In other words, taking into account the muscular dynamics, the CPU time of the inverse simulation was only doubled. Equations (5) and (6) indicate that the difference in computational cost between inverse simulation and dynamic optimization will increase with more complex systems. A strong effect will follow from the number of muscles N which contributes proportionally to equation (5) and quadratically to equation (6). Furthermore, with dynamic optimization more parameters generally require more iterations. The feasibility of the method of inverse dynamics including muscular dynamics for complex systems is shown in Happee and Van der Helm (1993b), where a shoulder model with 95 muscle elements was simulated in an acceptable CPU time. Experimental

data

In inverse simulations the displacements and the external forces are used as inputs (Fig. 2). These displacements and forces must be recorded in relevant experiments. Limb displacements can be recorded with video systems. Displacements and forces can be recorded with mechanical devices. Since the method of inverse dynamics requires multiple differentiations it is very important that the trajectory is sufficiently smooth. It is therefore important to collect movement data with a high resolution, low noise setup. In this respect mechanical measurement is favored against video. Furthermore, dedicated data processing will be required to generate velocity and acceleration signals of a sufficient smoothness. A general method would be to apply spline functions (Woltring, 1986). For periodic movements Fourier analysis could be employed where smoothness can be ensured by omitting higher frequencies and differentiation becomes a simple operation on the Fourier coefficients. For goal directed movements we proposed a parametric description of the movement trajectory which is based on a triphasic control (Fig. 5). Oprimality Inverse methods can be used to determine optimal inputs for a certain trajectory. However, inverse. methods cannot prove the optimality of a trajectory. The only sure thing about experimental trajectories is that they are realistic. This must be kept in mind whenever employing inverse dynamics in relation to optimality. The dynamic physiological constraint is computed from the state of the system which represents the movement history. The most important constraint is the limited rate of force decrease. When a force is no longer required, our method predicts compensation by antagonists. However, the dynamic constraint does not anticipate on forces required in the future. The criterion is optimized locally (per sample). Figure 4 demonstrates that, for a fast goal directed arm movement, the results of this local optimization agree well with the reference solution which has been derived by dynamic optimization of a static criterion. The criterion minimized consists of the sum of squared muscle forces [equation (2)]. With our method, any other criterion in the muscular forces can be minimized. In dynamic conditions, however, force will poorly reflect energy consumption (Herzog and Leonard, 1991). With the inverse muscle model variables like the active state and the shortening velocity of the contractile element can be computed. This allows a separation of the mechanical work of a muscle into energy released by the series elastic element and energy delivered by the contractile element. Such variables can be

960

Technical Note

included in the criterion to describe the metabolic cost of contraction (Happee and Van der Helm, 1993a). It is concluded that, with the method presented, muscular forces, neural inputs and muscular states can be simulated with a limited computational cost. This makes the method a useful tool, especially for more complex systems. REFERENCES

Audu, M. L. and Davy, D. T. (1985) The influence of muscle model complexity in musculoskeletal motion modelling. .I. biomech. Engng 107, 147-157. Caldwel, G. E. and Chapman, A. E. (1991) The general distribution problem: a physiological solution which includes antagonism. Hum. Mumt Sci. 10, 355-392. Crowninshield, R. D. and Brand, R. A. (1981) A physiologically based criterion of muscle force prediction in locomotion. J. Biomechanics 14, 793-801. _ Dul. J.. Townsend. M. A.. Shiavi. R. and Johnson. G. E. (1984a) Muscular synergism-L’on criteria for load sharing between synergistic muscles. J. Biomechanics 17, 663-673.

Dul, J., Johnson, G. E., Shiavi, R. and Townsend, M. A. (1984b) Muscular synergism-II: a minimum-fatigue criterion for load sharing between synergistic muscles. J. Biomechanics 17, 675-684.

Hannaford, B., Kim, W. S., Lee, S. H. and Stark, L. (1986) Neurological control of head movements: inverse modelling and electromyographic evidence. Math. Biosci. 78, 159-178. Happee, R. (1992a) Time optimality in the control of human movements. Biol. Cybern. 66, 357-366. Happee, R. (1992b) Goal-directed arm movement. I: analysis of EMG records in shoulder and elbow muscles. J. Electromyogr. Kinesiol. 2, 165-178.

Happee, R. (1993) Goal-directed arm movements. II: a kinematic model and its relation to EMG. J. Electromyogr. Kinesiol. 3, 13-23.

Happee, R. and Van der Helm, F. C. T. (1993a) Criteria representing the cost of muscular contraction, used to solve the load sharing problem in inverse dynamic optimizations. J. Biomechanics (submitted). Happee, R. and Van der Helm, F. C. T. (1993b) The control of shoulder muscles in goal directed movements, an inverse dynamic analysis. J. Biomechanics (submitted). Hatze, H. (1981) Myocybemetic control models of skeletal muscle, characteristics and applications. University of South Africa, ISBN 0 869812165. Herzog, W. and Leonard, T. R. (1991) Validation of optimization models that estimate the forces exerted by synergist muscles. J. Biomechanics 24 (Suppl. I), 31-39. Karlsson, D. and Peterson, B. (1992) Towards a model for force predictions in the human shoulder. J. Biomechanics 25, 189-199.

Pandy, M. G., Anderson, F. C. and Hull, D. G. (1992) A parameter optimization approach for the optimal con-

trol of large scale musculoskeletal systems. J. biomech. Engng 114,45@-460.

Patriarco, A. G., Mann, R. W., Simon, S. R. and Mansour, J. M. (1981) An evaluation of the approaches of optimization models in the prediction of muscle forces during human gait. J. Biomechanics 14, 513-525. Pedersen, D. R., Brand, R. A., Cheng, C. and Aurora, J. S. (1987) Direct comparison of muscle force predictions using linear and nonlinear programming. J. biomech. Engng 109, 192-198. Pedotti, A., Krishan, V. V. and Stark, L. (1978) Optimization of muscle-force sequencing in human locomotion. Math. Biosci. 38, 57-76.

Pierrynowski, M. R. and Morisson, J. B. (1985) A physiological model for the evaluation of muscular forces in human locomotion. Math. Biosci. 75, 69-101. Seif Naraghi, A. H. and Winters, J. M. (1990) Optimized strategies for scaling goal-directed dynamic limb movement. In Multiple Muscle Systems: Biomechanics and Movement Organization (Edited by Winters, W. M. and Woo, S. L. Y.). Springer, New York. Van der Helm, F. C. T. (1993a) A finite element musculoskeletal model of the shoulder mechanism. J. Biomechanics (accepted). Van der Helm. F. C. T. (199313)Analvsis of the kinematic and dynamic behaviour of the’shoulder mechanism. J. Biomechanics (accepted). Van der Helm, F. C. T., Veeger, H. E. J., Pronk, G. M., Van der Woude, L. H. V. and Rozendal, R. H. (1992) Geometry parameters for musculoskeletal modelling of the shoulder mechanism. J. Biomechanics 25, 129-144. Veeeer. H. E. J.. Van der Helm. F. C. T.. Van der Woude. LYH: V., Pronk, G. M. and Rozendal, R. H. (1991) Inertia and muscle contraction parameters for musculoskeletal modelling of the shoulder mechanism. J. Biomechanics 24, 615-629.

Winters, J. M. and Stark, L. (1985) Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models. IEEE Trans. biomed. Engng BME-32,826-839.

Winters, J. M. and Stark, L. (1987) Muscle models: what is gained and what is lost by varying model complexity. Biol. Cybern. 55,403420. Winters, J. M. and Stark, L. (1988) Estimated properties of synergistic muscles involved in movements of a variety of human joints. J. Biomechanics 21, 1027-1041. Yamaguchi, G. T., Sawa, A. G. U., Moran, D. W., Fessler, M. J. and Winters, J. M. (1990) A survey of human musculotendon actuator parameters. In Multiple Muscle Systems: Biomechanics and Movement Organization (Edited by Winters, W. M. and Woo, S. L. Y.). Springer, New York. Woltring, H. J. (1986) A Fortran package for generalized, cross-validatory spline smoothing and differentiation. Adv. Engng Software 8, 104-l 13. Zajac, F. E. (1989) Muscle and tendon: properties, models, scaling and application to biomechanics and motor control. Crit. Rev. biomed. Engng 17, 359-411.