Physiological prediction of muscle forces—I. Theoretical formulation

A’euroscience Vol. 40, No. 3, pp. 781-792, Printed in Great Britain

1991

0306-4522/91 $3.00 + 0.00 Pergamon Press plc 0 1991 IBRO

PHYSIOLOGICAL PREDICTION OF MUSCLE FORCES-I. THEORETICAL FORMULATION K. R.

KAUFMAN,*

K.-N. AN,? W. J. LITCHY~ and E. Y. S. Cmot$

*Motion Analysis Laboratory, Children’s Hospital, San Diego, CA 92123, U.S.A. tBiomechanics Laboratory, Department of Orthopedics and SSection of Electromyography, Division of Clinical Neurophysiology, Department of Neurology, Mayo Clinic/Mayo Foundation, Rochester, MN 55905, U.S.A. Ahatract-A physiological model for predicting muscle forces is described. Rigid-body mechanics and musculoskeletal physiology are used to describe the dynamics of the segment model and muscle model. Unknown muscle and joint contact forces outnumber the equilibrium equations resulting in an indeterminate problem. Mathematical optimization is utilized to resolve the indeterminacy. The modeling procedure relies entirely on established physiological principles. Data describing the muscle anatomy and body structures are included. A model defining the force-length-velocity-activation relationship of a muscle is adopted. The force a muscle produces is assumed to be proportional to its maximum stress, physiological cross-sectional area, activation, and its functional configurations including the muscle architecture, muscle length, contracting velocity, and passive tension. These factors are incorporated into inequality equations which limit the force for each muscle. Minimal muscular activation is forwarded as the optimization criterion for muscle force determination.

direct measurement of muscular forces is invasive and currently impractical, formulations to mathematically predict muscle forces typically involve modeling the musculoskeletal system as a number of interconnected rigid body segments. Kinematic and kinetic measurements are made to predict resultant actions between segments as a result of the motion. This is called the inverse dynamics problem.*’ Since there are often more muscles present than are required to produce any displacement pattern from kinematic data, the classical equations of kinetic analysis do not permit a unique solution of the muscular forces crossing the joints. There are more unknown forces than can be determined in the equipollence relations between resultants and individual muscle forces. This well recognized difficulty in actually determining the muscle forces is the so-called indeterminant problem.67 Two approaches have been used to deal with this problem. The first may be called the reduction method. The goal of this method is to reduce an initially indeterminate problem to one that is determinate. This method was used to get an approximate solution at the hip during gait6’ and to assign knee forces during dynamic situations. 61The second approach has been to use optimization methods in which the equipollence relations at any point are solved on the basis of optimizing some criteria, called the objective function. The underlying assumption behind the optimization method is that the central nervous system performs Since

the control of muscle action in an optimal manner, i.e. by minimizing some performance criteria or cost function, such as time, mechanical work, total energy expended, or tracking error, or by maximizing mechanical efficiency. The theory of optimization, which has been used extensively in control theory, has also been applied to motor behavior.2’*23*2s*3’~4i~42 The proper objective function is not known a priori, but must be theorized. The selection of different optimization criteria is based upon either physiological or empirical considerations. A particular objective function may be chosen because of a proposed physiological basis for belief in its validity to adequately predict the distribution of muscle force.s.3*26 The hypothesis that the force load is distributed among the muscles in order to minimize energy consumption4i or minimize individual muscle stress2’*22*25,26,28 is such an example. Alternatively, the choice of optimization criteria may be governed by empirical considerations. A particular objective function may, for instance, agree well with the measured temporal activity of the muscles or with expected quantitative levels of muscle force.” Unfortunately, the results of optimization methods are not always physiologically consistent. One way to improve the predictions of muscle forces with the optimization method is to formulate additional constraints. An inequality constraint that limits the stress in each muscle has resulted in improved predictions.3~25~26~29~64~68,74

§To whom correspondence should be addressed. Abbreviations: GM, m. gastrocnemius, caput mediale; ECSA, physiological cross-sectional area; SE, series elastic element; SM, m. semimembranosus.

In the present paper, an optimization approach is developed to predict individual muscle forces for a given activity. The methods used are an extension of previous muscle modeling efforts3 Physiological

COORDINATE

SYSTEMS

-IlIP

JOINT

ANKLE JOINT

Fig. 1. Mechanical model of the lower extremity. The exact procedures to define each of the anatomically based pelvic, femoral, tibial, and foot reference frames are described in the text.

characteristics of the muscle have been represented mathematically and included in the model development. These features include geometrical, morphological, and force-length-velocity-activation properties. Constraints on the muscle forces in the mathematical formulation of the optimization method are formulated using these physiological phenomena. A unique method for solving the individual muscle forces is presented.

MODEL

DEVELOPMENT

The approach used is to develop a biomechanical model to predict muscle forces. The model is hierarchical and incorporates rigid-body mechanics and musculoskeletal physiology. An anatomical model represents the length, moment arms, and velocities of the muscles. On the next level, the physiological aspects of muscle are incorporated. Anatomical model The lower extremity is modeled as a system shown in Fig. 1. This mechanical system is composed of four rigid body segments, Bi (i = 0, 1,2,3), corresponding to the pelvis (B,), the thigh (B,), the shank (B,), and the foot (B,). The foot is assumed to move with the shank as one integrated body. In order to establish a mathematically workable model, four Cartesian co-

ordinate systems OxqJl;_ii(k = 0. 1,2, 3) are established as shown in Fig. 1. A segmental axis system IS established for each body segment. These anatomically based axis systems are fixed in each body segment, B, (i = 0, 1,2,3), and move with it. These coordinate systems are used to define the attachment points of tendons and ligaments as well as joint orientation. Kinematics and kinetics. The intersegmental loading describes the external force system acting on the limb. The force and moment is estimated from kinematic and kinetic information by first finding the position and orientation of the body segment as a function of time. This is expressed in terms of the center of mass and the transformation matrix with respect to a fixed laboratory frame. The linear and angular velocities and accelerations are subsequently calculated by successive differentiation.” The kinematic data are then integrated with the forces and moments acting on the limb to derive the equations of motion for this system of interconnected rigid bodies. Through free-body analysis, a set of six force and moment equilibrium equations are obtained to describe the dynamics of the skeletal system. These equations represent the inverse dynamics problem which is solved based on Eulerian angles and Newtonian equations of motion to yield the intersegmental force, FP, and moment, MP. The details of these calculations are presented elsewhere.‘” The intersegmental resultants describe the net effect of the action of the muscles, ligaments, and the articular contact at the joint. Determination of the dynamic equilibrium of the musculoskeletal model requires calculation of the time-dependent external forces (intersegmental loading) which are then balanced by the time-dependent internal (muscle and joint) forces. This balance of external forces and internal forces can be formulated as:

