0
Research Papers

Flexing Computational Muscle: Modeling and Simulation of Musculotendon Dynamics

[+] Author and Article Information
Matthew Millard

e-mail: mjhmilla@stanford.edu

Thomas Uchida

e-mail: tkuchida@stanford.edu

Ajay Seth

e-mail: aseth@stanford.edu
Department of Bioengineering,
Stanford University,
Stanford, CA 94305

Scott L. Delp

Department of Bioengineering,
Department of Mechanical Engineering,
Stanford University,
Stanford, CA 94305
e-mail: delp@stanford.edu

Data available online at simtk.org.

Data available online at simtk.org.

1Corresponding author.

Contributed by the Bioengineering Division of ASME for publication in the JOURNAL OF BIOMECHANICAL ENGINEERING. Manuscript received November 16, 2012; final manuscript received January 7, 2013; accepted manuscript posted January 18, 2013; published online February 7, 2013. Editor: Victor H. Barocas.

J Biomech Eng 135(2), 021005 (Feb 07, 2013) (11 pages) Paper No: BIO-12-1561; doi: 10.1115/1.4023390 History: Received November 16, 2012; Revised January 07, 2013; Accepted January 18, 2013

Muscle-driven simulations of human and animal motion are widely used to complement physical experiments for studying movement dynamics. Musculotendon models are an essential component of muscle-driven simulations, yet neither the computational speed nor the biological accuracy of the simulated forces has been adequately evaluated. Here we compare the speed and accuracy of three musculotendon models: two with an elastic tendon (an equilibrium model and a damped equilibrium model) and one with a rigid tendon. Our simulation benchmarks demonstrate that the equilibrium and damped equilibrium models produce similar force profiles but have different computational speeds. At low activation, the damped equilibrium model is 29 times faster than the equilibrium model when using an explicit integrator and 3 times faster when using an implicit integrator; at high activation, the two models have similar simulation speeds. In the special case of simulating a muscle with a short tendon, the rigid-tendon model produces forces that match those generated by the elastic-tendon models, but simulates 2–54 times faster when an explicit integrator is used and 6–31 times faster when an implicit integrator is used. The equilibrium, damped equilibrium, and rigid-tendon models reproduce forces generated by maximally-activated biological muscle with mean absolute errors less than 8.9%, 8.9%, and 20.9% of the maximum isometric muscle force, respectively. When compared to forces generated by submaximally-activated biological muscle, the forces produced by the equilibrium, damped equilibrium, and rigid-tendon models have mean absolute errors less than 16.2%, 16.4%, and 18.5%, respectively. To encourage further development of musculotendon models, we provide implementations of each of these models in OpenSim version 3.1 and benchmark data online, enabling others to reproduce our results and test their models of musculotendon dynamics.

FIGURES IN THIS ARTICLE
<>
Copyright © 2013 by ASME
Your Session has timed out. Please sign back in to continue.

References

