Sensitivity of model predictions of muscle function to changes in moment arms and muscle–tendon properties: A Monte-Carlo analysis

Sensitivity of model predictions of muscle function to changes in moment arms and muscle–tendon properties: A Monte-Carlo analysis

Journal of Biomechanics 45 (2012) 1463–1471 Contents lists available at SciVerse ScienceDirect Journal of Biomechanics journal homepage: www.elsevie...

576KB Sizes 2 Downloads 40 Views

Journal of Biomechanics 45 (2012) 1463–1471

Contents lists available at SciVerse ScienceDirect

Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com

Sensitivity of model predictions of muscle function to changes in moment arms and muscle–tendon properties: A Monte-Carlo analysis David C. Ackland n, Yi-Chung Lin, Marcus G. Pandy Department of Mechanical Engineering, University of Melbourne, Parkville, Victoria 3010, Australia

a r t i c l e i n f o

abstract

Article history: Accepted 8 February 2012

Hill-type muscle models are commonly used in musculoskeletal models to estimate muscle forces during human movement. However, the sensitivity of model predictions of muscle function to changes in muscle moment arms and muscle–tendon properties is not well understood. In the present study, a three-dimensional muscle-actuated model of the body was used to evaluate the sensitivity of the function of the major lower limb muscles in accelerating the whole-body center of mass during gait. Monte-Carlo analyses were used to quantify the effects of entire distributions of perturbations in the moment arms and architectural properties of muscles. In most cases, varying the moment arm and architectural properties of a muscle affected the torque generated by that muscle about the joint(s) it spanned as well as the torques generated by adjacent muscles. Muscle function was most sensitive to changes in tendon slack length and least sensitive to changes in muscle moment arm. However, the sensitivity of muscle function to changes in moment arms and architectural properties was highly muscle-specific; muscle function was most sensitive in the cases of gastrocnemius and rectus femoris and insensitive in the cases of hamstrings and the medial sub-region of gluteus maximus. The sensitivity of a muscle’s function was influenced by the magnitude of the muscle’s force as well as the operating region of the muscle on its force-length curve. These findings have implications for the development of subject-specific models of the human musculoskeletal system. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Subject-specific model Muscle model Tendon slack length Optimum muscle-fiber length Peak isometric muscle force

1. Introduction The behavior of a musculoskeletal model is influenced by the values of the parameters used to describe muscle–tendon actuation. In a Hill-type muscle model, muscle force is influenced by the values of the mechanical and architectural properties of the muscle–tendon unit, including the peak isometric force and corresponding fiber length of the muscle and the slack length of the tendon (Zajac, 1989), all of which depend on the size, gender, age and neuromuscular condition of the subject (Thelen, 2003; Piazza, 2006). Muscle moment arms are also person-specific (Scheys et al., 2008) and can affect muscle recruitment and hence, muscle forces and joint moments (Herzog, 1992; Raikova and Prilutsky, 2001; Nagano and Komura, 2003). Understanding the sensitivity of model predictions of muscle function to changes in moment arms and muscle–tendon parameters is important because many of these parameters cannot be obtained by direct measurement in vivo, and also because their values can vary widely in the literature. For example, cadaveric measurements of the physiological cross-sectional area (PCSA) of

n

Corresponding author. Tel.: þ613 8344 0405; fax: þ 613 9347 8784. E-mail address: [email protected] (D.C. Ackland).

0021-9290/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2012.02.023

rectus femoris are reported to be between 13.5 cm2 and 43.0 cm2 (Friederich and Brand, 1990; Ward et al., 2009), whereas in vivo data obtained from magnetic resonance (MR) imaging suggest a PCSA closer to 50.0 cm2 and a nearly two-fold difference in muscle volume between males and females (Tate et al., 2006). Calculation of the peak isometric force of rectus femoris (product of muscle PCSA and specific tension) based on these measurements, and assuming a specific tension value of 350 kPa (Zajac, 1989), results in values ranging from 473 N to 1750 N (see also Pandy and Andriacchi (2010) for a review). To better understand model sensitivity to muscle–tendon parameter variations, changes in muscle forces and joint torques have been assessed in response to specific parameter perturbations during forward simulations of gait. Previous studies have found that muscle–tendon force and muscle-generated torque are most sensitive to changes in tendon slack length (Scovil and Ronsky, 2006; Redl et al., 2007). Muscle–tendon force may also be sensitive to changes in the peak isometric force and corresponding fiber length of muscle (De Groote et al., 2010), depending on the specific muscle under consideration (Xiao and Higginson, 2010). While the sensitivities of muscle force and joint torque to changes in muscle–tendon parameters have been quantified, little is known about how variations in muscle–tendon parameters

1464

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

