2001_Anderson_Pandy_JB.pdf

Journal of Biomechanics 34 (2001) 153}161

1999 ASB Pre-Doctoral Award

Static and dynamic optimization solutions for gait

are practically equivalent

Frank C. Anderson !, * , Marcus G. Pandy !,"

! Department of Mechanical Engineering, The University of Texas at Austin, Austin, Texas, USA

" Department of Kinesiology and Health Education, The University of Texas at Austin, Austin, Texas, USA

Accepted 3 July 2000

Abstract

The proposition that dynamic optimization provides better estimates of muscle forces during gait than static optimization is

examined by comparing a dynamic solution with two static solutions. A 23-degree-of-freedom musculoskeletal model actuated by 54

Hill-type musculotendon units was used to simulate one cycle of normal gait. The dynamic problem was to "nd the muscle excitations

which minimized metabolic energy per unit distance traveled, and which produced a repeatable gait cycle. In the dynamic problem,

activation dynamics was described by a "rst-order di!erential equation. The joint moments predicted by the dynamic solution were

used as input to the static problems. In each static problem, the problem was to "nd the muscle activations which minimized the sum

of muscle activations squared, and which generated the joint moments input from the dynamic solution. In the "rst static problem,

muscles were treated as ideal force generators; in the second, they were constrained by their force}length}velocity properties; and in

both, activation dynamics was neglected. In terms of predicted muscle forces and joint contact forces, the dynamic and static solutions

were remarkably similar. Also, activation dynamics and the force}length}velocity properties of muscle had little in#uence on the static

solutions. Thus, for normal gait, if one can accurately solve the inverse dynamics problem and if one seeks only to estimate muscle

forces, the use of dynamic optimization rather than static optimization is currently not justi"ed. Scenarios in which the use of dynamic

optimization is justi"ed are suggested.

(

Keywords: Dynamic optimization; Static optimization; Gait; Articular contact forces; Muscle forces

1. Introduction

and processing of body segmental kinematics (Patriarco

et al., 1981; Davy and Audu, 1987). Second, while it is

purportedly necessary to account for the underlying

physiological properties of muscle in order to obtain

more accurate estimates of muscle force (Hardt, 1978; An

et al., 1989; Kaufman et al., 1991; Prilutsky et al., 1997),

the time-independent nature of static optimization makes

it relatively di$cult to incorporate muscle physiology

properly. Third, the time-independence of the perfor-

mance criterion required by static optimization may not

permit the objectives of the motor task to be properly

characterized (Hardt, 1978; Pandy et al., 1995). Finally,

analyses based on an inverse dynamics approach may

not be appropriate for explaining muscle coordination

principles (Kautz et al., 2000; Zajac, 1993).

Dynamic optimization, a forward dynamics approach,

is not subject to these criticisms, and ostensibly can

provide more realistic estimates of muscle force. Dy-

namic optimization integrates system dynamics into the

Static optimization, an inverse dynamics approach,

has been used extensively to estimate in vivo muscle

forces during gait (Seireg and Arvikar, 1975; Pedotti et

al., 1978; Hardt, 1978; Crowninshield et al., 1978; Crow-

ninshield and Brand, 1981; Patriarco et al., 1981; RoK hrle

et al., 1984; Brand et al., 1986; Pedersen et al., 1997;

Glitsch and Baumann, 1997). Static models are computa-

tionally e$cient, have allowed full three-dimensional

motion, and have generally incorporated 30 or more

muscles per leg. However, static optimization has been

criticized on several points. First, the inverse dynamics

problem is highly dependent on the accurate collection

* Correspondence address. Biomechanical Engineering Division, De-

partment of Mechanical Engineering, Stanford University, Stanford,

California 94305-3030, USA. Fax: 650-725-1587.

E-mail address: fca@stanford.edu (F.C. Anderson).

PII: S 0 0 2 1 - 9 2 9 0 ( 0 0 ) 0 0 1 5 5 - X

154

F.C. Anderson, M.G. Pandy / Journal of Biomechanics 34 (2001) 153 } 161

solution process: quantities like muscle forces and the

performance criterion are treated as time-dependent state

variables whose behavior is governed by sets of di!eren-

tial equations (Hatze, 1976). Ideally, the di!erential equa-