~~,llM((i,xii),+s=~‘;

k=l,3,

(lb)

where Pp, tip = intersegmental force and moment due to the external force system; @’ = muscle force of ith muscle; 9, fii = joint constraint force and moment; & = unit force vector of ith muscle; and +, = location of ith muscle insertion with respect to the joint center. Muscle unit force vector. The lines of action of musculotendinous units may be modeled in two different ways.49 These two models are referred to as the straight line model and the centroid line model. The straight line model requires the points of muscle attachment on the body segments proximal and distal to the given joint and then assumes that the force transmission by the muscle acts along a straight line connecting the points of attachment. The centroid

Physiological prediction of muscle forces-I

line model requires the locus of transverse crosssectional centroids for any joint configuration and then assumes that muscle force is tangent to the three-dimensional locus of the centroid line. Although the centroid line model correctly models the line of muscle action as a curved path, it has a number of disadvantages. Firstly, a large amount of data are needed to represent a single muscle in any joint configuration. These data are available for only a few muscles. Secondly, it is difficult to define a transverse cross-section of an irregular anatomical structure, such as a muscle. This is particularly troublesome for muscles with broad attachments and unusual shapes, e.g. gluteus medius and gluteus minimus. Thirdly, the model is not easy to use when the joint configuration changes, especially if a shift occurs in the locus of the muscle’s line of action with respect to the joint center. Fourthly, curved centroid lines obtained from sectioned cadaver specimens or in vioo computed tomography or magnetic resonance images may not accurately represent actual data during muscle contractions. Muscles under tension probably follow a more direct course between the points of attachment. For these reasons, a straight line assumption is used to calculate the time-varying muscle force vectors. This approach is reasonable on the basis of a study in which a comparison between skeletal muscle length calculated using centroidal lines versus straight lines in force analysis of musculoskeletal systems was presented.49 For some muscles, a single straight-line model between the anatomical origin and insertion is not adequate. In this case, effective straight-line segments which pass through “effective” origin and/ or insertion points are used. Origin and insertion coordinate data’8*65 are used and transformed to the coordinate system in this study. The origins and insertions of the muscles are scaled to the limb dimensions of the individual whose motion patterns are under investigation.58 Each muscle is given a three-dimensional coordinate expressed in one of four Cartesian reference frames fixed to the pelvis, femur, tibia, or foot. The muscle origins and insertions, and the position of the knee joint center are in turn transformed to the coordinate system attached to the tibia. To quantify the force and moment effects of the muscles let i, = x,i + y, j + z,Z;

(2)

be the location of the origin and ii = xii + y, j + zii;

(3)

be the location of the insertion of the muscle relative to a localized coordinate system fixed at the joint center, where i, j, and c are unit vectors directed along the positive x, y, and z axes of the tibia, respectively. The muscle force, pi, acting on the tibia can now be expressed as P, = I#Jii,

(4)

783

where 1PiI is the magnitude of the muscle force. Thus, using the straight-line muscle model, * i, - fi

fi=g=(i,_ri,

where Zi is the unit muscle force vector for the ith muscle. Muscle moment arm. Moment is defined as the force times the perpendicular distance from the line of action for the force to the axis of rotation. The axis of rotation in this study is assumed to be the joint center. The time-varying muscle moment may be expressed as the moment generated by a unit force. Mathematically, ril = ii x i,

(6)

where I& is the moment associated with the ith muscle and ii is the position vector from the joint center of rotation to the muscle insertion. Muscle morphometrics

The essential features required to describe the morphology of muscle are its physiological crosssectional area (PCSA), architecture, fiber type, and stress limit. Physiological cross-sectional area. The moment arm of a muscle merely indicates the efficiency of the muscle for rotation of the body segment about that particular joint axis. The force that can be actively generated by each muscle is proportional to the PCSA.2.‘7 The PCSA is obtained by dividing the muscle volume by its true fiber length. The rationale of this definition is simply that the cross-sectional area of a muscle is proportional to the number of its fibers and the individual muscle fiber is the basic element which generates the active tension. The PCSA for lower extremity muscles is given in Table 1. Muscle architecture. Nicolaus Stenos3 described the geometry of muscles in 1667 with remarkable insight. His diagrams gave a clear picture of the way a series Table 1. Muscle architecture

Muscle RF VM VI VL TFL SM ST GRA SAR BFLH BFSH M GAST L GAST

PCSA* (cm*)

Fiber* length (cm)

Muscle* length (cm)

Index of architecture

12.73 21.13 22.30 30.57 6.50 16.87 5.35 1.77 1.70 7.44 5.39 16.22 16.22

6.60 7.03 6.83 6.57 14.00 6.27 15.80 27.70 45.50 8.53 13.93 3.53 5.07

31.60 33.53 32.87 32.43 54.00 26.23 31.70 33.53 50.30 34.20 27.07 24.83 21.67

0.209 0.210 0.208 0.203 0.259 0.239 0.498 0.826 0.905 0.249 0.515 0.142 0.234

*Source: Wickiewicz et al.,” except TFL lengths from Mikosz.r9 Biceps femoris and gastrocnemius PCSA from Wickiewicz et a/.*’ proportioned according to data of Boichut and Valentin.16

K. R. KAUFMANet ul.

784

A Tendon plate

I,=

Lf

-=

Fiber

indexof amhitfscture

htl Fig. 2. Muscle architecture. (A) Three-dimensional arrangement of muscle fibers and tendon plates. (B) Schematic representation of muscle architecture. (C) Important architectural parameters of muscle fiber length (&) and muscle belly length (L,) define the index of architecture, i,.

of muscle fibers join a tendon, and the way in which the tendon becomes progressively thicker as it accepts more and more fibers. Brand” confirmed that the length of fibers is surprisingly constant throughout each individual muscle. The generally fusiform shape of many muscles gives a false impression of unequal fiber lengths. Looking at muscles in situ, a few fibers might appear to extend from one end of muscle to the other, and some appear to reach only part way. However, when the muscle is displayed with the fibers aligned, most of the fibers are about equal in length and uniform in thickness. Furthermore, they all reach from the bone or tendon of origin to the tendon of insertion. A schematic representation of these concepts is given in Fig. 2. The muscle fibers he at an angle to the direction of induced motion between two tendon plates. Two measurements can be taken at muscle optimum length. The length of a single fiber is defined as the muscle fiber length. The distance from the proximal to distal tendon is defined as the muscle belly length. The ratio of the muscle fiber length to the muscle belly length at muscle optimum length is defined as the index of architecture.*’ Since muscles are obviously three-dimensional structures, the above classification represents a simplification. Yet, the index of architecture, i,, is an important parameter to describe the relative fiber length of muscles. All of these factors for the human lower limb musculature are presented in Table 1. The fiber length is calculated as the average length of the fiber bundles within each muscle corrected to a sarcomere length of 2.2pm. Muscle fiber types. In the early 17th century, anatomists noted a distinct difference between dark red and lighter “white” skeletal muscles. Since then, it has become a well established fact that mammalian skeletal muscle fibers can be subdivided into two distinct classes-Type I and II-based on their metabolic

and contractile characteristics.‘” Type 1 muscle tibers contract and relax slowly (slow twitch), utilize oxidative metabolism, are more resistant to fatigue and are recruited at lower force levels. Type II muscle tibers are distinguished by their ability to contract and relax rapidly (fast twitch), tend to be more anaerobic, are more fatigable, and are recruited at higher force levels. Individual muscles typically contain different proportions of these two fiber types. Data on the fiber type distribution of human muscles” are presented in Table 2. No data were available in the literature on the fiber type distribution for the tensor fasciae latae, semimembranosus, semitendinosus, and gracilis. Hence, due to the absence of data, it was assumed that these muscles had an equal proportion of fast and slow twitch fibers. Edgerton et al.‘2 also examined the fiber type populations of selected human leg muscles and demonstrated that the gastrocnemius contained 50% slow twitch fibers, vastus lateralis--32% and vastus intermedius47%. These proportions generally agree with those of Johnson rt a1.j” Muscle stress. The relationship between the force a muscle can generate and its cross-sectional area can be expressed as a constant, (T, the muscle stress limit. It is reasonable to assume the existence of such a constraint since the tensile strength of muscle fibers is limited. This constraint expresses conditions which will prevent tendon rupture and avulsion of the tendon-bone connection. Many authors have expressed the stress developed by a muscle in terms of the ratio of maximum isometric force to PCSA. It is clear, however, that this absolute muscle stress does not have much meaning unless the conditions under which it is measured are rigidly specified. Firstly, the force developed depends upon the length at which the muscle is studied. Secondly, the velocity of contraction will affect the total force produced. Thirdly, PCSA needs to be based on the current definition which accounts for obliquity of the muscle fibers. It is unfortunate that in some instances the tension per unit cross-sectional area has been estimated from the muscle length instead of the fiber length. This would usually lead to an overestimation of the intrinsic strength of the contractile material because, with few exceptions, the overall mammalian muscle length is greater than the lengths of its fibers.

Table 2. Muscle fiber distribution Muscle RF VM VI VL SAR BFLH BFSH M GAST L GAST

Type I (%)*

Type II (%)*

0.381 0.526 0.475 0.424 0.496 0.669 0.669 0.508 0.469

0.619 0.474 0.538 0.602 0.504 0.331 0.331 0.492 0.531

*Source: Johnson er ~1.~

Physioldgical prediction of muscle forces-I

785

tension developed in a tetanus is maximum at the optimum length and decreases with greater and lesser lengths. Several mathematical models have been proposed to describe the length-tension relationSome of these models ship 68.13,16,36,3812.45.63,68,80,83,85,88 involve torque-angle relationships.45,83 The joint torque-angle relationship is a complex function of the length-tension property of muscle and the muscle 100 N/cm*. The scatter of these data makes it imperamoment arms which vary with position. The muscle tive that any model developed to calculate joint and muscle forces has the flexibility to adapt to advances lever arm is a complicated and nonlinear function of the joint angle and thus clouds the true lengthin the field of muscle physiology. This is, fortunately, tension relationship. Other models consider the one of the strong points of describing the muscles as length-tension relationship to be either parabolic,8~8s dimensionless vectors, which is one of the goals of linear,7,68a third-degree this project. The obvious discrepancy in the value of gaussian,6s42asymmetric,36~63~88 polynomial,80 or viscoelastic.‘3 the absolute stress for human skeletal muscle obtained Despite all of these attempts to model the muscle by the above investigators is not germane to the active length-tension relationship, only a few mathepresent study. It is only important to know the upper matical models have attempted to account for the limit of stress which a muscle is capable of sustaining. effects of differing muscle architecture. While developFor this study, a value of 100 N/cm* will be used for ing the length-tension relationship for whole muscle, muscle maximum stress since it represents the upper Woittiez et a1.85accurately accounted for differing boundary reported by other authors. It is assumed muscle architecture. Overall, the force output was a that the maximum stress is the same for each muscle. function of the inclination of the muscle fibers. The The available evidence indicates that there is no relative fiber length of muscles, an important aspect marked difference in the intrinsic strength of the contractile material of fast and slow muscle fibers.24~37~76of architecture, was characterized by the index of architecture, i,. However, the mathematical function Muscle model utilized by Woittiez et al.85was a parabola. The real shape of the length-tension relationship is more comOne of the major problems in modeling musculoplex.“O,‘*Deviation from a parabola may occur and skeletal motion is that of describing the actuators, incorporation of a more complex shape into muscle i.e. the muscles, with appropriate models. The problem models is indicated.& The active length-tension curve is developing a muscle model that is phenomenoof whole muscle has a longer descending branch logically correct without being overwhelmingly comthan ascending branch, similar to that obtained with plex for practical applications. The most important single fiber preparations, if active muscle length is property of a muscle is the fact that it produces a force. measured.62 Asymmetry in the active force-length Intact muscle, as found in the body, consists essentially relationship is expected due to rotation of the muscle of two elements: (1) contractile muscle, made up from fibers8’ the aggregation of fleshy muscle fibers and (2) elastic Due to the lack of a muscle model which incorpor“passive” tissue, contributed predominantly by the ates the effects of differing muscle architecture and connective tissues. The passive component is comasymmetry of the length force property of muscles, posed of a series-elastic component and a parallelelastic component. The three-component model is a new mathematical model of the normalized active length-tension relationship was developed.56 A final useful in describing the properties of mammalian mathematical form for the length-tension relationmuscles. ship to be applied in muscle modeling is as follows: Actiue elements. When skeletal muscle is stimulated, it is rapidly activated and changes from a passive tissue (E + 1)P.%3”1 -A)1 - 1.0 * into a dynamic tissue capable of developing a force. P,(t, i,) = exp 0.3534 1 - i,) The contractile element models the active part of the muscle. The dynamics of this element are defined by for i, < 1; the length-tension and force-velocity relationships of muscle contractile tissue. i,) = exp[ - (2.727*ln(E + 1))2] for i, = 1, (7) Length-tension properties of contractile material. where P, is the normalized active muscle tension Blix,” using isolated frog muscles, first demonstrated (F/F_), i, is the muscle architecture, E is the muscle the relationship between the resting length of a muscle strain (I- 1,/l,), 1 is the instantaneous muscle length, and the tension it could generate when contracting. and 1, is the muscle optimum length. This muscle This data forms the Blix curve which demonstrates model predicts the normalized active muscle force, @,, the relationship of total muscle force, passive stretch as a function of the muscle strain, t, and the index of force and muscle contractile force to the length of the architecture, i,. The model was verified by comparing muscle. it to published data on muscle length-force relationIt has been demonstrated that in an isolated frog ships. The isometric twitch force at different lengths muscle fiber@ or intact skeletal muscle’2 the active