affect model predictions of muscle function during movement; specifically, the ability of the lower limb muscles to accelerate the center of mass (COM) of the body in the vertical, fore-aft and mediolateral directions, or alternatively, to generate vertical support, modulate forward progression and provide mediolateral stability, respectively. In multi-joint musculoskeletal systems, the sensitivity of muscle contributions to whole-body motion cannot be directly inferred from muscle-force or joint-torque sensitivity data because of the complex relationship between a muscle’s force and its contribution to the joint accelerations, which is determined primarily by the structure of the inverse mass matrix (Zajac and Gordon, 1989; Yu et al., 2011). The aim of the present study was to use a three-dimensional muscle-actuated model of the body in conjunction with the MonteCarlo method to investigate the sensitivity of model predictions of muscle function to changes in moment arms and muscle–tendon properties during gait. Unlike previous methods that have evaluated variations in a given dependent variable (e.g., muscle force) resulting from single percentage changes in an independent variable (e.g., tendon slack length) (Scovil and Ronsky, 2006; Redl et al., 2007), the Monte-Carlo method provides an approach for quantifying the influence on a given dependent variable of a broad distribution of values of a specific independent variable. Importantly, the Monte-Carlo method can be used to measure effects due to changes in multiple independent variables, for example, in determining how combined changes in both a muscle’s fiber length and tendon slack length may influence model predictions of muscle function.

2. Methods 2.1. Gait experiments Gait data were recorded for five healthy male subjects (age: 26 74 yrs; mass: 707 5 kg; height: 178 7 4 cm). Joint motion and ground reaction forces were measured simultaneously as each subject walked on level ground at his preferred speed. Three-dimensional locations of retro-reflective markers affixed to each subject’s body were measured using a 9-camera motion analysis system (Vicon, Oxford Metrics Ltd., Oxford), and ground force data were recorded using three instrumented force platforms (Advanced Mechanical Technology Inc., Watertown, MA) embedded in the laboratory floor. Video and analog force plate data were sampled at 120 Hz and 1080 Hz, respectively. Marker data and ground force data

were filtered with a low-pass, fourth-order Butterworth filter using a cut-off frequency of 4 Hz and 40 Hz, respectively. Details of the experimental protocol are given in Pandy et al. (2010).

Fig. 1. (A) Diagram illustrating a box plot of representative results obtained from a Monte-Carlo simulation. The box plot represents the statistical distribution of a dependent variable, where a red horizontal line represents 50th percentile data; upper and lower edges of the box represent 75th and 25th percentile data, respectively; upper and lower bars represent maximum and minimum data values, respectively; and red crosses represent outliers. (B) Sample box plot of data obtained from a hypothetical Monte-Carlo simulation. In the example shown, the dependent variable is RMS difference in muscle contribution to the vertical acceleration of the body’s center-of-mass (COM), expressed as a percentage of the muscle’s nominal contribution to the vertical COM acceleration. The COM acceleration is more sensitive to changes in optimal muscle-fiber length than peak isometric muscle force, since larger differences in COM acceleration result from changes in optimal muscle-fiber length (i.e., 50th percentile data of optimal muscle-fiber length are larger than that of peak isometric muscle force). In addition, COM acceleration is distributed more widely due to variations in optimal muscle-fiber length (i.e., box plot of optimal muscle-fiber length is longer than that of peak isometric muscle force), indicating that greater changes in COM acceleration result from variations in optimal musclefiber length. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Summary of the total number of simulations for each Monte-Carlo analysis. Monte-Carlo analyses were performed for each muscle–tendon unit by perturbing nominal values assumed for peak isometric muscle force, optimal muscle-fiber length, tendon slack length, a combination of all three muscle–tendon parameters, and muscle M

T

moment arm. Symbols appearing in the table are as follows: lo , optimal muscle-fiber length; ls , tendon slack length; and F M o , peak isometric muscle force. ‘‘Combined’’ represents changes made to all three muscle–tendon parameters simultaneously. Muscle acronyms are defined in the text.

Fore-aft direction o lm

SOL

GAS

1250

VAS

RF

HAMS

GMAXL

GMAXM

GMEDA

GMEDP

ILPSO

1330

1010

1020

1120

1050

1280

1200

1000

1580

lt F om Combined Moment arm

1680

1030

1290

1090

2520

1030

1010

1190

1040

1000

1090 1610 1160

1210 1470 1050

1280 1070 1020

1360 1010 1080

1020 1520 1070

1220 1520 1110

1490 1250 1040

1100 1380 1190

1010 1570 1500

1030 1680 1020

Vertical direction o lm

1110

1060

1380

1080

1350

1080

1400

1170

1170

1180

1100

1120

1080

1020

1170

1010

1130

1020

1610

1130

1380 1300 1200

1320 1050 1070

1040 1540 1410

1020 1050 1330

1410 1280 1400

1210 1190 1100

1020 1890 1540

1120 1150 1130

1090 1270 1000

1100 1090 1170

Mediolateral direction o 1310 lm s 1040 lt

1030

1010

1000

1010

1200

1170

1030

1480

1120

1040

1170

1110

1010

1180

1210

1040

1270

1160

F om Combined Moment arm

1030 1130 1010

1140 1030 1040

1050 1510 1050

1020 1070 1020

1210 1140 1100

1340 1350 1110

1100 1000 1470

1090 1020 1400

1190 1290 1060

s

s

lt F om Combined Moment arm

1010 1070 1290

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

1465

