A mathematical model for predicting relative muscle force with a perturbation analysis of selected muscle parameters

J. theor. Biol. (1989) 138, 213-233

A Mathematical Model for Predicting Relative Muscle Force with a Perturbation Analysis of Selected Muscle Parameters B. A. YOUNG? AND M. A. ASHLEY~

? Department of Biological Sciences, University of Calgary, Calgary, Alberta T 2 N 1N4, Canada and ~. Department of Biological Sciences, NA U Box 5640, Northern Arizona University, Flagstaff, Arizona 86011, U.S.A. (Received 29 April 1988, and accepted in revised form 24 January 1989) A mathematical model is presented that predicts relative muscle forces using a minimum of experimentally derived input data. Tests of this model against literature values for maximum muscle force of four cat hindlimb muscles show a maximum error of only 5%. A perturbation analysis using this model demonstrates its sensitivity and applicability, as well as the congruence between this model and previous theoretical discussions of muscle function. Introduction

Throughout the long history of morphology, researchers have focused a disproportionate share of their attention upon the musculo-skeletal system. The relative ease of examination of the musculo-skeletal system, as compared to other anatomical and physiological systems, and the applicability of mechanical analysis as a tool to understand function were major factors in maintaining this emphasis. Even with the proliferation of modern physiological techniques to examine muscle function (e.g. histochemistry, electromyography, etc.), the mechanical analysis o f force lever systems has remained an important tool in functional morphology. Biomechanical analyses o f muscle-bone force lever systems enable researchers to explore many aspects of organismal function including the biological roles, functional integration, and phylogenetic specialization o f the musculo-skeletal system. In order to extrapolate muscle function from one organism to a closely related organism, or to examine functional consequences of the structural organization of muscle (e.g. pinnation angle, cross sectional area, etc.), it is necessary to predict muscle force values. Although an ideal model would estimate the actual force produced by a muscle, this is an almost impossible task due to the range o f variation between individuals and within the physiological state of an individual. Relative muscle force values, with few exceptions (e.g. the estimation of absolute joint excursion), provide as much functional information as absolute muscle force values. Typically, models attempting to predict relative muscle force center on one of three approaches: electromyography, physiology, or biomechanics. While some ? To whom all correspondence should be sent. 1:Current address: Department of Cell and Developmental Biology, University of California, lrvine, CA 92717, U.S.A. 213 0022-5193/89/100213+21 $03.00/0 © 1989 Academic Press Limited

214

B.A. YOUNG AND M. A. ASHLEY

models incorporate aspects of all three, this division is helpful in understanding the different approaches to the problem of relative muscle force. Models based on electromyography (for recent examples see Weijs & Dantuma, 1975; Hof& Van den Berg, 1977; Kardong et al., 1986) all assume that muscle force can be correlated to the intensity and amplitude of EMG signals. This correlation is not widely accepted; Loeb & Gans (1986) state that: "Variously processed EMG signals are frequently offered as a direct indication of the force output of the muscle, which is most unlikely at best". Criticisms leveled at these models center on their failure to incorporate the muscle's length-tension curve, biomechanical properties, and the inherent variation of EMG signals (including potential variations and errors produced during integration). Notwithstanding these criticisms, many studies have relied upon EMG data to predict muscle force. The second group of models rely on physiological and anatomical properties of muscle to estimate relative muscle force. The most accurate model of this type is the work of Poweli et al. (1984) on the hindlimb of the guinea pig. They estimated maximum tetanic tension by multiplying the cross-sectional area by the specific tension of the muscle. The formula used by Powell et al. (1984) for cross-sectional area incorporated muscle mass, angle of pinnation, fiber length and muscle density. Pinnation angle is an important functional factor (Gans & Bock, 1965; Gans, 1982) that can be directly measured (Sacks & Roy, 1982). Fiber length and muscle density can be determined experimentally (Mendez & Keys, 1960; Murphy & Beardsley, 1974; Sacks & Roy, 1982). The specific tension used by Powell et al. (1984) was taken from the literature (Spector et al., 1980; Roy et al., 1982). The model of Powell et al. (1984) appears to be quite accurate in that the predicted muscle forces exhibit a strong linear relation to the measured muscle forces. While this model will be o f benefit within this experimental system, at this time the degree to which it can be applied beyond the guinea pig remains to be seen. The major problem in applying this model to a biomechanical analysis (note that this was not the goal of Powell and his associates) is that it produces only maximum force, not a more complete picture of dynamic changes in force output. Additionally, the usefulness of this model for prediction is limited by its reliance upon experimentally derived input data. The biomechanical models, the third major group, rely predominantly upon the mechanical location of the muscle to determine its relative force. Cross-sectional area, assumed to be directly proportional to force, is estimated by dividing muscle mass by muscle length. Although this simple formula forms the basis for most biomechanical models, it can be improved upon by correcting for fiber angle (Gans & Bock, 1965; Gans, 1982; Gans et al., 1985). The action of a muscle upon a lever system is dependent upon the force it produces, and its moment arm. The latter is defined as the distance from the perpendicular point of application of the muscle to the axis of movement of the lever. This value can be obtained by direct measurements (McClearn, 1985) and must be included within any valid biomechanical model. To be realistic, the model must vector the muscle force into three dimensions (as was done by Hiiemae, 1971; Riley, 1985; Sinclair & Alexander, 1987; Weijs et al., 1987). Although few biomechanical models go beyond it, this is not a valid