The absolute strength of human skeletal muscle has been the subject of interest for three-quarters of a century. In 1910, Fick34 reported that human skeletal muscles exert an absolute stress of 100 N/cm*. Since then, diverse opinions of maximal muscle effort have been reported in the literature. These estimates”4’1).**,25,39,43.47.48,60,73,87 have ranged from 10 to

u

I>

K. R. KAUFMAN et (Il.

786

MEDIAL

SEMIMEMBRANOSUS

GASTROCNEMIUS

Ia - 0.703

Ia - o.se*

-0.2 (4

-0.1 MUSCLE

0;o

-0.2

0.1

-0.1 MUSCLE

STRAIN

0.0

0.1

STRAIN

(W

Fig. 3. Normalized active length-tension data of the medial gastrocnemius (A) and semimembranosus (B) based on the data of Woittiez,r5 compared to the length-tension relationship predicted by the muscle model given in the text.

was measured in situ for the m. gastrocnemius, caput mediale (GM) and the m. semimembranosus (SM) of 13 male Wistar ratss5 The raw data and predicted relationship [equation (711are shown in Fig. 3. A linear regression was performed on the observed and predicted data. The predicted active forcelength relations of these muscles fit well (r= 0.950for GM; r = 0.951 for SM) with the experimental data. Force-velocity properties of contractile material. There have been many attempts to mathematically relate the velocity of shortening of an isotonically contracting skeletal muscle to the force it develops (isotonic force-velocity curve). It has been shown that as the speed of shortening increases, the muscle force decreases. Perhaps the best known relationship between force and velocity of muscular contraction is based on Hill’s equation4 which shows that a hyperbolic relation exists between muscle tensile force and velocity. The higher the load, the slower the contraction velocity. The higher the velocity, the lower the tension. This is in direct contrast to the viscoelastic behavior of a passive material, for which higher velocity of deformation calls for higher forces that cause the deformation. Therefore, the active contraction of a muscle has no resemblance to the viscoelasticity of a passive material. If a load greater than the maximum isometric force, F, is applied to a muscle, the muscle will lengthen at a velocity related to the load, F. The surprise turns out to be that the steady speed of lengthening is much smaller than would be expected from an extrapolation of Hill’s equation to the negative velocity region. In fact, - dF/dv, the negative slope of the force-velocity curve, is about six times greater for slow lengthening than for slow shortening.54 The velocity of lengthening increases with load until the normalized force ratio, F/F,, reaches a critical value at which the muscle “gives”, or increases length rapidly, when the load is raised above the critical value. Values suggested by the literature range from no change45 to increases of 20 to ~oyo~35.51.54.57.63

Several mathematical relationships have been presented to account for the whole range of negative velocities (shortening or concentric contractions) as well as positive velocities (lengthening or eccentric contractions). Audu and Davy6 and Otte# presented discontinuous hyperbolic relationships. Hatze42 published a continuous relationship as follows: P”(rj) = 0.1433[0.1074 + exp - 1.409 sinh(3.2 rj + 1.6)1-l,

(8)

where P, is the normalized active muscle tension due to force-velocity relationship (F/F,) and rj is not-malized contractile element velocity. This relationship is shown in Fig. 4. The normalized contractile element velocity may be expressed as: rj=&, 6,X

(9)

where i is the muscle strain rate (muscle lengths/s) and i,,, is the maximum muscle strain rate of the contractile element (muscle lengths/s). The muscle strain rate, 1, may be calculated from the following formula: i = p”;;? I?

(10)

where ?O,iis the relative velocity of the muscle origin with respect to the muscle insertion, and ii is the unit muscle vector of ith muscle. Both pO,iand i, need to be in the reference frame of the muscle insertion coordinates. There are conflicting reports in the literature regarding the maximum strain rate of the contractile element, i,,. One report6* states that the maximum velocity is obtained by multiplying the mean fiber length by a constant value of 2.5 s-l, while another reportM states that the maximum velocity is 2.5 times the rest length of the muscle. White and Winters0 state that the maximum velocity for the gastrocnemius contracting isotonically under zero load conditions is 10 fiber lengths/s. They further propose that i,,, for

Phy~io~o~cal prediction of muscle forces--I ,.4

2 0

787

~CONCE~TRIC-~_ECCENTR~C-

1.2-

5

-1.0

-0.0

-04

-0.4

NORMALIZED

Fig. 4. Nomadic

muscle force-vekxity

-0.2

MUSCLE

0.0

0.2

0.4

,

VELOCITY

relationship for concentric and eccentric ~ntractions

each muscle modeled should then be derived from their percentage fast twitch content relative to the gastrocnemius. The scaling of the maximal contraction velocity is appropriate, since it has been demons~t~ that the maximum shortening velocity of muscle is a distinct function of the muscle fiber composition. The force-velocity relationship was examined in the human knee extensors and muscle biopsies were obtained from the vastus lateralis of 25 male subjects.” It was demonst~ted that the individual maximal speeds of contraction were si~ifi~ntly related to the percentage of fast twitch muscle fibers.% Other authors have demonstrated that for various muscles, the maximal contraction velocity increases with increased fast twitch muscle fiber percentage.6gJ6*79*s2 These results support the general correlation, although it need not be a cause and effect ~lationship, between a muscle’s intrinsic speed and its myosin ATPase activity. When all of these data are put together, and weighted based on merit, the following material property approximations emerge: slow--&,, = 2.0 muscle lengths/s; fast---Z,, = 7.5 muscle lengths/s. Thus, i, for each muscle modeled will then be derived from its muscle fiber type distribution. Woittiez et ~1.‘~studied muscles of differing architecture but equal optimum length. They found that normalizing the force-velocity relationships of muscles with respect to maximal active force and maximal velocity of shortening deletes all e&cts of architecture. Therefore, equation (8) does not need to be modified for differing muscle architecture. Passive elements. During either an isotonic or an isometric contraction, a portion of the total work expended is utilized to deform or to stretch the non~nt~ctile elements of the muscle. Thus, only part of the energy liberated during a contraction is employed to move a load or to develop tension in a muscle. The contractile elements also must act against the viscosity of the sarcoplasm and the connective tissue elements, including the tendons. When an active

(Hatze4*).

human skeletal muscle is forcibly stretched, the muscle stores elastic energy, part of which can be recovered during the subsequent active shortening of the entire muscle.i9 Series elastic component. Evidence for the existence of a lightly damped elastic component in series with the contractile component is based primarily on the biphasic response of a muscle when it is released during isometric contraction and allowed to shorten freely against a load that is less than the isometric tension at the time of release. The initial rapid phase of shortening is completed within a few milliseconds and has been attributed to rapid release of energy stored in the elastic component after the change in load. The second slow phase is the result of shortening of the contractile component in accordance with its force-velocity properties. In the body, most of this elasticity is located in the tendon and in the actin and myosin cross-bridges of the muscle.” In attempting to evaluate the properties of the series elastic component that have been determined from measurements on whole muscle it is important to remember that the proportion of tendon in the total length of the muscle varies greatly from one muscle to another. The series elastic element (SE) is comparable to a nonlinear spring and has been modeled as a polynomial” or with an exponential relationship.“ Parallel elastic component. The passive elastic elements are due to the inguence of st~ctures which lie in parallel with the force producing sarcomeres. Passive muscle held at lengths greater than the rest length develop force. The tension exerted by the parallel elastic element is independent of the contractile element and represents tension in the connective tissue of the muscle. This element therefore provides a passive component of tension which must be added to the active tension of the contractile element to obtain the total tension exerted by the muscle unit. The parallel elastic element can be modeled with an exponential force-strain function, which also depends

78X

K. R.

KAUFMAN et al.

MUSCLE STRAIN

Fig. 5. Normalized passive length-tension relationship for muscles of differing geometry. The index of architecture, i,, describes the geometry and is defined as the muscle fiber length (&) divided by the muscle belly length (L,). on the muscle architecture. 84,85 The passive forcestrain curve for differing muscle architecture is shown in Fig. 5. Modeling

the musculature:

physiological

factors.

With this understanding of muscle physiology, a mathematical model of muscular contraction can be developed. Active components. The contractile element is the active component in the model. The tension developed by the contractile element (Fm) is a function of three variables: the neuromuscular activation (a), the length-tension relationship @,), and the force velocity relationship @J. The word “activation” is used in this paper to describe the number of active motor units and their frequencies of firing in response to stimulation through motor nerve fibers. The activation, a, varies between 0 and 1. A separation of the activation, length-tension relationship, and forcevelocity relationship, respectively, is well founded.9,‘o Thus, P,, = a .P,.@,.

(11)

The F, curve is assumed to be the same among actuators. Passive components. Passive muscle forces develop at lengths greater than the rest length of the muscle. It has been assumed that this normalized force arising from elastic structures internal to the muscle fiber depends on normalized muscle length and that this relationship is the same among actuators.” Obviously, in any one muscle and at any one time, the amount of connective tissue is constant and consequently the passive tension does not change with activation level. A potentially significant element in any muscle modeling attempt is the SE. The question that needs to be answered is-how signScant is the tendon compliance of the leg muscles to the analysis? This finding will have considerable impact on the solution method used for the force optimization. Certainly, the compliance of the tendons of the horse or kangaroo is

high enough to make them significant energy storage elements, but is the same true of human tendons? The calcaneal or Achilles tendon provides a worst case example because of its tendon length to muscle length ratio. The medial gastrocnemius has the lowest index of architecture (i. = 0.142) which means that it will have the largest change in force resulting from a change in muscle strain. The strain of a human tendon during normal use is approximately 3%.5 Substituting these conditions into the force-length relationship [equation (7)], the change in force is -0.09 6 AF/F, i 0.14. Likewise, the velocity difference between the ends of the tendon is found from the time derivative of the elongation. Expanding the analysis above, the ‘us muscle change in velocity for the medial gaatra is 0.002 < v_/v, < 0.007. These figures indicate that tendon elasticity may be a significant energy storage element; however, this is a worst case analysis and would apply to only a few of the muscles considered in this study. The majority of the leg muscles have much shorter or less slender tendons than the Achilles tendon, so this small subset of significant elastic elements does not warrant the difficulties of the requisite optimization procedures to include the series elastic element. Muscle force constraint. The total force developed by the simpli&d muscle model can now be expressed as the sum of active and passive tensions: B=&+&,

(12)

where F is the total normalized force in the muscle. Because muscles only exert force in tension and not in comprehension, a proper muscle force solution requires imposition of a muscle force constraint. This constraint requires that any individual mu&e must have a force greater than or equal to zero and less than an upper limit; thus, O
i=

l,m.

(13)

Physiological prediction of muscle force-1 The upper bounded solution space for the muscle force, Fy, is determined by the muscle stress limit, u, the PCSA, the normalized muscle active force correction, &, for the length-tension and force-velocity relationships, the normalized muscle passive force characteristics, Pp, and the muscle activation, a. This constraint is considered as a means of incorporating physiological information into the musculoskeletal model.

789

The approach to the optimization problem taken in this paper is to minimize the amount of neuromuscular activation in order to obtain a unique solution. The underlying assumption is that the load sharing between the muscles is unique during learned motor activities, and that the neurological control of muscle action is governed by certain physiological criteria that guarantee efficient muscle actions. The objective function should correspond to these physiological criteria. Intuitively, it seems reasonable that the control process of human motion involves the optimization of Optimal criteria body function. This general notion was first suggested The unique solution for distributing muscle forces by the Weber brothers (1836), who stated that man is then obtained by minimizing the neuromuscular walks in “the way that affords us the highest energy activation, a. Thus, the objective function is expenditure for the longest time and with the best Minimize LX. (14) results”.27 An optimization criterion based on the minimization of muscle energy consumption,4’ using In this approach, constraint equations for dynamic a thermodynamic relationship between input chemical equilibrium [equation (l)] and muscle force [equation power and output mechanical power, failed to result in (13)], as well as the objective function [equation (14)] muscle force predictions which differed significantly consist of linear combinations of the unknown varifrom other optimal criteria. However, promising ables. Thus, the well-known linear programming results have been obtained when the endurance time algorithm’* can be used for solution, and the solution of the activity is maximized and, hence, muscular is guaranteed to converge on a global minimum. fatigue is minimized. 30*31 The present work adopts a similar hypothesis that the neuromuscular system minimizes the amount of neuromuscular activation to DISCUSSION achieve muscular contractions. An important aspect Given a particular activity to be performed with a of the use of this criterion is that the objective function large number of potential force-carrying structures has a physiological basis. Thus, the unique feature at a joint, an individual may exercise considerable in this paper is the incorporation of physiological discretion in generating force in the actively conphenomena into the optimal criterion and muscle trolled muscles. The distribution problem at a joint is, force constraints and the suggestion of how the human typically, an indeterminate problem, since the numbody activates individual muscles during activity. ber of muscles and articular contact regions available The length-tension relationship of pennifomr to transmit force across a joint, in many cases, muscle could be accounted for by considering the exceeds the minimum number of equations required trigonometric relationship between the long axis of to generate a determinate solution. This condition the muscle fibers and the long axis of the muscle belly. results in an underdetermined set of equations. An In this situation, the contributions of the contractile infinite number of possible solutions exist for the elements are transmitted along the fiber centroids and indeterminate equations. One method of solution the fibers can be assumed to be symmetrically distriwithout simplification is that of seeking an optimum buted about the common tendon of the penniform muscle. The total muscle force can then be taken to solution. Previously, the optimization of muscular force, be the summation of forces from the individual fibers, power, or energy tended to assign all of the activation effectively resulting in a cancellation of forces transto a single optimal muscle. This would lead to fatigue verse to the long axis of the muscle. However, pennate in that muscle which would not be accounted for by fibers tend to swing about their origin during contracthe optimization method. In addition, electromyotion. The path traveled by the insertion is generally greater than the distance the fibers shorten. The graphic evidence shows that muscles tend to act synergistically as a group.66 Hence, if there is no actual travel of the tendon will have to be calculated known physiological evidence incorporated into the from the fiber length, the fiber contraction coefficient, the initial angle of pennation, and the change in optimization method, the results are not physiologically viable. Patriarco et af.@ concluded that distance between the parallel attachment surfaces of the fibers.‘4,62 Also, the effective force component “the imposition of additional physiologically-based contributed by a fiber to the central tendon will constraints [is] essential to distinguishing the role decrease with increasing angle of pennation. Finally, of individual muscles”. Several investigators have attempted to combine mathematical formulations it is necessry to average or integrate the forces exerted with consideration of physiological function.29.4’,42,59,70 by the individual fibers during a contraction and combine them into a consolidated length-tension relationWhile their approaches have varied considerably, ship. Consequently, theoretical considerations of the their efforts represent an improvement over previous phenomenon of pennation will lead to complexities modeling efforts.

K. R.

7YO

1(AUFMAN

that may be much greater than originally anticipated. It is for these reasons that a mathematical relationship of the normalized active length-tension was used which consolidated the effects of differing muscle

architecture and as~me~y of the length force property of muscles. For simplicity, the SE was considered to be a stiff spring. All force was then transmitted directly to other elements of the muscle model. This contention is supported by Bawa et al.” who found the stiffness of the SE to be much larger than the parallel elastic element. However, in the future, a more detailed analysis of muscle function may want to examine the tendon contribution more closely. The word “activation” has been used in this paper to describe both the number of active units and the degree of their activity, since the mechanical response is determined by both of these factors. The available evidence suggests that the patterns of motor unit activity are well suited to the demands place on the motor system. Joyce et aI. showed that with different activation levels the basic shape of the force curve remains approximately the same. Phillips and Petrofsky69 documented the effect of activation level on the force-velocity level, showing that both peak load and peak velocity were modified. Hence, less than fully activated muscle, i.e. tl < 1, was assumed to have a contractile element relationship that is identical in shape to fully activated muscle but reduced in ampIitude. At the neurophysiolo~cal level, future refinements to the model could include basic muscle phenomena such as (a) recruitment patterns and (b) decay of force with fatigue. The contractile element is the only active component in the model. Its behavior was assumed to depend nonlinearly on its length, velocity of movement and degree of activation. Furthermore, it is known that most mammalian muscles are composed of different fiber types based on their contractile and histochemical

1’1 tri.

properties. It is thought that during voluntary muscle contractions, muscle units are recruited in order of

diminishing fatigue resistance. The total force output of a muscle would then be the sum of the force output from each individual muscle fiber type. In any case, this model provides a starting point for theoretical investigations of muscle force predictions which may lead to the use of better models with more accurate results in the future. Several elements could be incorporated to improve the model’s predictive ability. It would be naive to suggest that the model developed in this paper accurately describes the complex processes for the recruitment of the musculoskeletal system. However, it has been an attempt to incorporate physiological data in the solution of muscle force distribution. As suggested, individual muscle forces can be predicted using optimization techniques. The optimi~tion approach, in contrast with the use of greatly simplified, mechanical system models that permit easier solutions, maintains the model correspondence with the real physical system. The virtue of this approach allows the establishment of future study in clearly identifiable directions.

CONCLUSION

In summary, it is important to select muscle prediction criteria based on sound physiological principles. Accordingly, known physiologicat principles have been incorporated into a mathematical model for predicting muscle forces. A new optimization approach of minimizing muscle activation is suggested. The methodology proposed in this paper represents a rational method for dealing with the complex redundant problems encountered while attempting to solve mu~uloskeletal systems biomecha~caIIy. Acknowledgements-Supported in part by NIH Grants CA 23751, AR 26287, and the Bush Foundation.

1. Alexander R. McN. and Vernon A. (1975) The dimensions of knee and ankle muscles and the forces they exert. J. Hum. Movement Studies 1, 115-123. 2. An K. N., Hui F. C., Morrey B. F., Linscheid R. L. and Chao E. Y. S. (1981) Muscles across the elbow joint: a biomechanical analysis. J. Biomech. 14, 659-669. 3. An K. N., Kwak B. M., Chao E. Y. and Morrey B. F. (1984) Determination of muscle and joint forces: a new technique to solve the inde~~na~ problem. Tram ASME 106, 364-367. 4. An K. N., Kaufman K. R. and Chao E. Y. S. (1989) Physioio~~l ~nsid~ations of muscle force through the elbow joint. J. Biomech. 22, 1249-1256. 5. An K. N. (1989) Internal publication, Biomechanics Laboratory, Mayo Clinic, Rochester, MN. 6. Audu M. L. and Davy D. T. (1985) The influence of muscle model complexity in musculoskeletal motion modeling. J. Biomech. Engng 107, 147-156. 7. Bahill A. T. (1981) Bioengineering: Biomedical, Medical and Clinical Engineering. Prentice-Hall, En&wood Cliffs, NJ. 8. Rahier A. S. (1968) Modeling of mammalian skeletal muscle. IEEE Tram Biomed. Engng BME lS, 249-257. 9. Bahler A. S., Fales J. T. and Zierler K. L. (1%7) The active state of rn~rn~~an skeletal muscle. J. gen. Pllrysiol. 50, 2239-2253. 10. Bahier A. S., Fales J. T. and Zierler K. L. (1968) The dynamic properties of mammalian skeletal mu&e. J. gen. Physiol. 51, 369-384. 11. Bahler A. T. (1967) Series elastic component of mammalian skeletal muscle. Am. J. Physiol. 212, 1560-1564. 12. Banus M. G. and Z&in A. M. (1938) The relation of isometric tension to length in skeletal muscle. J. cell. camp. Physiol. 12. 403-420.

Physiological prediction of muscle forces-1

791

13. Bawa P., Mannard A. and Stein R. B. (1976) Predictions and experimental tests of a v&o-elastic muscle model using elastic and inertial loads. Biol. Cyberner. 22, 139-145. 14. Benninghoff A. and Rollhatlser H. (1952) Zur inneren mechanik des gefiederten Muskels. PJliigers Arch. ges. Physiol. 254, 527-548.

15. Blix M. (1894) Die Lange und die Spannung des Muskels. &and. Arch. Physiol. 5, 149-206. 16. Boichut D. A. and Valentini Fr. A. (1983) Human locomotion analvsis. Determination of muscular forces and nervous orders. Inr. J. Bio-med. Comp. 14, 217-230. 17. Brand P. W., Beach R. B. and Thompson D. E. (1981) Relative tension and potential excursion of muscles in the forearm and hand. J. Hand Surg. 6, 209219. 18. Brand R. A., Crowninshield R. D., Wittstock C. E., Pedersen D. R. and Clark C. R. (1982) A model of lower extremity muscular anatomy. J. Biomech. Engng 104, 304-310. 19. Cavagana C. A., Dusman B. and Margaria R. (1968) Positive work done by a previously stretched muscle. J. appl. Physiol. 24, 21-32. 20. Chao E. Y. S. (1971) Determination

of applied forces in linkage systems with known displacements: with special application to biomechanics. Ph.D. dissertation, University of Iowa, Iowa City. 21. Chao E. Y. and An K. N. (1978) Determination of internal forces in the human hand. J. Engng Mech. Div. Am. Sot. Civ. Engrs 104, 255-272. 22. Chao E. Y. and An K. N. (1978) Graphical intepretation of the solution to the redundant problem in biomechanics. J. Biomech. Bngng 100, 159-167. 23. Chow C. K. and Jacobson D. H. (1971) Studies in human movement locomotion via optimal programming. Math. Biosci. 10, 239-306. 24. Close R. I. (1972) Dynamic properties of mammalian skeletal muscles. Physiol. Rev. 52, 129-197. 25. Crowninshield R. D. (1978) Use of optimization techniques to predict muscle forces. J. Biomech. Engng 100, 88-92. 26. Crowninshield R. D. and Brand R. A. (1981) A physiologically based criterion of muscle force prediction in locomotion. J. Biomech. 14, 793-801. 27. Crowninshield R. D. and Brand R. A. (1981) The prediction of forces in joint structures: distribution of intersegmental resultants. Exercise Sports Sci. Rev. 9, 159-181. 28. Dantzig G. B. (1963) Linear Programming and Extensions. Princeton University Press, Princeton.

29. Davy D. T. and Audu M. L. (1987) A dynamic optimization technique for predicting muscle forces in the swing phase of gait. J. Biomech. 20, 187-201. 30. Dul J. (1983) Development of a minimum-fatigue optimization technique for predicting individual muscle forces during human posture and movement with application to the ankle musculature during standing and walking. Ph.D. dissertation, Vanderbilt University, Nashville, TN. 31. Dul .I., Johnson G. E., Shiavi R. and Townsend M. A. (1984) Muscular synergism II. A minimum-fatigue criterion for load sharing between synergistic muscles. J. Biomech. 17, 675-684. 32. Edgerton V. R., Smith J. L. and Simpson D. R. (1975) Muscle fibre type populations of human leg muscles. Histochem. J. 7, 259-266. 33. Enderle J. D. and Wolfe J. W. (1987) Time-optimal control of saccadic eye movement. IEEE Trans Biomed. Engng 6, 113-120. 34. Fick R. (1910) Hanabuch der Anatomie des Menschen, Vol. 2, Verlag von Gustav Fischer, Jena.

35. Flitney F. W. and Hirst D. G. (1978) Cross-bridge detachment and sarcomere ‘give’ during stretch of acting frog’s muscle. J. Physiol. 276, 449-465. 36. Freivalds A. (1985) Incorporation of active elements into the articulated total body model. Publication No. AAMRL-TR-85-061, U.S. Air Force/AFSC Aeronautical Systems Division, Wright-Patterson AFB, Ohio. 37. Froese E. A. and Houston M. E. (1985) Torque-velocity characteristics and muscle fiber type in human vastus lateralis. J. appl. Physiol. 59, 309314. 38. Fung Y. C. (1981) Biomechanics:

Mechanical Properties of Living Tissues, pp. 302-327. Springer, New York. 39. Gans C. and Bock W. J. (1965) The functional significance of muscle architecture-a theoretical analysis. In Ergebnisse der Anatomie und Enrwicklungsgeschichte, pp. 115-142, Springer, Berlin. 40. Gordon A. M., Huxley A. F. and Julian F. J. (1966) The variation in isometric tension with sarcomere length in vertebrate muscle fibres. J. Physiol. 184, 170-192. 41. Hardt D. E. (1978) A minimum energy solution for muscle force control during walking. Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, MA. 42. Hatze H. (1981) Myocybernetic Conrrol Models of Skeletal Muscle. University of South Africa, Pretoria. 43. Haxton H. A. (1944) Absolute muscle force in ankle flexors of man. J. Physiol. 103, 267. 44. Hill A. V. (1938) The heat of shortening and the dynamic constants of muscle. Proc. R. Sot. Lond. 126B, 136-195. 45. Hof A. L. and Van Den Berg J. W. (1981) EMG to force processing II: estimation of parameters of the Hill muscle model for the human triceps surae by means of a calfergometer. J. Biomech. 14, 759-770. 46. Huijing P. A. and Woittiez R. D. (1985) Notes on planimetric and three-dimensional muscle models. Neth. J. Zool. 35, 521-525. 47. Ikai M. and Fukunaga T. (1968) Calculation of muscle strength per unit cross-sectional area of human muscle by means of ultrasonic measurements. Ini. Z. Angewandie Physiol. 26, 26-32. 48. Inman V. T. and Ralston H. J. (1968) Human Limbs and Their Substitutes (eds Kloostee _ I P. E. and Wilson P. D.).,. pp. 296-3 17. McGraw-Hill, New York. 49. Jensen R. H. and Davy D. T. (1975) An investigation of muscle lines of action about the hip: a centroid line approach vs the straight line approach. J. Biomech. 8, 103-110. 50. Johnson M. A., Polgar J., Weightman D. and Appleton D. (1973) Data on the distribution of fibre types in thirty-six human muscles: an&topsy study. J. Neurosci. 18, 11l-129: 51. Jovce G. C. and Rack P. M. H. (1969) Isotonic lengthening . I - and shortening movements of cat soleus muscle. J. Physiol. 20.4, 4755491. 52. Joyce G. C., Rack P. M. H. and Westbury D. R. (1969) The mechanical properties of cat soleus muscle during controlled lengthening and shortening movement. J. Physiol. 284, 461-474. I

K. R.

792

KAUFMAN

et ul.

53. Kardel T. (1986) A specimen of observation upon the muscles. In Nicoluus Steno 1638.-1686 (eds Paulsen J. E. and Snorrason E.), pp. 97-134. Nordisk Insulinlaboratorium, Denmark. 54. Katz B. (1939) The relation between force and speed in muscular contraction. J. Physiol. %, 45-64 55. Kaufman K. R., An K. N. and Chao E. Y. S. (1991) A mathematical model for the dynamics of the knee joint during isokinetic exercise. /. Biomech. Engng (in press). 56. Kaufman K. R., An K. N. and Chao E. Y. S. (1989) Incorporation of muscle architecture into the muscle length-tension relationship. J. Biomech. 22, 943-949. 51. Komi P. V. (1973) Measurement of the force-velocity relationship in human muscle under concentric and eccentric contractions. In Medicine and Sport Vol. 8: Biomechanics III, pp. 224-229. Karger, Basel. 58. Lew W. D. and Lewis J. L. (1977) An anthropometric scaling method with application to the knee joint. J. Biomech. 10, 171-181.

59. Mikosz R. P. (1985) Mathematical model for the study of forces in the human knee joint during locomotion. Ph.D. dissertation, University of Illinois, Chicago, IL. 60. Morris C. B. (1948) The measurements of the strength of muscle relative to the cross-section. Res. Q. 19, 2955303. 61. Morrison J. B. (1967) The forces transmitted by the human knee joint during activity. Ph.D. dissertation, University of Strathclyde, U.K. 62. Muhl 2. F. (1982) Active length--tension relation and the effect of muscle pinnation on fiber lengthening. J. Morphol. 173, 285-292.

63. Otten E. (1987) A myocybernetic mode1 of the jaw system of the rat. Orofacial Research Group, Bloemsingel 10, 9712 KZ Groningen, The Netherlands. 64. 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. Biomech. 14, 513-525. 65. Patriarco A. G. (1982) The prediction of individual muscle forces during human movement. Ph.D. dissertation, Harvard University, Cambridge, MA. 66. Paul J. P. (1972) Comparison of EMG signals from leg muscles with the corresponding force actions calculated from walkpath measurements. In Proceedings on Human Locomozor Engineering. University of Sussex. Inst. of Mech. Engng. London. 67. Paul J. P. (1965) Bioengineering studies of forces transmitted by joints (II). In Engineering Analysis in Biomechanics and Related Bio-engineering Topics, pp. 369-380. Pergamon Press, Oxford. 68. Pedotti A., Krishnan V. V. and Stark L. (1978) Optimization of muscle-force sequencing in human locomotion. Math. Biosci. 38, 57-76.

69. Phillips C. A. and Petrofsky J. S. (1980) Velocity of contraction of skeletal muscle as a function of activation and fiber composition: a mathematical model. J. Biomech. 13, 549-558. 70. Pierrynowski M. R. and Morrison J. B. (1985) A physiological model for the evaluation of muscular forces in human locomotion: theoretical aspects. Math. Biosci. 75, 699101. 71. Rack P. M. H. and Westbury D. R. (1974) The short range stiffness of active mammalian muscle and its effect on mechanical properties. J. Physiol. 240, 331-350. 72. Ramsey R. W. and Street S. F. (1940) The isometric length-tension diagram of isolated skeletal muscle fibers of the dog. J. cell. camp. Physiol. 15, 1l-34. 73. Rechlinghausern N. (1920) Gliedermechanik und Luhmungsprothesen. Springer, Berlin. 14. Schultz A., Haderspeck K., Warwick D. and Portillo D. (1983) Use of lumbar trunk muscles in isometric performance of mechanically complex standing tasks. J. Ortho. Res. 1, 77-91. 75. Seireg A. and Arvikar R. J. (1973) A mathematical model for evaluation of forces in lower extremities of the musculoskeletal system. J. Biomech. 6, 3 133326. 76. Spector S. A., Gardiner P. F., Zernicke R. F., Roy R. R. and Edgerton V. R. (1980) Muscle architecture and force-velocity characteristics of cat soleus and medial gastrocnemius: implications for motor control. J. Neurophysiol. 44, 951-960.

77. Stephenson D. G. and Williams D. A. (1982) Effects of sarcomere length on the forcepCa relation in fast- and slow-twitch skinned muscle fibres from the rat. J. PhysioL 333, 637-653. 78. Thorstensson A., Grimby G. and Karlson J. (1976) Force-velocity relations and fiber composition in human knee extensor muscles. J. appl. Physiol. 40, 12-16. 79. Wells J. B. (1965) Comparison of mechanical properties between slow and fast mammalian muscles. J. Physiol. 178, 252.-269.

80. White S. C. and Winter D. A. (1986) The prediction of muscle force using EMG and a muscle model. Proceedings of the North American Congress on Biomechanics, Montreal (Quebec), Canada, 25-27 August, pp. 67-68. 81. Wickiewicz T. L., Roy R. R., Powell P. L. and Edgerton V. R. (1983) Muscle architecture of the human lower limb. Clin. Ortho. Rel. Res. 179, 275-283.

82. Winters J. M. (1985) Generalized analysis and design of antagonistic muscle models, effect of nonlinear muscle-joint properties on the control of fundamental movements. Ph.D. dissertation, University of California, Berkeley, CA. 83. Winters J. M. and Stark L. (1985) Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models. IEEE Tram Biomed. Engng BME 32, 826-839. 84. Woittiez R. D., Huijing P. A. and Rozendal R. H. (1983) Influence of muscle architecture on the length-force diagram: a model and its verification. Pjlcgers Arch. ges. physiol. 397, 73-74. 85. Woittiez R. D., Huijing P. A., Boom H. B. K. and Rozendal R. H. (1984) A three-dimensional muscle model: a quantified relation between form and function of skeletal muscles. J. Morphol. 182, 95-l 13. 86. Woltring H. J. (1986) A Fortran package for generalized, cross-validatory spline smoothing and differentiation. In Advances in Engineering Software. C.M.L. Publications, Southampton, U.K. 87. Yamada H. (1970) Sfrength of Biological Maferials, pp. 93-97. Williams & Wilkins, Baltimore, MD. 88 Zajac F. E., Topp E. L. and Stevenson P. J. (1986) Musculotendon actuator models for use in computer studies and design of neuromuscular stimulation systems. Proceedings of %h Annual Conference Rehabilitation Engineering Society, Minneapolis, MN, 23327 June. (Accepted 19 July 1990)