Fig. 2. Variations in muscle forces obtained from Monte–Carlo analyses of the stance phase of gait. Plots are given for individual lower limb muscles obtained from simulations that varied peak isometric muscle force, optimal muscle-fiber length, tendon slack length, a combination of all three muscle–tendon parameters, and muscle moment arm. Gray shaded regions represent the range ( 7 one standard-deviation) of muscle force values, whereas black solid lines represents the nominal muscle force values. Muscle symbols as defined in the text. 2.2. Calculations of leg-muscle function A three-dimensional musculoskeletal model was used to study the sensitivity of changes in muscle moment arms and muscle–tendon parameters to model predictions of (i) muscle forces, (ii) muscle contributions to net joint torques, and (iii) muscle contributions to COM accelerations. The model was identical in structure with that described by Anderson and Pandy (1999). The skeleton was represented as a 10-segment, 23 degree-of-freedom linkage system. The head, arms and torso were modeled as a single rigid body, which articulated with the pelvis via a balland-socket back joint. Each hip was modeled as a ball-and-socket joint, each knee as a hinge, each ankle-subtalar complex as a universal joint, and each metatarsal joint as a hinge. Segment lengths, muscle attachment sites and inertial properties of the generic model were scaled to the average anthropometry of the five subjects. Fiftyfour muscle–tendon units actuated the model, each unit represented as a Hill-type muscle in series with tendon. The force–length and force–velocity properties of each muscle–tendon actuator are given in Anderson and Pandy (1999). A nominal simulation of walking was generated by applying the average joint motion and ground reaction forces measured for the five subjects to the scaled-generic model. Inverse dynamics was used to calculate the net moments developed about the back, hip, knee and ankle joints for the stance phase of walking. The net joint moments were then decomposed into muscle forces by solving a static optimization

problem that minimized the sum of the squares of all muscle activations subject to the physiological bounds imposed by each muscle’s force–length–velocity properties (Anderson and Pandy, 2001b). This performance criterion was used because it has been shown to produce lower limb muscle forces for walking which are consistent in sequence and timing with measured EMG activity (Anderson and Pandy, 2001a). A pseudo-inverse decomposition method was used to determine the nominal contribution of each muscle force to the COM acceleration (Lin et al., 2011). First, the contribution of each muscle force to the ground reaction force was computed. Five foot–ground contact points were defined on the sole of the foot in the model, four points located at the perimeter of the hind-foot and one at the distal tip of the toes segment. At each time instant, the dynamical equations of motion were solved for each muscle’s contribution to the joint angular accelerations and the ground reaction force acting at each of the five foot–ground contact points. Once each muscle’s contribution to the ground reaction force was found, Newton’s 2nd law of motion was used to determine its contribution to the accelerations of the COM in the vertical, fore-aft and mediolateral directions. 2.3. Monte-Carlo analyses The Monte-Carlo method was used to evaluate the sensitivity of calculations of muscle function to changes in the values of muscle moment arm, peak isometric

1466

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

muscle force, optimal muscle-fiber length, and tendon slack length. The following muscles were selected for analysis as they have been shown to contribute most significantly to the acceleration of the COM in normal walking (Liu et al., 2006; Pandy et al., 2010): soleus (SOL), gastrocnemius (GAS), vasti (VAS), rectus femoris (RF), hamstrings (HAMS), medial (GMAXM) and lateral (GMAXL) sub-regions of gluteus maximus, anterior (GMEDA) and posterior (GMEDP) sub-regions of gluteus medius, and iliopsoas (ILPSO). For each actuator and for each of the three muscle–tendon parameters, a random parameter value was chosen between þ 10% and  10% of the nominal value of the parameter. The static optimization problem was re-solved and muscle contributions to COM acceleration then computed in the vertical, fore-aft and mediolateral directions for the stance phase of gait. Similarly, a new time history of moment arm values was created by randomly perturbing each muscle’s moment arm between þ10% and  10% of its nominal value throughout stance. The optimization problem was again solved, and muscle contributions to COM acceleration computed. For each optimization solution, the root-mean-square (RMS) difference was calculated between the muscle contribution associated with the perturbed value of the parameter and that associated with the nominal value. A convergence criterion was defined as a stopping rule for each Monte-Carlo analysis (Valero-Cuevas et al., 2003). A given percentage of the simulations was examined to ensure that the mean and coefficient of variation (i.e., 100  standard deviation/mean) were within a given convergence band. Specifically, the convergence criterion was reached when the mean and coefficient of variation across the final 10% of the simulations were within 2% of the overall mean and coefficient of variation. The total number of simulations performed for each Monte-Carlo analysis is given in Table 1. To quantify the sensitivity of the model to combined changes in muscle– tendon parameters, an additional Monte-Carlo analysis was performed by randomly varying peak isometric muscle force, optimal muscle-fiber length, and tendon slack length simultaneously; that is, by simultaneously perturbing each parameter between þ 10% and  10% of its nominal value. The results were displayed as box plots (see Fig. 1).