Zajac, F. E., Neptune, R. R., and Kautz, S. A., 2002, “Biomechanics and Muscle Coordination of Human Walking: Part I: Introduction to Concepts, Power Transfer, Dynamics and Simulations,” Gait Posture, 16(3), pp. 215–232. [CrossRef] [PubMed]
Zajac, F. E., Neptune, R. R., and Kautz, S. A., 2003, “Biomechanics and Muscle Coordination of Human Walking: Part II: Lessons From Dynamical Simulations and Clinical Implications,” Gait Posture, 17(1), pp. 1–17. [CrossRef] [PubMed]
Anderson, F. C., and Pandy, M. G., 2001, “Dynamic Optimization of Human Walking,” ASME J. Biomech. Eng., 123(5), pp. 381–390. [CrossRef]
Ackermann, M., and van den Bogert, A. J., 2010, “Optimality Principles for Model-Based Prediction of Human Gait,” J. Biomech., 43(6), pp. 1055–1060. [CrossRef] [PubMed]
Arnold, E. M., and Delp, S. L., 2011, “Fibre Operating Lengths of Human Lower Limb Muscles During Walking,” Philos. T. R. Soc. B, 366(1570), pp. 1530–1539. [CrossRef]
Liu, M. Q., Anderson, F. C., Pandy, M. G., and Delp, S. L., 2006, “Muscles That Support the Body Also Modulate Forward Progression During Walking,” J. Biomech., 39(14), pp. 2623–2630. [CrossRef] [PubMed]
Hamner, S. R., Seth, A., and Delp, S. L., 2010, “Muscle Contributions to Propulsion and Support During Running,” J. Biomech., 43(14), pp. 2709–2716. [CrossRef] [PubMed]
Neptune, R. R., and Sasaki, K., 2005, “Ankle Plantar Flexor Force Production Is an Important Determinant of the Preferred Walk-to-Run Transition Speed,” J. Exp. Biol., 208(5), pp. 799–808. [CrossRef] [PubMed]
Selbie, W. S., and Caldwell, G. E., 1996, “A Simulation Study of Vertical Jumping From Different Starting Postures,” J. Biomech., 29(9), pp. 1137–1146. [CrossRef] [PubMed]
Neptune, R. R., and Hull, M. L., 1999, “A Theoretical Analysis of Preferred Pedaling Rate Selection in Endurance Cycling,” J. Biomech., 32(4), pp. 409–415. [CrossRef] [PubMed]
Neptune, R. R., and van den Bogert, A. J., 1998, “Standard Mechanical Energy Analyses Do Not Correlate With Muscle Work in Cycling,” J. Biomech., 31(3), pp. 239–245. [CrossRef] [PubMed]
van der Krogt, M. M., Delp, S. L., and Schwartz, M. H., 2012, “How Robust Is Human Gait to Muscle Weakness?,” Gait Posture, 36(1), pp. 113–119. [CrossRef] [PubMed]
Steele, K. M., Seth, A., Hicks, J. L., Schwartz, M. S., and Delp, S. L., 2010, “Muscle Contributions to Support and Progression During Single-Limb Stance in Crouch Gait,” J. Biomech., 43(11), pp. 2099–2105. [CrossRef] [PubMed]
Crabtree, C. A., and Higginson, J. S., 2009, “Modeling Neuromuscular Effects of Ankle Foot Orthoses (AFOs) in Computer Simulations of Gait,” Gait Posture, 29(1), pp. 65–70. [CrossRef] [PubMed]
Hicks, J. L., Schwartz, M. H., Arnold, A. S., and Delp, S. L., 2008, “Crouched Postures Reduce the Capacity of Muscles to Extend the Hip and Knee During the Single-Limb Stance Phase of Gait,” J. Biomech., 41(5), pp. 960–967. [CrossRef] [PubMed]
Fregly, B. J., Reinbolt, J. A., Rooney, K. L., Mitchell, K. H., and Chmielewski, T. L., 2007, “Design of Patient-Specific Gait Modifications for Knee Osteoarthritis Rehabilitation,” IEEE Trans. Biomed. Eng., 54(9), pp. 1687–1695. [CrossRef] [PubMed]
Riener, R., and Fuhr, T., 1998, “Patient-Driven Control of FES-Supported Standing Up: A Simulation Study,” IEEE Trans. Rehab. Eng., 6(2), pp. 113–124. [CrossRef]
Delp, S. L., Loan, J. P., Hoy, M. G., Zajac, F. E., Topp, E. L., and Rosen, J. M., 1990, “An Interactive Graphics-Based Model of the Lower Extremity to Study Orthopaedic Surgical Procedures,” IEEE Trans. Biomed. Eng., 37(8), pp. 757–767. [CrossRef] [PubMed]
Rasmussen, J., Tørholm, S., and de Zee, M., 2009, “Computational Analysis of the Influence of Seat Pan Inclination and Friction on Muscle Activity and Spinal Joint Forces,” Int. J. Ind. Ergonom., 39(1), pp. 52–57. [CrossRef]
Eisenberg, E., Hill, T. L., and Chen, Y., 1980, “Cross-Bridge Model of Muscle Contraction. Quantitative Analysis,” Biophys. J., 29(2), pp. 195–227. [CrossRef] [PubMed]
Zahalak, G. I., and Ma, S.-P., 1990, “Muscle Activation and Contraction: Constitutive Relations Based Directly on Cross-Bridge Kinetics,” ASME J. Biomech. Eng., 112(1), pp. 52–62. [CrossRef]
Haselgrove, J. C., and Huxley, H. E., 1973, “X-ray Evidence for Radial Cross-Bridge Movement and for the Sliding Filament Model in Actively Contracting Skeletal Muscle,” J. Mol. Biol., 77(4), pp. 549–568. [CrossRef] [PubMed]
Zajac, F. E., 1989, “Muscle and Tendon: Properties, Models, Scaling, and Application to Biomechanics and Motor Control,” Crit. Rev. Biomed. Eng., 17(4), pp. 359–411. Available at: http://europepmc.org/abstract/MED/2676342
Epstein, M., and Herzog, W., 1998, Theoretical Models of Skeletal Muscle: Biological and Mathematical Considerations, Wiley, New York.
Winters, J. M., and Stark, L., 1987, “Muscle Models: What Is Gained and What Is Lost by Varying Model Complexity,” Biol. Cybern., 55(6), pp. 403–420. [CrossRef] [PubMed]
Shi, P., and McPhee, J., 2000, “Dynamics of Flexible Multibody Systems Using Virtual Work and Linear Graph Theory,” Multibody Syst. Dyn., 4(4), pp. 355–381. [CrossRef]
Stoianovici, D., and Hurmuzlu, Y., 1996, “A Critical Study of the Applicability of Rigid-Body Collision Theory,” ASME J. Appl. Mech., 63(2), pp. 307–316. [CrossRef]
Bowden, F. P., and Leben, L., 1939, “The Nature of Sliding and the Analysis of Friction,” Proc. R. Soc. Lon. Ser. A, 169(938), pp. 371–391. [CrossRef]
Krylow, A. M., and Sandercock, T. G., 1997, “Dynamic Force Responses of Muscle Involving Eccentric Contraction,” J. Biomech., 30(1), pp. 27–33. [CrossRef] [PubMed]
Perreault, E. J., Heckman, C. J., and Sandercock, T. G., 2003, “Hill Muscle Model Errors During Movement Are Greatest Within the Physiologically Relevant Range of Motor Unit Firing Rates,” J. Biomech., 36(2), pp. 211–218. [CrossRef] [PubMed]
Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., Guendelman, E., and Thelen, D. G., 2007, “OpenSim: Open-Source Software to Create and Analyze Dynamic Simulations of Movement,” IEEE Trans. Biomed. Eng., 54(11), pp. 1940–1950. [CrossRef] [PubMed]
Seth, A., Sherman, M., Reinbolt, J. A., and Delp, S. L., 2011, “OpenSim: A Musculoskeletal Modeling and Simulation Framework for In Silico Investigations and Exchange,” Procedia IUTAM, Symposium on Human Body Dynamics, 2, pp. 212–232. [CrossRef]
Reinbolt, J. A., Seth, A., and Delp, S. L., 2011, “Simulation of Human Movement: Applications Using OpenSim,” Procedia IUTAM, Symposium on Human Body Dynamics, 2, pp. 186–198. [CrossRef]
Matsubara, I., and Elliott, G. F., 1972, “X-ray Diffraction Studies on Skinned Single Fibres of Frog Skeletal Muscle,” J. Mol. Biol., 72(3), pp. 657–669. [CrossRef] [PubMed]
Randhawa, A., Jackman, M. E., and Wakeling, J. M., 2012, “Muscle Gearing During Isotonic and Isokinetic Movements in the Ankle Plantarflexors,” Eur. J. Appl. Physiol., 113(2), pp. 437–447. [CrossRef] [PubMed]
Brainerd, E. L., and Azizi, E., 2005, “Muscle Fiber Angle, Segment Bulging and Architectural Gear Ratio in Segmented Musculature,” J. Exp. Biol., 208(17), pp. 3249–3261. [CrossRef] [PubMed]
Arnold, E. M., Ward, S. R., Lieber, R. L., and Delp, S. L., 2010, “A Model of the Lower Limb for Analysis of Human Movement,” Ann. Biomed. Eng., 38(2), pp. 269–279. [CrossRef] [PubMed]
Winters, J. M., 1995, “An Improved Muscle-Reflex Actuator for Use in Large-Scale Neuromusculoskeletal Models,” Ann. Biomed. Eng., 23(4), pp. 359–374. [CrossRef] [PubMed]
Thelen, D. G., 2003, “Adjustment of Muscle Mechanics Model Parameters to Simulate Dynamic Contractions in Older Adults,” ASME J. Biomech. Eng., 125(1), pp. 70–77. [CrossRef]
Cheng, E. J., Brown, I. E., and Loeb, G. E., 2000, “Virtual Muscle: A Computational Approach to Understanding the Effects of Muscle Properties on Motor Control,” J. Neurosci. Meth., 101(2), pp. 117–130. [CrossRef]
Magnusson, S. P., Aagaard, P., Rosager, S., Dyhre-Poulsen, P., and Kjaer, M., 2001, “Load–Displacement Properties of the Human Triceps Surae Aponeurosis In Vivo,” J. Physiol., 531(1), pp. 277–288. [CrossRef] [PubMed]
Maganaris, C. N., and Paul, J. P., 2002, “Tensile Properties of the In Vivo Human Gastrocnemius Tendon,” J. Biomech., 35(12), pp. 1639–1646. [CrossRef] [PubMed]
Winters, T. M., Takahashi, M., Lieber, R. L., and Ward, S. R., 2011, “Whole Muscle Length-Tension Relationships Are Accurately Modeled as Scaled Sarcomeres in Rabbit Hindlimb Muscles,” J. Biomech., 44(1), pp. 109–115. [CrossRef] [PubMed]
Gollapudi, S. K., and Lin, D. C., 2009, “Experimental Determination of Sarcomere Force–Length Relationship in Type-I Human Skeletal Muscle Fibers,” J. Biomech., 42(13), pp. 2011–2016. [CrossRef] [PubMed]
Mashima, H., 1984, “Force-Velocity Relation and Contractility in Striated Muscles,” Jap. J. Physiol., 34(1), pp. 1–17. [CrossRef]
Joyce, G. C., Rack, P. M. H., and Westbury, D. R., 1969, “The Mechanical Properties of Cat Soleus Muscle During Controlled Lengthening and Shortening Movements,” J. Physiol., 204(2), pp. 461–474. Available at: http://jp.physoc.org/content/204/2/461.abstract
Mortenson, M. E., 2006, Geometric Modeling, 3rd ed. Industrial Press, New York.
van den Bogert, A. J., Blana, D., and Heinrich, D., 2011, “Implicit Methods for Efficient Musculoskeletal Simulation and Optimal Control,” Procedia IUTAM, Symposium on Human Body Dynamics, 2, pp. 297–316. [CrossRef]
Hogan, N., 1985, “The Mechanics of Multi-Joint Posture and Movement Control,” Biol. Cybern., 52(5), pp. 315–331. [CrossRef] [PubMed]
Vinnars, E., Bergstöm, J., and Fürst, P., 1975, “Influence of the Postoperative State on the Intracellular Free Amino Acids in Human Muscle Tissue,” Ann. Surg., 182(6), pp. 665–671. [CrossRef] [PubMed]
Hairer, E., Nørsett, S. P., and Wanner, G., 1987, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed. Springer-Verlag, Berlin.
Sherman, M. A., Seth, A., and Delp, S. L., 2011, “Simbody: Multibody Dynamics for Biomedical Research,” Procedia IUTAM, Symposium on Human Body Dynamics, 2, pp. 241–261. [CrossRef]
Talmadge, R. J., Roy, R. R., Caiozzo, V. J., and Edgerton, V. R., 2002, “Mechanical Properties of Rat Soleus After Long-Term Spinal Cord Transection,” J. Appl. Physiol., 93(4), pp. 1487–1497. [CrossRef] [PubMed]
Nelder, J. A., and Mead, R., 1965, “A Simplex Method for Function Minimization,” Comput. J., 7(4), pp. 308–313. [CrossRef]
Scott, S. H., Brown, I. E., and Loeb, G. E., 1996, “Mechanics of Feline Soleus: I. Effect of Fascicle Length and Velocity on Force Output,” J. Muscle Res. Cell. M., 17(2), pp. 207–219. [CrossRef]
He, J., Levine, W. S., and Loeb, G. E., 1991, “Feedback Gains for Correcting Small Perturbations to Standing Posture,” IEEE Trans. Auto. Contr., 36(3), pp. 322–332. [CrossRef]
Hairer, E., and Wanner, G., 1991, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed. Springer-Verlag, Berlin.
Yamaguchi, G. T., and Zajac, F. E., 1989, “A Planar Model of the Knee Joint to Characterize the Knee Extensor Mechanism,” J. Biomech., 22(1), pp. 1–10. [CrossRef] [PubMed]
Anderson, F. C., and Pandy, M. G., 1999, “A Dynamic Optimization Solution for Vertical Jumping in Three Dimensions,” Comp. Meth. Biomech. Biomed. Eng., 2(3), pp. 201–231. [CrossRef]
Herzog, W., and Leonard, T. R., 2002, “Force Enhancement Following Stretching of Skeletal Muscle: A New Mechanism,” J. Exp. Biol., 205(9), pp. 1275–1283. Available at: http://jeb.biologists.org/content/205/9/1275.short
Ranatunga, K. W., 1982, “Temperature-Dependence of Shortening Velocity and Rate of Isometric Tension Development in Rat Skeletal Muscle,” J. Physiol., 329, pp. 465–483. Available at: http://jp.physoc.org/content/329/1/465.short
Westerblad, H., Allen, D. G., Bruton, J. D., Andrade, F. H., and Lännergren, J., 1998, “Mechanisms Underlying the Reduction of Isometric Force in Skeletal Muscle Fatigue,” Acta Physiol. Scand., 162(3), pp. 253–260. [CrossRef] [PubMed]