tions accurately represent the underlying physiological

properties of the system, and so the time histories of

muscle forces predicted are consistent with those forces

which could naturally arise during movement.

Unfortunately, dynamic optimization incurs so much

computational expense that relatively few dynamic solu-

tions for gait have been found (Chow and Jacobson,

1971; Davy and Audu, 1987; Yamaguchi and Zajac, 1990;

Tashman et al., 1995). Further, for gait, this approach has

required that the dynamic models be simpli"ed to such

an extent that it has been di$cult to ascertain whether its

computational expense is justi"ed. While Davy and

Audu (1987) solved static and dynamic optimization

problems for gait as part of the same study and found

di!erences between the two solutions, it is di$cult to

draw conclusions based on these di!erences because (1)

the authors simulated only the swing phase of gait, (2) the

model was con"ned to the sagittal plane and included

only nine muscles, and (3) the joint moments used as

input to the static optimization problem were di!erent

from those predicted by the dynamic optimization solu-

tion.

The aim of this work is to assess whether or not

dynamic optimization provides more realistic estimates

of muscle force than static optimization. In particular, we

address the following questions:

When estimating in vivo muscle and joint contact

forces during gait

2. Methods

2.1. Musculoskeletal model

The body was modeled as a 10-segment, 23-degree-of-

freedom linkage (Fig. 1). The inertial properties of the

segments were based on anthropometric measures ob-

tained from "ve healthy adult males (age 26$3 yr,

height 177$3 cm, and mass 70.1$7.8 kg) (McConville

et al., 1980). Interactions of the feet with the ground

were modeled using a series of spring-damper units

distributed under the sole of each foot (Anderson and

Pandy, 1999).

The model was controlled by a total of 54 musculoten-

dinous actuators (Fig. 2). Each leg was controlled by 24

actuators. Relative movements of the pelvis and upper

body were controlled by six abdominal and back ac-

tuators. Each actuator was modeled as a three-element,

Hill-type muscle in series with tendon (Zajac, 1989). The

muscle parameters, as well as the origin and insertion

sites, were based on data reported by Delp (1990). The

actions of ligaments were used to prevent anatomically

infeasible joint angles from arising during a simulation.

Each ligament was modeled as an exponential curve

which provided a restoring torque about a joint axis as

a function of joint angle (Davy and Audu, 1987). For

(1) Is the large computational cost of dynamic optimiza-

tion justi"ed?

(2) Is it important to account for muscle physiology,

namely activation dynamics and the force}length}

velocity properties of muscle?

The basis for the assessment is a comparison of a

previously attained dynamic optimization solution for

gait (Anderson, 1999; Anderson and Pandy, submitted)

with two analogous static optimization solutions, one in

which the force}length}velocity properties of muscle

were taken into account and one in which they were

neglected. The musculoskeletal model used allowed

three-dimensional motion of the body parts and simula-

tion of a full gait cycle. Furthermore, the model was

actuated by a number of muscles which approaches that

used in earlier static optimization models. Signi"cantly,

the joint torques predicted by the dynamic optimization

solution were used as input to the static optimization

problems. In this way, the experimental errors endemic

to solving the inverse dynamics problem were eliminated,

and the results of the static and dynamic optimization

solutions could be compared more directly.

Fig. 1. Model of the body. The "rst six-degrees-of-freedom were used to

de"ne the position and orientation of the pelvis relative to the ground.

The remaining nine segments branched out in an open chain from the

pelvis. The head, arms, and torso were represented as a single rigid body

which articulated with the pelvis via a three-dof ball-and-socket joint

located at approximately the third lumbar vertebra. Each hip was

modeled as a three-dof ball-and-socket joint, each knee as a one-dof

hinge joint, each ankle-subtalar joint as a universal joint with a single

joint center, and each metatarsal joint as a hinge joint. The directions of

the knee, ankle, subtalar, and metatarsal joint axes were anatomical

and based on in vivo and cadaveric measurements (Inman, 1976;

Anderson, 1999; Anderson and Pandy, 1999).

F.C. Anderson, M.G. Pandy / Journal of Biomechanics 34 (2001) 153 } 161

155

served as controls in the dynamic optimization problem:

, 2 , u

1 m

(0), m "1, 2 , 54,

(1)