3. Results The calculated values of muscle forces were highly sensitive to variations in tendon slack length, particularly in the cases of GAS and RF (Fig. 2). By comparison, muscle forces were relatively insensitive to variations in peak isometric muscle force, optimal muscle-fiber length and moment arm. Varying values of all three muscle–tendon parameters simultaneously resulted in similar magnitude changes in muscle forces to those obtained from varying tendon slack length alone. In general, varying the moment arm or properties of a given muscle resulted in changes in the contribution of that muscle to the torque generated about the joint spanned by the muscle; however, the joint torques generated by adjacent, non-perturbed muscles also changed substantially (Fig. 3). For example, varying the peak isometric muscle force, optimal muscle-fiber length, tendon slack length, or moment arm of ILPSO led to a greater variation in the hip flexion torque generated by RF than that generated by ILPSO. Similarly, varying the peak isometric muscle force, optimal muscle-fiber length, or tendon slack length of SOL led to a greater variation in the ankle plantarflexion torque generated by PFEV than that generated by SOL. Muscle function was substantially more sensitive to variations in tendon slack length than to variations in peak isometric muscle force, optimal muscle-fiber length, or moment arm (Fig. 4 and Table 2). When the three muscle–tendon parameters were varied simultaneously, the magnitude changes in muscle function were similar to those obtained by varying tendon slack length alone. The sensitivity of muscle function to variations in muscle–tendon parameters was muscle-specific. Changes in muscle function resulting from variations in tendon slack length of GAS and RF were at least three times larger than those due to variations in peak isometric muscle force, optimal muscle-fiber length, and moment arm of any other muscle (Fig. 4). In contrast, calculations of muscle function were insensitive to variations in muscle–tendon parameters and the moment arm of HAMS (Fig. 4), which generated forces transiently during early stance (Fig. 2). The function of GMEDA was also insensitive to variations in muscle–tendon parameters and muscle moment arm.

Fig. 3. Sensitivity of muscle-generated joint torques to changes in muscle–tendon properties and moment arms for the stance phase of gait. Box plots of RMS differences in a muscle’s contribution to a net joint torque (expressed as a percentage of the muscle’s nominal contribution to the net joint torque) are given for SOL (ankle plantarflexion), RF (knee extension), GMAXL (hip extension), GEMDP (hip abduction) and ILPSO (hip flexion) for simulations that varied peak isometric muscle force, optimal muscle-fiber length, tendon slack length, a combination of all three muscle–tendon parameters, and muscle moment arm. For each muscle, boxes are plotted for that muscle (see labels at the top of each diagram) along with the remaining two most significant muscle contributors to RMS differences in joint torque. Muscle symbols appearing in the charts are defined in the text and as follows: DFIN, tibialis anterior and extensor hallucis longus combined; ADLB, adductor longus brevis; PFEV, peroneus brevis and peroneus longus combined.

For SOL, GAS, VAS, RF and HAMS, a substantially greater distribution of muscle-fiber lengths was observed when tendon slack length was varied compared to optimal muscle-fiber length (Fig. 5(A)). Most notably, variations in the tendon slack lengths of

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

1467

Fig. 4. Sensitivity of muscle function to changes in muscle–tendon properties and moment arms for the stance phase of gait. Box plots of RMS differences in a muscle’s contribution to the vertical center-of-mass (COM) acceleration (expressed as a percentage of the muscle’s nominal contribution to the vertical COM acceleration) are given for individual lower limb muscles obtained from simulations that varied peak isometric muscle force, optimal muscle-fiber length, tendon slack length, a combination of all three muscle–tendon parameters, and muscle moment arm.

SOL and GAS resulted in a distribution of muscle-fiber lengths spanning the majority of the ascending and descending regions of the muscles’ force–length curves (Fig. 5(B)), reflecting a high sensitivity of muscle function to changes in tendon slack length. In contrast, variations in both the tendon slack lengths and optimal fiber lengths of GMAXL, GMAXM, GMEDA, GMEDP and ILPSO resulted in a much smaller distribution of muscle-fiber lengths. The fiber lengths in a number of these muscles, particularly ILPSO, were confined to just the upper ascending and plateau regions of the muscles’ force–length curves, reflecting a relatively low sensitivity of muscle function to changes in the tendon slack lengths and optimal muscle-fiber lengths of these actuators.

4. Discussion Hill-type muscle models are often combined with articulated rigid-body models of the skeleton to predict muscle and joint loading during motor tasks such as walking (Anderson and Pandy, 2001a; Viceconti et al., 2006; Erdemir et al., 2007; Correa et al., 2010; Steele et al., 2010). However, it is not clear how the functional role of a given muscle is affected by variations in the values assumed for the moment arms and architectural properties of the muscle. In this study, we used a Monte-Carlo approach to evaluate the sensitivity of the functional roles of each major lower limb muscle to changes in muscle moment arms and muscle–tendon