Figures

Grahic Jump Location
Fig. 1

Simplified geometric representation of muscle fibers and tendon for musculotendon modeling. Muscle fibers are assumed to be straight, parallel, of equal length, coplanar, and attached to tendon at a pennation angle (α). As the muscle shortens, the distance h remains constant and the pennation angle increases. Adapted from Zajac [23].

Grahic Jump Location
Fig. 2

Muscle-driven simulations use a model of musculotendon contraction dynamics to determine muscle lengths (ℓM), velocities (vM), and forces (fM) from neural excitations (u), generalized coordinates (q), and generalized speeds (q·). A model of activation dynamics determines muscle activations (a) from neural excitations (u). A musculoskeletal model determines musculotendon lengths and velocities (ℓMT and vMT) from the generalized coordinates and speeds (q and q·). The model of musculotendon contraction dynamics uses the results of the activation dynamics and the musculoskeletal model to produce a forward simulation of muscle length (ℓM), velocity (vM), and force (fM).

Grahic Jump Location
Fig. 3

Schematic of the equilibrium musculotendon model (a), tendon-force–length curve (b), active- and passive-force–length curves (c), and force–velocity curve (d). Experimental data for the tendon-force–length curve in panel (b) are illustrated as 95% confidence intervals [41,42]. Data points in panels (c) and (d) denote experimental data for the force–length [43,44] and force–velocity [45,46] curves. The default curves used in the musculotendon models are shown in comparison to these experimental data. Note that f˜M and f˜T represent, respectively, muscle force and tendon force normalized by foM.