RELATIVE

MUSCLE

FORCE

PREDICTION

215

stopping point; muscles produce movement by exerting a torque upon a bony lever system. There is a tendency to reduce the dynamics of the muscle-joint system to a single plane (i.e. M~ller, 1987) and to grossly simplify the muscle attachments; these temptations must be resisted as they remove all biological reality from the analysis. The purpose of this paper is to describe a biomechanical model of muscle function, and to explore the implications of this model. This model is designed as a simple, yet hopefully realistic, tool for the prediction of relative forces produced by any musculo-skeletal system. In an effort to make this model accessible to any researcher, we have designed it so as to require a minimum of experimentally derived input data. It must be noted that since this model predicts relative force it can only be used when comparing two muscles. Although the required inputs have been kept to a minimum, this model offers many advantages over previous techniques. Specifically, our model has been designed to work in a dynamic, rather than static, framework; to look at the functional system with a multidimensional perspective; and to carry the analysis through to the calculations of torque and torsion. This mathematical model is amenable to computer modelling, enabling it to be used not only to predict relative forces for existing muscle systems, but also to perform perturbation analyses of force-lever systems. To use this model, the following data must be available: muscle mass, range of joint movement in all three planes, lever arm distances in all three planes, pinnation angle (where applicable), and a basic knowledge of the physiology of the system. This model incorporates the length-tension curve of the muscle as well as some aspects of the intensity and timing of muscle stimulation; although these parameters would ideally be directly measured, values from the literature may be used without appreciable sacrifice of precision. The remainder of this paper will be devoted to a step-by-step description of this model, the testing of its accuracy using literature sources, and providing a preliminary example of this model's capacity to perform perturbation analyses. We specify underlying assumptions in the relevant sections of the step-by-step description of the mathematics involved. In regards to the mathematics, we agree with Otten (1983) who stated that: " . . . the mathematics is usually quite simple and straightforward and that it is laid out in the sub-routines of the model. The complexity of the model is produced by the combination of the sub-routines". Materials and Methods RELATIVE

MUSCLE

FORCE

MODEL

(1) The points of attachment The initial step of this model involves the depiction of the muscle origin and insertion as single points (this has become a standard procedure, see for example Iordansky, 1970; Riley, 1985). We believe that this single point can be satisfactorily determined by simply taking the midpoints of several measurements all passing

216

B. A . Y O U N G

AND

M.

A. ASHLEY