where u i m

is the i th muscle excitation node for muscle

(0) is the activation level of muscle m at t "0in

the simulation. The dependence of muscle activation

level on excitation level was described by a "rst-order

di!erential equation

5 "( u

2 ! ua )/

q 3*4%

#( u ! a )/

q &!--

, a " a

( t ), u " u

( t ),

(2)

where a

( t ) is the activation level of muscle m at time t ,

( t ) is the excitation level of muscle m at time t , and

q 3*4%

(200 ms) are the rise and decay

constants of muscle activation (Pandy et al., 1992).

The dynamic optimization problem was minimally

constrained: only those constraints necessary to ensure

a repeatable gait cycle were imposed. Speci"cally, ter-

minal boundary constraints were applied to the joint

angular displacements, joint angular velocities, muscle

activation levels, and muscle excitation levels so that the

states at the end of the simulation were approximately

equal to those at t "0. Because only one-half of the gait

cycle was simulated, terminal states for the left-hand side

of the body were constrained to equal the initial states for

the right-hand side of the body, and vice-versa:

(22 ms) and

q &!--

Fig. 2. Muscles in the model. Abbreviations used for the muscles are as

follows: (ERCSPN) erector spinae; (EXTOBL) external abdominal ob-

liques; (INTOBL) internal abdominal obliques; (ILPSO) iliopsoas;

(ADLB) adductor longus brevis; (ADM) adductor magnus; (GMEDA)

anterior gluteus medius and anterior gluteus minimus; (GMEDP) pos-

terior gluteus medius and posterior gluteus minimus; (GMAXM) me-

dial gluteus maximus; (GMAXL) lateral gluteus maximus; (TFL) tensor

fasciae latae; (SAR) sartorius; (GRA) gracilis; (HAMS) semimem-

branosus, semitendinosus, and biceps femoris long head; (RF) rectus

femoris; (VAS) vastus medialis, vastus intermedius, and vastus lateralis;

(BFSH) biceps femoris short head; (GAS) gastrocnemius; (SOL) soleus;

(PFEV) peroneus brevis and peroneus longus; (DFEV) peroneus tertius

and extensor digitorum; (DFIN) tibialis anterior and extensor hallucis

longus; (PFIN) tibialis posterior, #exor digitorum longus, and #exor

hallucis longus. Muscles included in the model but not shown in the

diagram are: (PIRI) piriformis; (PECT) pectinius; (FDH) #exor

digitorum longus/brevis and #exor hallucis longus/brevis and; (EDH)

extensor digitorum longus/brevis and extensor hallucis longus/brevis.

( t

)" q

(0), i "4, 2 , 23,

(3)

5 #-

( t

)" q

5 i

(0), i "4, 2 , 23,

(4)

# m

( t

)" a

(0), m "1, 2 , 54,

(5)

u 15,#-

" u m

, m "1, 2 , 54,

(6)

( t ) are the i th generalized coordinate

and generalized velocity of the model, respectively, a

( t ) and q

5 i

details concerning the model refer to Anderson (1999)

and Anderson and Pandy (1999).

( t )

are the

15th and 1st control node for muscle m , respectively. The

superscript cl indicates the contralateral quantity. For

example, in Eq. (5) the contralateral quantity is the ac-

tivation level of the contralateral muscle. Note that con-

straints were not placed on the generalized coordinates

or velocities relating to the translation of the pelvis (i.e,

i "1, 2 , 3 in Eqs. (3) and (4)).

There is much experimental evidence to suggest that

people select walking speeds which minimize the meta-

bolic energy expended per unit distance traveled (Ral-

ston, 1958,1976; Corcoran and Brengelmann, 1970;

Finley and Cody, 1970). The performance criterion for

the dynamic optimization problem was therefore ex-

pressed as follows

1 m

and u

2.2. Dynamic optimization problem

A "xed-"nal-time dynamic optimization problem was

solved for one cycle of normal gait. Because the gait cycle

was assumed to be bilaterally symmetric, it was only

necessary to simulate one-half of the gait cycle. The initial

states were obtained by averaging kinematic and force

platform data obtained from "ve subjects. The "nal time

for the dynamic optimization problem was "xed to the

average time taken by the "ve subjects to complete half of

a gait cycle, which was 0.56 s.

The excitation history for each muscle was para-

meterized by 15 nodal points which were allowed to vary

continuously between zero and one, zero indicating no

excitation and one indicating maximal excitation. The

nodal points representing the muscle excitations, to-

gether with the initial activation level for each muscle,

t &

B Q # 5 +

m /0

( A Q

# M Q

# S Q

#=Q

)

d t

J "

(7)

( t

)! X

(0)

m and a

where q

is the activation level of muscle m , and u

156

F.C. Anderson, M.G. Pandy / Journal of Biomechanics 34 (2001) 153 } 161