Grahic Jump Location
Fig. 4

Schematic of the musculotendon actuator used in the elastic-tendon and rigid-tendon computational benchmarks. In the simulations, the initial length of the musculotendon actuator was set to c1 = ℓsT+ℓoMcosαo and then the length was varied by c2sin2πt.

Grahic Jump Location
Fig. 5

Normalized tendon force profiles of the equilibrium (left) and damped equilibrium (right) musculotendon models. Each muscle underwent constant-activation, sinusoidal-displacement simulations of 1-s duration with activations varying from 0 (0.01 for the equilibrium model) to 1 in increments of 0.1. The average absolute difference between normalized force profiles generated by the models was less than 0.3% at low activation and less than 2.5% at high activation.

Grahic Jump Location
Fig. 6

Time required to generate the simulation results shown in Fig. 5 using explicit (left) and implicit (right) integrators. Wall clock time was divided by the amount of time simulated to obtain the real-time fraction; values below 1 indicate that the simulations completed faster than real time.

Grahic Jump Location
Fig. 7

Differences in force profiles generated by rigid-tendon and damped equilibrium musculotendon models as a function of tendon slack length (ℓsT) normalized by the optimal fiber length (ℓoM). The force profiles were obtained using maximal-activation, sinusoidal-displacement simulations. The mean (solid curve) and standard deviation (shaded area) of the absolute difference in force profiles over the duration of each simulation show that errors in the rigid-tendon model increase rapidly as normalized tendon slack length increases.