through (or near) the center of the muscle's attachment. There is no doubt that treating muscle attachment as a single point is a simplification; however, we do not feel that this greatly affects the accuracy of the model. (2) Measurement of muscle mass The accurate measurement of muscle mass is essential for this model. We advocate the usage of dry muscle mass (where the muscle is first cleaned of connective tissue, then oven dried) since it is far more reliable than is wet mass. While formulae exist to convert dry muscle mass to "physiological" mass (Jenkins & Weijs, 1979) in a model predicting only relative forces this step serves no function. The input value for muscle mass will ideally be the mean of a set of measurements. (3) Measurement of lever arms We define the lever arm as the distance from the point of attachment o f the muscle to the center of the joint. This differs from the moment arm in that the muscle has not been solved to a perpendicular. Lever arm distance does not change with joint movement or muscle activity. Input values for the lever arms of both the origin and insertion must be measured in all three directions (cranio-caudal, dorso-ventral, and medio-lateral). (4) Measurement of tendon length Where applicable, the length of all tendinous elements in series with the muscle must be measured. Failure to incorporate tendon length greatly inflates the estimated length of the muscle. We assume that the length o f the tendinous elements is constant, even though a small amount of stretch within tendons has been documented (Wainwright et al., 1982). (5) Correction of lever arms In the "typical" bony lever system two of the lever arms will be orientated at 90 degrees to the main axis o f the bony segment, with the third lever arm laying along the main axis of the bony segment. If this situation holds, the only time correction of the lever arm is mandated is when the origin and insertion points are on opposite sides of the joint (as discussed below). In a two joint system, the lever arms about the origin may not be orientated at 90 degrees to the intervening bony segment (due to the movement o f the first joint). As Fig. 1 illustrates, this difference can be corrected for as follows: Corrected Lever Arm = cos lABS ( 9 0 - Angle)] x Lever Arm. In this way, the lever arms in all three directions can be corrected. (6) Correction of joint angles Joint angles are assumed to be measured between the long axes o f the bony segments forming the joint. Since the lever arms are orientated at 90 degrees to

RELATIVE

MUSCLE

FORCE

217

PREDICTION

J /J a

B

B

D