parameters. We found that muscle function calculations were most dependent on tendon slack length and least affected by changes in moment arm. The sensitivity of muscle function to variations in tendon slack length was similar to that resulting from simultaneous variations in tendon slack length, peak isometric muscle force and optimal muscle-fiber length. The sensitivity of muscle function was also muscle-specific and dependent on the magnitude of the muscle’s force as well as the operating region of the muscle on its force–length curve. Perturbing the architectural properties or moment arm of an individual muscle resulted in changes in the contribution of that muscle to the torque generated about the joint spanned by the muscle; however, the torques generated by other non-perturbed muscles also changed substantially as a result of the parameter variations (Fig. 3). For example, perturbing the peak isometric force or moment arm of ILPSO led to greater variations in the hip-flexion torque of RF than that of ILPSO. Increases in the forces of muscles other than those perturbed during simulations were a consequence of the objective function assumed in the optimization problem solved for walking (i.e., minimization of the sum of the squares of all muscle activations). For a given muscle group to achieve a desired joint torque, activating an agonist muscle with a larger peak isometric force and/or a larger moment arm may have been more ‘cost-efficient’ than increasing the activation of a perturbed muscle with a smaller peak isometric force and/or a smaller moment arm.

1468

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

Table 2 Summary of Monte-Carlo simulation results illustrating the sensitivity of muscle function to changes in muscle–tendon properties and moment arms for the stance phase of gait. For each simulation, the mean and range of RMS differences in center-of-mass (COM) acceleration (as a percentage of the nominal acceleration) in the vertical, fore-aft, and mediolateral directions are given. See Table 1 for symbol definitions. SOL Fore-aft COM acceleration Mean o lm Max Min

GAS

VAS

RF

HAMS

GMAXL

GMAXM

GMEDA

GMEDP

ILPSO

(%) 0.4 0.9 0.0

0.6 1.3 0.0

0.5 1.6 0.0

2.6 6.8 0.0

0.2 0.5 0.0

2.1 4.1 0.0

0.9 2.1 0.0

0.2 0.5 0.0

1.0 2.7 0.0

2.6 5.6 0.0

lst

Mean Max Min

8.7 33.0 0.0

17.1 83.7 0.0

3.5 17.1 0.0

14.3 45.2 0.0

2.0 9.2 0.0

2.7 5.6 0.0

1.1 3.0 0.0

0.3 0.6 0.0

1.3 3.9 0.0

2.8 5.9 0.0

o Fm

Mean Max Min

0.9 2.0 0.0

1.6 3.4 0.0

0.2 0.4 0.0

2.5 5.3 0.0

0.5 1.1 0.0

1.1 2.3 0.0

1.6 3.3 0.0

0.6 1.3 0.0

1.9 3.8 0.0

1.4 3.0 0.0

Combined

Mean Max Min

9.6 73.3 0.1

17.7 166.1 0.1

3.2 23.1 0.0

14.9 47.6 0.5

2.5 14.8 0.0

3.5 11.9 0.1

2.1 8.7 0.0

0.7 2.4 0.0

2.6 8.9 0.1

3.2 14.8 0.1

Moment arm

Mean Max Min

1.3 2.5 0.1

0.6 1.6 0.1

1.8 3.9 0.0

1.7 3.5 0.3

0.4 0.8 0.0

0.7 1.5 0.0

0.8 1.6 0.0

0.6 1.3 0.0

0.3 0.8 0.0

1.3 2.6 0.2

(%) 0.4 0.8 0.0

1.0 4.9 0.0

0.5 1.5 0.0

2.1 5.0 0.0

0.1 0.2 0.0

2.4 6.7 0.0

0.8 1.8 0.0

0.1 0.4 0.0

0.7 2.2 0.0

1.6 3.5 0.0

Vertical COM acceleration Mean lom Max Min lst

Mean Max Min

7.2 30.8 0.0

18.4 84.3 0.0

2.8 15.1 0.0

11.8 40.2 0.0

0.8 4.6 0.0

3.2 9.9 0.0

1.0 2.7 0.0

0.2 0.5 0.0

0.7 2.0 0.0

1.8 3.7 0.0

o Fm

Mean Max Min

0.8 1.8 0.0

2.0 4.3 0.0

0.2 0.5 0.0

2.3 4.7 0.0

0.2 0.4 0.0

1.7 3.6 0.0

1.4 2.9 0.0

0.5 1.1 0.0

2.6 5.4 0.0

1.2 2.7 0.0

Combined

Mean Max Min

8.0 56.6 0.0

20.1 154.5 0.3

3.1 21.9 0.0

11.9 42.9 0.1

1.0 8.3 0.0

4.5 21.9 0.1

2.0 7.2 0.0

0.5 1.8 0.0

3.0 9.2 0.1

2.9 29.8 0.1

Moment arm

Mean Max Min

1.4 2.9 0.1

0.5 1.2 0.1

1.7 3.8 0.0

1.5 3.1 0.2

0.2 0.8 0.0

0.8 1.5 0.0

0.8 2.3 0.0

0.3 0.7 0.0

0.6 1.5 0.0

1.0 1.9 0.1

Mediolateral COM acceleration (%) Mean 0.4 1.0 Max 0.8 2.0 lom Min 0.0 0.0

0.4 1.4 0.0

0.6 1.4 0.0

0.1 0.3 0.0

10.3 23.0 0.1