Grahic Jump Location
Fig. 8

Computational speed of the damped equilibrium (black curve) and rigid-tendon (gray curve) musculotendon models as functions of tendon slack length (ℓsT) normalized by the optimal fiber length (ℓoM) using explicit (left) and implicit (right) integrators. Wall clock time was divided by the amount of time simulated to obtain the real-time fraction; values below 1 indicate that the simulations completed faster than real time.

Grahic Jump Location
Fig. 9

Waveform of the change in musculotendon length [29] used in the maximal-activation biological benchmark, shown here with a maximum length change of 1.0 mm. This waveform was scaled by the appropriate maximum length change, added to ℓsT+ℓoMcosαo-2 mm, and then used to prescribe the displacement of the free end of the tendon.

Grahic Jump Location
Fig. 10

Comparison of experimental [29] (gray curve) and simulated (black curve) force profiles from maximally-activated muscle undergoing length changes of maximum amplitude 0.05 (top), 0.1, 0.25, 0.5, 1.0, and 2.0 mm (bottom). The maximum isometric force used in the model was foM = 1.27 N and the optimal fiber length was ℓoM = 17.1 mm.

Grahic Jump Location
Fig. 11

Waveform of the change in musculotendon length [30] used in the submaximal-activation biological benchmark, shown here with a maximum length change of 1.0 mm. This waveform was scaled by the appropriate maximum length change (either 1.0 mm or 8.0 mm), added to ℓsT+ℓoMcosαo-4 mm, and then used to prescribe the displacement of the free end of the tendon.

Grahic Jump Location
Fig. 12

Comparison of experimental [30] (gray curve) and simulated (black curve) force profiles from submaximally-activated muscle undergoing length changes of maximum amplitude 1.0 mm (left) and 8.0 mm (right). Constant-frequency stimulation rates of 10 (top), 20, and 30 Hz (bottom) were applied to the biological muscle. Force profiles measured experimentally during isometric trials were used to determine the corresponding activation signals that must be applied to the damped equilibrium musculotendon model. The maximum isometric force used in the model was foM = 25.1 N and the optimal fiber length was ℓoM = 30 mm.

Grahic Jump Location
Fig. 13

Comparison of experimental [30] (gray curve) and simulated (black curve) force profiles from submaximally-activated muscle undergoing length changes of maximum amplitude 1.0 mm (left) and 8.0 mm (right). Random stimulation signals were applied to the biological muscle, with mean frequencies of 10 (top), 20, and 30 Hz (bottom), as described by Perreault et al. [30]. Force profiles measured experimentally during isometric trials were used to determine the corresponding activation signals that must be applied to the damped equilibrium musculotendon model. The maximum isometric force used in the model was foM = 25.1 N and the optimal fiber length was ℓoM = 30 mm.

Tables

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In