b ("

1

2

3

C/C-.. - ....

4

FIG. 1. Illustration of some of the procedures used to calculate muscle length. A = origin lever a n n ; B = intervening segment; C = insertion lever arm; D = corrected origin lever a n n ; a = origin joint angle; b = insertion joint angle; c = corrected insertion joint angle. The measured anatomical system is depicted in (1). Since the origin lever a n n is not orientated at 90 degrees to the intervening segment, it must be corrected. The corrected lever arm, calculated as cos [ABS ( 9 0 - a ) ] A , is represented by D in (2). Since the intervening segment does not lie in the direction of these two lever arms, it is of no importance in determining muscle length in this direction. The omission o f the intervening s e g m e n t lets us conceptually shift the corrected lever a n n ( D ) to the insertion end o f the intervening s e g m e n t as represented in (3). The insertion joint angle (b) must now be corrected to a c c o m m o d a t e the 90 degree orientation of the insertion lever arm (D), the corrected insertion angle is calculated as ABS ( 9 0 - b ) . To correct for the lever arms being on opposite sides of the intervening segment we must add 90 ° to the corrected insertion angle, which is represented as (c) in (4). Once the stage has been reached, muscle length in this plane can be calculated as: x/ D 2 + C 2 - 2 D C cos (c).

these main axes, the joint angle must be corrected (Fig. 1). The joint angle is corrected simply by subtracting it from this 90 degree element, then taking the absolute value of the result. While the absolute values of the joint movement will change, the relative values will not. Corrected Joint Angle = ABS (90- Joint Angle). (7) Calculation o f directional muscle length Once the lever arm and joint angles have been corrected, the muscle length in a given direction can be calculated by using the law of cosines (Fig. 1). For this example, A = insertion lever arm, B = origin lever arm, c =joint angle in this plane; muscle length would then equal: Length = x/ A 2 + B 2 - 2 A B cos (c). Since the range of the joint in this plane is known, the joint angle (c) can be expressed as a range of values from the initial to the terminal joint angles. In this way, a distinct planar muscle length is calculated for each joint position. If no joint

218

a.A.

YOUNG

AND

M . A. A S H L E Y

movement occurs in this plane then the muscle length in this plane will be a constant. This model requires that the number of intervals produced for joint movement be held constant for all joints in the system being studied. This same approach can be taken with a two-joint muscle. In this case, however, the length of the intervening bony segment must be added to the calculated muscle length that lies in the same plane as the intervening segment. If the intervening segment is not orientated along a plane, its length in each direction can be calculated and this length added to the respective muscle lengths. If the two lever arms in a given direction lie on the opposite sides of the main axis of the insertion segment (e.g. the origin is dorsal to the axis while the insertion is ventral to the axis) then an additional correction is required. The joint angle used in the calculation of directional muscle length would be the corrected joint angle plus 90 degrees. The movements of the joint of origin are used to correct the lever arms, while the movements of the joint of insertion are used to determine muscle length. With this system the relative timing of the movement of these joints can be easily incorporated by varying the rate of change of their respective values. (8) Calculation of total muscle length Given the length of the muscle in all three directions, the uncorrected muscle length can be calculated by using a form of the pythagorean theorem. For example, if A = cranio-caudal muscle length, B = dorso-ventral muscle length, and C = medio-lateral muscle length, then total uncorrected muscle length will equal: Length =

x/A 2 + B 2 + C 2 .

This total muscle length can be corrected by subtracting the length of all tendinous elements in series with the muscle. The corrected muscle length will be calculated over as many intervals as were used in the calculation of muscle length in each direction. (9) Calculation of muscle force Muscle cross-sectional area is estimated by dividing muscle mass by total muscle length. Traditionally, it has been assumed that cross-sectional area is directly proportional to muscle force. Were this assumption valid, no muscle would show a bell-shaped length-tension curve; all muscles would show simply a straight line with negative slope when force is plotted against length. Our model incorporates a multiplicative correction factor to approximate a lengthtension curve or other aspects of the physiological basis of force production, including such factors as neural activation. The correction factor we use is the natural log (e) raised to the x power. The values for x may be varied by the researcher to produce curves of different shapes. For instance, if the initial value of x is -0.50, and x is first increased by 0.05 each interval for ten intervals, and then decreased by 0.05 each interval for ten intervals, a perfect bell-shaped curve will result. This correction factor curve will multiply the calculated cross-sectional area by 1.00 at its peak, and by 0.60 at its extremes to approximate the true length-tension values

RELATIVE MUSCLE FORCE PREDICTION

219

for the muscle. The value of x should not exceed 0-0, but clearly almost any shape may be produced by simply varying the initial value and the size and sign of the increment value. In this way, any form of length-tension curve, neural stimulation pattern, or other physiological property is incorporated into the model. The number of curve intervals should equal the number of joint angle intervals used. Thus, for a parallel-fibered muscle, total muscle force may be calculated by: Force -

(Muscle Mass) x (e x) (Total Muscle L e n g t h ) - (Tendon Length)"

The calculations are more complex for pinnate muscles (Fig. 2). Force is assumed to be proportional to the number of muscle fibers, which is assumed to be proportional to the cross-sectional area of the muscle. The cross-sectional area of an equivalent parallel-fibered muscle may be approximated by an ellipse, with the

A

E

FIG. 2. Measurements and calculations used in relation to pinnate muscles. A = muscle length; B = width of minor axis of origin; C = width of major axis of origin; D = length of muscle fiber; E = length of central tendon spanned by average muscle fiber; a (alpha) = angle of pinnation. If this were a pinnate muscle, its cross-sectional area would be equal to ~r(C/2)(B/2). In the case of the pinnate muscle illustrated, its cross-sectional area is equal to ~-(A/2)(B/2). The difference in relative number of muscle fibers is found by {(~'(C/2)(B/2)[sin (a)])/Tr(A/2)(B/2)}. The muscle fiber length (D) is equal to (C/2)/sin (a). The length of the tendon spanned by this fiber is found by D cos a). Since the length of the tendon spanned by the muscle fiber is a constant proportion of the total muscle length a constant (K) can be defined: K = E/A. Lastly, these equations can be combined to allow for a dynamic modelling of pinnation angle: cos (a) = (A x K)/D.

220

B.A.

YOUNG

AND

M.

A. ASHLEY

major axis being the longest measurement taken in section (1), and the minor axis being the shortest. The cross-sectional area of this muscle would then be: Area = lr x (Major Axis/2) x (Minor Axis/2). If the muscle in question is bipinnate rather than unipinnate, the calculated area should be divided by two. The area of insertion of a pinnate muscle on a central tendon is approximated by a second ellipse having the same minor axis as the ellipse described above, and a major axis equal to the muscle length. The area of this second ellipse would be proportional to the number of muscle fibers if the fibers were oriented at 90 degrees to the tendon. When calculating the relative increase in number of muscle fibers, pinnation angle must be incorporated: Relative Fiber # -

(Insertion Area Pinnate)x sin(Pinnation Angle) Cross-Sectional Area if Parallel

This relative increase in the number of muscle fibers represents a portion of the advantage of being pinnate; for a bipinnate muscle, this quantity should be multiplied by two (since fibers insert on both sides of the central tendon). This advantage is countered by a very real disadvantage: only a fraction of the force produced by the muscle fibers is exerted along the line of the tendon, in effect "wasting" a percentage of the total available muscle force. The percentage of force "wasted" will increase with increasing angles of pinnation. Because angle of pinnation changes with muscle length, it must be modelled dynamically to reflect changes in effective muscle force. Changes in muscle length can be equated to changes in the length of the central tendon associated with the contracting muscle fibers. The length of the section of the central tendon spanned by a single muscle fiber (Fig. 2) can be calculated from pinnation angle and muscle fiber length: Spanned Tendon Length = Fiber Length x cos (Pinnation Angle). Individual muscle fiber length is calculated for a unipinnate muscle as: Fiber L e n g t h -

Maximum Muscle Width sin (Pinnation Angle)

and for a bipinnate muscle as: Fiber L e n g t h -

½Maximum Muscle Width sin (Pinnation Angle)

This calculated length represents the maximum muscle fiber length at rest. Since pinnate muscles are fusiform in appearance this maximum muscle fiber length can not be taken as representing the "typical" length of the fibers in the muscle. To account for the shorter fiber lengths at the ends of the muscle an "average" fiber length is calculated by taking this maximum length and multiplying by 0.80. The portion of the central tendon spanned by a single muscle fiber is related to the total length of the muscle by a proportionality constant, K: K=

Spanned Tendon Length Resting Length of Muscle"

RELATIVE

MUSCLE

FORCE

PREDICTION

221

Since the fibers of a pinnate muscle do not appreciably change in length during contraction (the pinnation angle changing instead), the angle of pinnation may be calculated for any muscle length: cos (Pinnation Angle) -

(Length of Muscle) x K (Muscle Fiber Length)"

Total muscle force may now be calculated for any muscle length by the formula: Force = cos (Pinnation Angle) x

Mass x [(R + Mass)/Mass] x Q Muscle Length

This equation is similar to the one used to calculate force for a parallel fibered muscle except it incorporates two aspects of the advantage of being pinnate (the relative increase in the number of fibers (R), and the packing advantage (Q) which is the number of pinnation axes squared) as well as the major disadvantage of being pinnate, that being the angle of pinnation. While muscle mass, R, and Q are held constant, pinnation angle will change with each new calculated muscle length. (10) Calculation of muscle force in each direction Since total muscle force is inversely proportional to total muscle length, it follows that the force in each direction will be inversely proportional to the muscle length in that direction. Muscle force in each direction may be calculated by dividing the muscle length in that direction (section 7) by the total muscle length (section 8), then multiplying the quotient by total muscle force (section 9): Directional Force =

Directional Length x Total Force. Total Length

This equation uses uncorrected muscle length (including tendon length), and is repeated for all intervals in all three directions.

(11) Calculation of the angle of insertion Insertion angle calculation requires values for insertion joint movement, directional planar muscle length (section 7), and the corrected lever arm (section 5): sin (Insertion Angle) -

(Corrected Lever Arm) x sin (Joint Movement) Directional Muscle Length

(12) Calculation of muscle torque Muscle force is in and of itself not a useful value in analyses of movement, since muscles produce movement primarily by exerting torques upon bony elements. Torque is defined as the component of muscle force perpendicular to the main axis of movement multiplied by the distance from the muscle insertion to the pivot point. Calculation of torque requires directional muscle force (section 10), the angle of

222

S.A.

YOUNG

AND

M. A. A S H L E Y

muscle insertion, and the insertion lever arm lying parallel to the long axis of the insertion element. Torque is given by: Torque = sin (Insertion Angle) x (Directional Force) x (Lever Arm). The muscle will also produce a torque along the axis oriented at 90 degrees to this axis. The magnitude of this second torque is: Torque = cos (Insertion Angle) x (Directional Force) x (Lever Arm). (13) Calculation of torsion In addition to producing torque along the long axis of the bone, the directional muscle force will produce torsions in the two other planes if the insertion is offset from the bone's long axis. The formula given above for the calculation of torque may be used to calculate torsions with the substitutions of the remaining two lever arms (in the two other directions) for the lever arm along the main axis of the bone. TESTING PROCEDURES FOR THE PROPOSED MODEL The accuracy of the proposed model was tested by comparing its output with literature values. Since the cat is the major research animal used in muscle physiology, we tested our model using four muscles from the cat hindlimb. Specifically, we used this model to predict the maximum forces produced by tibialis anterior (TA; a parallel fibered flexor of the foot), extensor digitorium longus (EDL; a unipinnate flexor of the foot that spans the knee joint), soleus (Sol; a unipinnate extensor of the foot), and medial gastrocnemius (MG; a multipinnate extensor of the foot that spans the knee joint). The testing of only maximum muscle forces was a limitation placed upon us by the literature, since no detailed comparative studies of torque are available. For consistency, all measurements of lever arm distance, muscle mass, tendon length, tibial length, and pinnation angle were taken from both legs of two preserved specimens (Table 1). We recognize that this small sample size and the use of TABLE 1

Input values for testing the proposed model Parameter D - V insertion lever A - P insertion lever M - L insertion lever D - V origin lever A - P origin lever M - L origin lever Muscle mass (g) Length of tendons Angle of pinnation Maximum width o f origin Minimum width of origin

EDL

TA

0.37 6"58 0-01 0.42 0.21 0.51 1.6 5.74 10 0.45 0.35

0.06 2-21 0-843 8.98 0.508 0.388 3.26 3.987 ----

MG 0.273 1.87 0.01 0.508 0.778 0.528 5.15 4.07 30 2-00 0-45

Sol 0.273 1.87 0-01 7.56 0.614 0.87 1.94 3'2 10 0.75 0.34

RELATIVE

MUSCLE

FORCE

PREDICTION

223

preserved material could possibly affect the output of the model. It is our opinion that the possible error due to these factors is not of-a magnitude to obscure the accuracy o f the model. For this test M G was treated as a bipinnate muscle. As no three-dimensional description of the kinematics of the cat hindlimb are available, we were restricted to examining joint movements in one plane. All values for joint movement were taken from Rasmussen et al. (1978); maximum force values for M G and Sol were taken from Walmsley et al. (1978), and those for TA and E D L were taken from Goslow et al. (1977). Once again it should be noted that the nature of this testing was also determined by the inherent requirement of our model to compare two muscles. P R O C E D U R E S FOR P E R T U R B A T I O N A N A L Y S I S

To further test this proposed model, and to explore some of the functional implications of the musculo-skeletal system, a perturbation analysis was performed. A baseline system was developed using the data presented in Table 2. Different properties of this system were then systematically varied (a single component at a time) and the resulting change in force or torque charted. Thus this perturbation analysis is actually creating a series of separate muscles, each of which has a different state of the selected variable, the relative forces of which are then compared. TABLE 2 Baseline values for perturbation analysis D - V insertion lever arm A - P insertion lever arm M - L insertion lever arm D - V origin lever arm A - P origin lever arm M - L origin lever arm Muscle mass T en don length Pinnation angle (where applicable) M a x i m u m width of origin (where applicable) M i n i m u m width of origin (where applicable) D - V range of joint movement A - P range of joint m o v e m e n t M - L range of joint m o v e m e n t

0-2 cm 5.0 cm 0.3 cm 4.5 cm 0.6 cm 0.4 cm 4.0g 2.5 cm 25* 1.2 cm 0.5 cm 60*-70* 30*-35* 40°_45 °

Except where noted, all perturbation analyses were conducted using both a parallel fibered and bipinnate muscle. A total of eight perturbations analyses were conducted, with the following parameters: (i) variations in pinnation, contrasting parallel fibered with unipinnate and bipinnate; (ii) effects of alteration in pinnation angle, varying the pinnation angle of a bipinnate muscle from 10 to 48 degrees; (iii) varying muscle mass from 2.0 to 6.0 g; (iv) alterations in muscle length produced by changing the tendon length from 0.5 to 4.5 cm; (v) changes in range of joint motion, by changing the joint ranges 20 degrees above and below the resting level, using a five