1.0 2.3 0.0

0.2 0.5 0.0

1.0 2.4 0.0

1.8 3.9 0.0

lst

Mean Max Min

6.4 31.7 0.0

13.5 46.5 0.0

2.6 13.9 0.0

3.5 13.8 0.0

1.0 4.4 0.0

15.4 30.5 0.0

1.2 3.4 0.0

0.3 0.7 0.0

1.3 3.5 0.0

2.0 4.0 0.0

o Fm

Mean Max Min

0.8 1.9 0.0

1.9 4.1 0.0

0.2 0.4 0.0

0.8 1.7 0.0

0.2 0.4 0.0

2.5 4.9 0.0

1.8 3.7 0.0

0.6 1.2 0.0

2.1 4.6 0.0

1.1 2.3 0.0

Combined

Mean Max Min

7.1 53.4 0.1

15.8 89.8 0.3

3.0 20.4 0.0

4.0 18.4 0.1

1.1 7.5 0.0

17.0 49.1 0.4

2.5 9.1 0.0

0.7 2.6 0.0

2.9 9.4 0.1

2.8 28.8 0.1

Moment arm

Mean Max Min

1.4 2.9 0.1

0.4 0.6 0.1

1.7 3.6 0.0

1.2 2.5 0.1

0.2 0.9 0.0

4.8 9.7 0.3

0.9 1.9 0.0

1.2 2.4 0.0

0.5 1.1 0.1

1.1 2.3 0.1

The sensitivity of muscle function to changes in muscle–tendon parameters was muscle-specific. For example, muscle function was highly sensitive to muscle–tendon parameter changes in the cases of GAS and RF, but relatively insensitive to such changes in the cases of HAMS and GMEDA (Fig. 4). In general, however, a muscle’s function was most sensitive to changes in the values of tendon slack length. This finding extends the results of previous investigations which have shown that calculations of muscle forces and joint torques for walking are sensitive to changes in the values assumed for tendon slack length (Redl et al., 2007; De Groote et al., 2010). Unfortunately, in practice, direct measurement of tendon slack length is problematic due to the difficulty in identifying the aponeurotic part of the tendon

from the muscle belly, and because the tendon may be elongated when post-mortem measurements are made. When developing a musculoskeletal model for simulations of human movement, care should be taken in applying tendon slack length values for specific parameter-sensitive muscles such as GAS and RF, since the assumed parameter values may greatly influence the overall behavior of the model and hence calculations of muscle function. Changes in tendon slack length can substantially affect the operating region of a muscle on its force–length curve, and therefore the force-generating potential of the muscle. However, tendon slack length for HAMS was changed considerably during the simulations (Fig. 5(A)), with a small corresponding effect on muscle function

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

1469

Fig. 5. (A) Range of normalized muscle-fiber lengths obtained from Monte–Carlo simulations of the stance phase of gait. Orange shaded regions represent data for simulations that varied tendon slack-length; gray shaded regions represent data for simulations that varied optimal muscle-fiber length; and black solid lines represent the nominal muscle-fiber lengths. (B) Mean, maximum and minimum normalized muscle-fiber lengths of selected muscles obtained from Monte-Carlo analyses of the stance phase of gait, shown relative to each muscle’s force–length curve. Black circles and bars represent mean and maximum/minimum normalized fiber lengths, respectively, for the nominal case; orange squares and bars represent mean and maximum/minimum normalized fiber lengths for simulations that varied tendon slack-length; gray diamonds and bars represent mean and maximum/minimum normalized fiber lengths for simulations that varied optimal muscle-fiber length. Muscle symbols are as M follows: L, muscle-fiber length; LM o , optimal muscle-fiber length; F, muscle force; F o , peak isometric muscle force. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(Fig. 4 and Table 2). This is because muscle function was relatively insensitive to changes in muscle–tendon parameters and moment arms of muscles that were either weakly active or active for only a short period of time. Because HAMS was only briefly active during early stance, the functional role of this muscle was not greatly affected by changes in its architectural properties or moment arms (Figs. 2 and 4). Muscle moment arms are commonly measured using MR imaging in vivo (Wilson et al., 1999; Wilson and Sheehan, 2009) or by applying the tendon excursion technique ex vivo (Nemeth and Ohlsen, 1985; Ackland et al., 2010). Although the importance of accurately representing muscle paths in musculoskeletal models is widely recognized (Scheys et al., 2008; Correa et al., 2011), the sensitivity of lower limb muscle function to variations in moment arms has not been fully documented. Our results indicate that tendon slack length is a more significant determinant of the functional role of

a muscle–tendon unit than its moment arm. The effects of moment arm changes on a muscle’s force, the torque exerted by the muscle about a joint, and the functional roles of the muscle during walking were similar to the effects resulting from changes in the muscle’s peak isometric force and optimal-fiber length. It should be noted in the interpretation of these results that muscles were treated as ideal force generators (i.e., the force–length–velocity properties were neglected) in assessing model sensitivity to changes in moment arm. This approach eliminated muscle-fiber length change as a confounding variable in these analyses, since it was not necessary to alter the muscle’s path in order to perturb its moment arm. Nonetheless, modeling a muscle as an ideal force generator has been shown to have a small influence on the calculated values of muscle forces during gait (Anderson and Pandy, 2001b) (Fig. 2, compare nominal muscle force curves in plots of ‘Optimal muscle-fiber length’ versus plots of ‘Muscle moment arm’).

