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

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. 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. 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. 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. 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