where the numerator is a computation of the total meta-

bolic energy consumed and the denominator is the

change in position of the center of mass in the direction of

progression. The equations used to estimate metabolic

energy from the internal states of the muscles were de-

veloped based largely on the work of Davy and Audu

(1987), Hatze and Buys (1977), and Mommaerts (1969).

B Q is the basal metabolic rate of the whole body. For

muscle m , A Q

) is the activation level of muscle m at a dis-

crete time step designated by t

( t

. The simulation interval

(i.e., 0.0) t )0.56) was broken into 172 segments, and

a static optimization problem was solved at each t

Activation dynamics was neglected.

To assess the importance of incorporating muscle

physiology, one static optimization problem was solved

neglecting the force}length}velocity properties of muscle,

and a second was solved with these properties included.

In both cases, at each discrete time step, t

is the activation heat rate, M Q

is the

maintenance heat rate, S Q

is the shortening heat rate, and

, the force

output of each muscle was computed from its activation

level, a

is the mechanical work performed (Anderson, 1999).

Hyper-extension of the joints was penalized by appen-

ding the following function to the performance criterion:

). In the `non-physiologicala case, each muscle

was assumed to be an ideal force generator

( t

)" w

t &

1 +

j /1

-*' j

( t )

d t ,

(8)

( t

)" a

( t

) F

o m

(11)

its maximum isometric force. In the `physiologicala case,

the force generated by a muscle was constrained by its

force}length}velocity properties:

( t

) is the force generated by muscle m and F

o m

( t ) is the torque applied by the ligaments

about joint axis j at time t , and w (0.001) is a weight

factor. The performance criterion was therefore aug-

mented

-*' j

( t

)" a

( t

) f ( F

o m

, l

, v

(12)

@ " P

t &

B Q # 5 +

m /0

( A Q

# M Q

# S Q

#=Q

)

d t

) is the general

force}length}velocity surface for muscle assumed in the

musculoskeletal model, and l

o m

, l

, v

# /

( t

)! X

(0)

the

shortening velocity of muscle m (Zajac, 1989). Muscle

lengths and shortening velocities were computed from

the joint angular displacements and velocities predicted

by the dynamic optimization solution, with tendon com-

pliance and muscle pennation angle taken into account.

The force outputs of the muscles were constrained to

produce the muscular joint moments predicted by the

dynamic optimization solution

is the length and v

(9)

Thus, the dynamic optimization problem was to "nd

values of the controls (Eq. (1)) such that the constraints

(Eqs. (3)}(6)) were satis"ed and the performance criterion

(Eq. (9)) was minimized. Note that the dynamic equations

for the skeleton (not shown), for the muscle forces (not

shown), and for the muscle activations (Eq. (2)) were also

constraints in the dynamic problem and were enforced

implicitly by numerically integrating these equations for-

ward in time during the simulation. The dynamic optim-

ization problem was solved using a parallel parameter

optimization algorithm (Pandy et al., 1992; Anderson

et al., 1995).

5 +

m /1

( t

) r

m , j

( t

)"¹

( t

), j "7, 2 , 23, i "1, 2 , 172,

(13)

)is

the moment arm of muscle m about the j th joint axis, and

( t

) is the force generated by muscle m , r

m , j

( t

2.3. Static optimization problems

) is the muscular joint moment acting about the j th

joint axis as predicted by the dynamic optimization solu-

tion. Note that Eq. (13) begins with j "7, the

#exion}extension axis of the back. Because the "rst six

degrees-of-freedom concerned the displacement and ori-

entation of the pelvis relative to the ground (see Fig. 1),

no muscular moments were exerted about these general-

ized coordinates.

The performance criterion was to minimize activation