1470

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

Our results also suggest that the sensitivity of muscle function to changes in muscle–tendon parameters and moment arms is task-dependent. Muscles operating on the plateau of their force– length curve are generally less sensitive to variations in muscle architectural properties than those operating in the ascending or descending regions, since muscle–tendon parameter variations result in smaller changes in muscle force values in the plateau region. Thus, greater range of joint motion in a task such as running might push the operating region of a muscle off the plateau and into the ascending or descending region of its force–length curve, thereby increasing the sensitivity of the muscle’s function to variations in its architectural properties. Different motor tasks are also likely to generate different muscle activation patterns, and hence muscle forces, further affecting the sensitivities of these muscles to changes in their muscle–tendon parameters and moment arms. For example, HAMS generates force only during early stance in walking, but is active over the entire stance phase of running (Besier et al., 2009; Schache et al., 2012), indicating that calculations of HAMS’ function during running may be more sensitive to changes in the muscle’s architectural properties. There are limitations of our analysis which should be considered when interpreting the results. First, we investigated only the stance phase of the gait cycle because during this period the forces developed by the muscles are largest. The model may have responded differently to parameter variations had we analyzed other discrete periods of the gait cycle, such as terminal swing, or investigated other motor tasks, such as running. Second, we examined the sensitivity of a musculoskeletal model to parameter variations within 710% of their nominal values. While others have perturbed muscle parameters by up to 50% of their nominal values (Scovil and Ronsky, 2006), we found that tendon slack lengths exceeding þ10% had a tendency to push some muscles off their force–length curves. We believe the chosen range of perturbations is sufficient to demonstrate the relative changes in calculations of muscle function resulting from variations in the moment arms and architectural properties assumed for muscles. Third, results of the Monte-Carlo analyses were expressed as the RMS difference in COM acceleration, and were not normalized to the magnitude of the specific parameter changes (Scovil and Ronsky, 2006; Redl et al., 2007). However, we found that normalizing results of muscle–tendon parameter perturbations by the magnitude of the perturbed parameter did not significantly affect the muscle-specific trends in parameter sensitivity. Finally, the use of different objective functions may have influenced our predictions of model sensitivity to muscle–tendon parameters. While repeating our analyses with a different exponent of the cost-function (i.e., minimizing the sum of cubed muscle activations) had a negligible effect on both model predictions of muscle function and model parameter sensitivity, further analyses may be required to quantify model sensitivity when alternative performance criteria are used (e.g., minimizing the cost of transport).

Conflict of interest statement The authors do not have any financial or personal relationships with other people or organizations that could inappropriately influence their work.

Acknowledgments This work was funded by Australian Discovery Projects grants DP0878705 and DP1095366 and by an Australian Linkage Project grant LP0990369. Partial funding from a VESKI Innovation Fellowship awarded to M.G.P is also gratefully acknowledged.

References Ackland, D.C., Roshan-Zamir, S., Richardson, M., Pandy, M.G., 2010. Moment arms of the shoulder musculature after reverse total shoulder arthroplasty. Journal of Bone & Joint Surgery 92, 1221–1230. Anderson, F.C., Pandy, M.G., 1999. A dynamic optimization solution for vertical jumping in three dimensions. Computer Methods in Biomechanics and Biomedical Engineering 2, 201–231. Anderson, F.C., Pandy, M.G., 2001a. Dynamic optimization of human walking. Journal of Biomechanical Engineering 123, 381–390. Anderson, F.C., Pandy, M.G., 2001b. Static and dynamic optimization solutions for gait are practically equivalent. Journal of Biomechanics 34, 153–161. Besier, T.F., Fredericson, M., Gold, G.E., Beaupre, G.S., Delp, S.L., 2009. Knee muscle forces during walking and running in patellofemoral pain patients and painfree controls. Journal of Biomechanics 42, 898–905. Correa, T.A., Baker, R., Kerr Graham, H., Pandy, M.G., 2011. Accuracy of generic musculoskeletal models in predicting the functional roles of muscles in human gait. Journal of Biomechanics 44, 2096–2105. Correa, T.A., Crossley, K.M., Kim, H.J., Pandy, M.G., 2010. Contributions of individual muscles to hip joint contact force in normal walking. Journal of Biomechanics 43, 1618–1622. De Groote, F., Van Campen, A., Jonkers, I., De Schutter, J., 2010. Sensitivity of dynamic simulations of gait and dynamometer experiments to hill muscle model parameters of knee flexors and extensors. Journal of Biomechanics 43, 1876–1883. Erdemir, A., Mclean, S., Herzog, W., Van Den Bogert, A.J., 2007. Model-based estimation of muscle forces exerted during movements. Clinical Biomechanics 22, 131–154. Friederich, J.A., Brand, R.A., 1990. Muscle fiber architecture in the human lower limb. Journal of Biomechanics 23, 91–95. Herzog, W., 1992. Sensitivity of muscle force estimations to changes in muscle input parameters using nonlinear optimization approaches. Journal of Biomechanical Engineering 114, 267–268. Lin, Y.-C., Kim, H.-J., Pandy, M.G., 2011. A computationally efficient method for assessing muscle function during human locomotion. International Journal of Numerical Methods in Biomedical Engineering 27, 436–449. Liu, M.Q., Anderson, F.C., Pandy, M.G., Delp, S.L., 2006. Muscles that support the body also modulate forward progression during walking. Journal of Biomechanics 39, 2623–2630. Nagano, A., Komura, T., 2003. Longer moment arm results in smaller joint moment development, power and work outputs in fast motions. Journal of Biomechanics 36, 1675–1681. Nemeth, G., Ohlsen, H., 1985. In vivo moment arm lengths for hip extensor muscles at different angles of hip flexion. Journal of Biomechanics 18, 129–140. Pandy, M.G., Andriacchi, T.P., 2010. Muscle and joint function in human locomotion. Annual Review of Biomedical Engineering 12, 401–433. Pandy, M.G., Lin, Y.C., Kim, H.J., 2010. Muscle coordination of mediolateral balance in normal walking. Journal of Biomechanics 43, 2055–2064. Piazza, S.J., 2006. Muscle-driven forward dynamic simulations for the study of normal and pathological gait. Journal of NeuroEngineering and Rehabilitation 3, 5. Raikova, R.T., Prilutsky, B.I., 2001. Sensitivity of predicted muscle forces to parameters of the optimization-based human leg model revealed by analytical and numerical analyses. Journal of Biomechanics 34, 1243–1255. Redl, C., Gfoehler, M., Pandy, M.G., 2007. Sensitivity of muscle force estimates to variations in muscle–tendon properties. Human Movement Science 26, 306–319. Schache, A.G., Blanch, P.D., Brown, N.A., Pandy, M.G., 2012. Mechanics of the human hamstring muscles during sprinting. Medicine & Science in Sports & Exercise 44, 647–658. Scheys, L., Van Campenhout, A., Spaepen, A., Suetens, P., Jonkers, I., 2008. Personalized MR-based musculoskeletal models compared to rescaled generic models in the presence of increased femoral anteversion: effect on hip moment arm lengths. Gait & Posture 28, 358–365. Scovil, C.Y., Ronsky, J.L., 2006. Sensitivity of a Hill-based muscle model to perturbations in model parameters. Journal of Biomechanics 39, 2055–2063. Steele, K.M., Seth, A., Hicks, J.L., Schwartz, M.S., Delp, S.L., 2010. Muscle contributions to support and progression during single-limb stance in crouch gait. Journal of Biomechanics 43, 2099–2105. Tate, C.M., Williams, G.N., Barrance, P.J., Buchanan, T.S., 2006. Lower extremity muscle morphology in young athletes: an MRI-based analysis. Medicine & Science in Sports & Exercise 38, 122–128. Thelen, D.G., 2003. Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. Journal of Biomechanical Engineering 125, 70–77. Valero-Cuevas, F.J., Johanson, M.E., Towles, J.D., 2003. Towards a realistic biomechanical model of the thumb: the choice of kinematic description may be more critical than the solution method or the variability/uncertainty of musculoskeletal parameters. Journal of Biomechanics 36, 1019–1030. Viceconti, M.T.D., Taddei, F., Martelli, S., Clapworthy, G.J., Van Sint Jan, S., 2006. Biomechanics modeling of the musculoskeletal apparatus: status and key issues. Proceedings of the IEEE 94, 725–739. Ward, S.R., Eng, C.M., Smallwood, L.H., Lieber, R.L., 2009. Are current measurements of lower extremity muscle architecture accurate? Clinical Orthopaedics and Related 467, 1074–1082. Wilson, D.L., Zhu, Q., Duerk, J.L., Mansour, J.M., Kilgore, K., Crago, P.E., 1999. Estimation of tendon moment arms from three-dimensional magnetic resonance images. Annals of Biomedical Engineering 27, 247–256.

D.C. Ackland et al. / Journal of Biomechanics 45 (2012) 1463–1471

Wilson, N.A., Sheehan, F.T., 2009. Dynamic in vivo 3-dimensional moment arms of the individual quadriceps components. Journal of Biomechanics 42, 1891–1897. Xiao, M., Higginson, J., 2010. Sensitivity of estimated muscle force in forward simulation of normal walking. Journal of Applied Biomechanics 26, 142–149. Yu, J., Ackland, D.C., Pandy, M.G., 2011. Shoulder muscle function depends on elbow joint position: an illustration of dynamic coupling in the upper limb. Journal of Biomechanics 44, 1859–1868.

1471

Zajac, F.E., 1989. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biotechnology 17, 359–411. Zajac, F.E., Gordon, M.E., 1989. Determining muscle’s force and action in multi-articular movement. Exercise and Sport Sciences Reviews 17, 187–230.