squared, summed across all muscles (Kaufman et al.,

1991; Crowninshield and Brand, 1981):

( t

So that the static and dynamic solutions could be

compared more directly, the joint angular displacements,

joint angular velocities, and muscular joint moments

predicted by the dynamic optimization solution were

input into the static optimization problem, as though

these quantities came from experimental measures. In

this way, imprecision in the collection and processing of

experimental data and errors in the estimation of body

anthropometry were eliminated.

Following the example of Kaufman et al. (1991), in

order to account for the force}length}velocity properties

of each muscle, the muscle activation levels were used as

the controls in the static optimization problems:

" 5 +

m /1

( a

( t

))

, i "1, 2 , 172.

(14)

In the non-physiological case (Eq. (11)), note that

muscle activation is equal to muscle stress multiplied by

( t

), m "1, 2 , 54, i "1, 2 , 172,

(10)

where a

where F

where ¹

where the function f ( F

where F

F.C. Anderson, M.G. Pandy / Journal of Biomechanics 34 (2001) 153 } 161

157

some proportionality constant

0 m

" k

PCSA ,

(15)

/PCSA is muscle stress, PCSA is physiological

cross-sectional area, and k is a constant. In the physiolo-

gical case (Eq. (12)), muscle activation is approximately

equal to muscle stress multiplied by some propor-

tionality constant, provided a muscle operates on the #at

portion of its force-length curve and at small contraction

velocities:

) + k

PCSA

when l

+ l

0 m

and

f ( F

0 m

, l

, v

+0.0.

(16)

is the optimal muscle-"ber length.

Thus, the non-physiological static optimization prob-

lem was to "nd an activation level for each muscle in the

model such that the performance criterion (Eq. (14)) was

minimized, the joint torques predicted by the dynamic

optimization solution (Eq. (13)) were generated at each

joint, and each muscle was treated as an ideal force

generator (Eq. (11)). The physiological static optimiza-

tion problem was the same, except that each muscle was

constrained to act in accordance with its

force}length}velocity properties (Eq. (12)). Both static

optimization problems were solved using a gradient-

based sequential quadratic programming algorithm

(Pandy et al., 1992).

0 m

Fig. 3. The net muscle moments applied about the axes of the hip joint

are shown as predicted by the dynamic optimization solution (thick

black lines) and as measured by Crowninshield et al. (1978) (circles),

Cappozzo et al. (1975) (triangles), Patriarco et al. (1981) (crosses), and

Inman et al. (1981) (squares).

3. Results

The joint moments predicted by the dynamic optim-

ization solution (and used as input to the static optimiza-

tion problem) were similar to those reported in the

literature (Figs. 3 and 4) (Cappozzo et al., 1975; Crownin-

shield et al., 1978; Patriarco et al., 1981; Inman et al.,

1981). Joint moments about the axes of the back, subta-

lar, and metatarsal joints were also predicted by the

dynamic solution but are not reported here. Discontinui-

ties in the joint torques are present because the terminal

boundary constraints imposed in the dynamic optimiza-

tion solution were not precisely met.

In general, the forces predicted by the dynamic and

static solutions were very similar for many muscles (Figs

5 and 6: ILPSO, GMAXL, GMEDA, VAS, GAS, SOL,

TFL, BFSH, PFEV). However, there were two notable

exceptions. Following heelstrike, the dynamic solution

predicted signi"cant activity in the medial portion of

gluteus maximus, whereas the static solutions predicted

very little activity in this muscle (Fig. 5: GMAXM,

0}30%). On the other hand, the static solutions predicted

signi"cant activity in hamstrings, whereas the dynamic

Fig. 4. The net muscle moments applied about the axes of the knee and

ankle joints are shown as predicted by the dynamic optimization

solution (thick black lines) and as measured by Cappozzo et al. (1975)

(triangles) and Inman et al. (1981) (squares).

solution predicted much less activity in this muscle

(Fig. 5: HAMS, 0}30%).

The resultant joint contact forces at the hip, knee, and

ankle predicted by the dynamic and both static optimiza-

tion solutions were very similar to each other and some-

what similar to contact forces reported by Crowninshield

and Brand (1981) and Hardt (1978) (Fig. 7). For both the

where F

where l

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika: