Pulling Out the Stops: Designing a Muscle for a Tug-of-War Competition

Objective

The force-producing properties of muscle are complex, highly nonlinear, and can have substantial effects on movement (see McMahon, 1984 for a review). For simplicity, lumped-parameter, dimensionless muscle models that are capable of representing a variety of muscles with different architectures are commonly used in the dynamic simulation of movement (Zajac, 1989). In this exercise, you will explore the differential equations that describe muscle activation and muscle–tendon contraction dynamics when using a Hill-type muscle model. You will use OpenSim to implement a simple muscle–tendon model and conduct simulations to investigate how various model parameters affect the dynamic response of the actuator. The lab will conclude with a Virtual Muscle Tug-of-War in which you will design an optimal muscle and compete against other user-generated muscles. May the best muscle win!

By working through this lab, you will:

• Become familiar with the differential equations describing muscle activation and muscle–tendon contraction dynamics.
• Learn how to model and simulate muscle–tendon dynamics using OpenSim.
• Become comfortable with modifying existing code that models muscle activation and muscle–tendon mechanics.
• Explore the effect of various model parameters and simulation conditions on the dynamic response of muscle.
• Design your own optimal muscle–tendon actuator to compete in a virtual muscle tug-of-war.

Acknowledgements

The original lab was designed by Jeff Reinbolt, B.J. Fregly, Kate Saul, Darryl Thelen, Silvia Blemker, Clay Anderson, and Scott Delp. The lab was refined by Hoa Hoang and Daniel Jacobs.

Model

The model in this exercise consists of a cube with a single translational degree of freedom along the Z-axis. A Tug_of_War model has been included in the OpenSim distribution (Models/Tug_of_War/) that uses two Thelen 2003 Muscles to move the cube. In this exercise, we will use two Millard 2012 Muscles instead of the Thelen muscles. The muscles are arranged to pull on opposite sides of the cube. The cube has a mass of 20 kg and sides of length 0.1 m, and the distance between the fixed ground supports is 0.7 m. Thus, each muscle–tendon actuator is 0.3 m long when the block is centered.

Background

Millard 2012 Muscle Model

Musculotendon actuators consist of an active contractile element, a passive elastic element, and an elastic tendon. The maximum active force a muscle can develop varies nonlinearly with its length, represented by the active-force–length curve $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big)\end{array} //]]>$, peaking at a force of $//<![CDATA[ \begin{array}{l}f_o^{\,\mathrm{M}}\end{array} //]]>$ at a length of $//<![CDATA[ \begin{array}{l}\ell_o^{\,\mathrm{M}}\end{array} //]]>$ (the tilde is used to denote forces, velocities, muscle lengths, and tendon lengths that are normalized by $//<![CDATA[ \begin{array}{l}f_o^{\,\mathrm{M}}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}v_{\mathrm{max}}^{\,\mathrm{M}}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\ell_o^{\,\mathrm{M}}\end{array} //]]>$, and $//<![CDATA[ \begin{array}{l}\ell_\mathrm{s}^{\mathrm{T}}\end{array} //]]>$, respectively). During non-isometric contractions, the force developed by muscle varies nonlinearly with its rate of lengthening, which is represented by the force–velocity curve $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{V}} \big( \tilde{v}^{\,\mathrm{M}} \big)\end{array} //]]>$. Force is also developed when muscle is stretched beyond a threshold length, regardless of whether the muscle is activated, which is represented by the passive-force–length curve $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{PE}} \big( \tilde{l}^{\,\mathrm{M}} \big)\end{array} //]]>$. Muscle force ($//<![CDATA[ \begin{array}{l}f^{\,\mathrm{M}}\end{array} //]]>$) is computed using these curves as follows:

$$//<![CDATA[ \begin{array}{l}\displaystyle f^{\,\mathrm{M}} = f_o^{\,\mathrm{M}} \left( a \mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) \mathbf{f}^{\mathrm{V}} \big( \tilde{v}^{\,\mathrm{M}} \big) + \mathbf{f}^{\mathrm{PE}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) \right)\end{array} //]]>$$

where $//<![CDATA[ \begin{array}{l}a\end{array} //]]>$ is the muscle activation, which ranges from $//<![CDATA[ \begin{array}{l}a_{\mathrm{min}}\end{array} //]]>$ to 1.

Muscle attaches to bone through tendon. Since a long tendon may stretch appreciably beyond its slack length ($//<![CDATA[ \begin{array}{l}\ell_{\mathrm{s}}^{\mathrm{T}}\end{array} //]]>$) when under tension, tendon is modeled as a nonlinear elastic element developing force according to the tendon-force–length curve $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{T}} \big( \tilde{\ell}^{\,\mathrm{T}} \big)\end{array} //]]>$. Muscle fibers attach to tendon at a pennation angle ($//<![CDATA[ \begin{array}{l}\alpha\end{array} //]]>$), scaling the force they transmit to the tendon. If the tendon is assumed to be elastic and the mass of the muscle is assumed to be negligible, then the muscle–tendon force equilibrium equation (i.e., $//<![CDATA[ \begin{array}{l}f^{\,\mathrm{M}} \cos \alpha - f^{\,T} = 0\end{array} //]]>$) takes the following form:

$$//<![CDATA[ \begin{array}{l}\displaystyle f_o^{\,\mathrm{M}} \left( a \mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) \mathbf{f}^{\mathrm{V}} \big( \tilde{v}^{\,\mathrm{M}} \big) + \mathbf{f}^{\mathrm{PE}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) \right) \cos \alpha - f_o^{\,\mathrm{M}} \mathbf{f}^{\mathrm{T}} \big( \tilde{\ell}^{\,\mathrm{T}} \big) = 0\qquad(1)\end{array} //]]>$$

We have developed default force curves for the musculotendon model that have been fit to experimental data (see Millard et al., 2013 for details). These curves can be adjusted to model muscle and tendon whose characteristics deviate from these default patterns.

In a forward-dynamic simulation, the force generated by a musculotendon actuator is calculated from the length ($//<![CDATA[ \begin{array}{l}\ell^{\mathrm{M}}\end{array} //]]>$), velocity ($//<![CDATA[ \begin{array}{l}v^{\mathrm{M}}\end{array} //]]>$), and activation ($//<![CDATA[ \begin{array}{l}a\end{array} //]]>$) of the muscle. The length $//<![CDATA[ \begin{array}{l}\ell^{\mathrm{M}}\end{array} //]]>$ is obtained by integrating the velocity $//<![CDATA[ \begin{array}{l}v^{\mathrm{M}}\end{array} //]]>$, which is obtained by solving Eq. (1) for the normalized muscle velocity ($//<![CDATA[ \begin{array}{l}\tilde{v}^{\,\mathrm{M}}\end{array} //]]>$), thereby providing an ordinary differential equation (see Zajac, 1989):

$$//<![CDATA[ \begin{array}{l}\displaystyle \tilde{v}^{\,\mathrm{M}} = \mathbf{f}_{\mathrm{inv}}^{\mathrm{V}} \left( \frac{ \mathbf{f}^{\mathrm{T}} \big( \tilde{\ell}^{\,\mathrm{T}} \big) / \cos \alpha \ -\ \mathbf{f}^{\mathrm{PE}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) }{ a \mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) } \right)\qquad(2)\end{array} //]]>$$

where $//<![CDATA[ \begin{array}{l}\mathbf{f}_{\mathrm{inv}}^{\mathrm{V}}\end{array} //]]>$ is the inverse of the force–velocity curve. Although Eq. (1) is devoid of numerical singularities, Eq. (2) has four: as $//<![CDATA[ \begin{array}{l}\alpha \rightarrow 90^{\circ}\end{array} //]]>$, as $//<![CDATA[ \begin{array}{l}a \rightarrow 0\end{array} //]]>$, as $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) \rightarrow 0\end{array} //]]>$, and as $//<![CDATA[ \begin{array}{l}\partial \mathbf{f}^{\mathrm{V}} \big( \tilde{v}^{\,\mathrm{M}} \big) / \partial \tilde{v}^{\,\mathrm{M}} \rightarrow 0\end{array} //]]>$. Since these conditions are often encountered during a simulation, the quantities causing singularities in Eq. (2) are altered so that the singularities are approached but never reached: $//<![CDATA[ \begin{array}{l}\alpha < 90^{\circ}\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}a > 0\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\mathbf{f}^{\mathrm{L}} \big( \tilde{\ell}^{\,\mathrm{M}} \big) > 0\end{array} //]]>$, and $//<![CDATA[ \begin{array}{l}\partial \mathbf{f}^{\mathrm{V}} \big( \tilde{v}^{\,\mathrm{M}} \big) / \partial \tilde{v}^{\,\mathrm{M}} > 0\end{array} //]]>$.

Without modifying the formulation of the equilibrium model, the muscle is able to reach unrealistically short lengths (Thelen, 2003) and cannot be simulated when fully deactivated. We use a unilateral constraint on muscle length to prevent the muscle from becoming unrealistically short:

$$//<![CDATA[ \begin{align} \tilde{v}^{\,\mathrm{M}} = \begin{cases} 0 &: \tilde{\ell}^{\,\mathrm{M}} \leq \tilde{\ell}_{\mathrm{min}}^{\,\mathrm{M}} \;\;\text{and}\;\; \tilde{v}^{\,\mathrm{M*}} < 0\\ \tilde{v}^{\,\mathrm{M*}} &: \text{otherwise} \end{cases} \end{align} //]]>$$

where $//<![CDATA[ \begin{array}{l}\tilde{v}^{\,\mathrm{M*}}\end{array} //]]>$ is a candidate value for $//<![CDATA[ \begin{array}{l}\tilde{v}^{\,\mathrm{M}}\end{array} //]]>$ computed using Eq. (2).

Activation dynamics

A muscle can neither generate force nor relax instantaneously. The development of force is a complex sequence of events that begins with the firing of motor units and culminates in the formation of actin–myosin cross-bridges within the myofibrils of the muscle. When the motor units of a muscle depolarize, action potentials are elicited in the fibers of the muscle and cause calcium ions to be released from the sarcoplasmic reticulum. The increase in calcium ion concentrations then initiates the cross-bridge formation between the actin and myosin filaments. In isolated muscle twitch experiments, the delay between a motor unit action potential and the development of peak force has been observed to vary from as little as 5 milliseconds for fast ocular muscles to as much as 40 or 50 milliseconds for muscles comprised of higher percentages of slow-twitch fibers. The relaxation of muscle depends on the re-uptake of calcium ions into the sarcoplasmic reticulum. This re-uptake is a slower process than the calcium ion release, so the time required for muscle force to fall can be considerably longer than the time for it to develop.

The activation dynamics of muscle can be modeled with a first-order differential equation. This equation relates the rate of change of muscle activation (i.e., the concentration of calcium ions within the muscle) to the muscle excitation (i.e., the firing of motor units):

$$//<![CDATA[ \begin{array}{l}\displaystyle \frac{da}{dt} = \frac{u-a}{\tau(a,u)}\qquad(1)\end{array} //]]>$$

where $//<![CDATA[ \begin{array}{l}u\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}a\end{array} //]]>$ are the excitation and activation signals, respectively. In the model, activation is allowed to vary continuously between 0 (no contraction) and 1 (full contraction). In the body, the activation of a muscle is a function of the number of motor units recruited and the firing frequency of these motor units. Some models of excitation–contraction coupling distinguish these two control mechanisms (Hatze, 1976), but it is often not computationally feasible to use such models when conducting complex dynamic simulations. In a simulation, the muscle excitation signal is assumed to represent the net effect of both motor neuron recruitment and firing frequency. Like muscle activation, the excitation signal is also allowed to vary continuously between 0 (no excitation) and 1 (full excitation). The activation and deactivation time constants can be assumed to be 10 and 40 ms, respectively (Zajac, 1989; Winters, 1990).

The activation model presented by Thelen (2003) closely follows the activation dynamic model found in Winters (1995 – Eq. 2, line 2 and Eq. 3), where the time derivative of activation ($//<![CDATA[ \begin{array}{l}da/dt\end{array} //]]>$) is equal to the difference between excitation ($//<![CDATA[ \begin{array}{l}u\end{array} //]]>$) and activation ($//<![CDATA[ \begin{array}{l}a\end{array} //]]>$) scaled by a variable time constant ($//<![CDATA[ \begin{array}{l}\tau(a,u)\end{array} //]]>$). The primary difference between the activation model presented by Thelen and the one presented by Winters lies in their expressions for $//<![CDATA[ \begin{array}{l}\tau(a,u)\end{array} //]]>$, which differ only in the values of the two coefficients shown in parentheses in Eqs. 2 and 3:

$$//<![CDATA[ \begin{array}{l}\tau(a,u) = \begin{cases} t_{\mathrm{act}} \left( 0.5+1.5a \right) &: u > a\qquad(2)\\ t_{\mathrm{deact}}\ / \left( 0.5+1.5a \right) &: u \leq a\qquad(3)\end{cases}\end{array} //]]>$$

Note that the model of activation dynamics presented in Eq. 1 does not respect a lower bound for activation. Equilibrium muscle models (commonly used to model muscle in lumped-parameter musculoskeletal simulations) have a singularity in their state equations when activation is zero, making the above activation dynamic model unsuitable for simulations using equilibrium muscle models. Equation 1 can be made to respect a lower bound on activation by first constraining both activation and excitation to remain between the minimum activation level (0.01 by default) and the maximum activation level (1.0 by default).

 

Design Challenge

The competition is a single-elimination tournament between pairs of muscles. Each match is a one-second forward dynamic simulation where each muscle starts with minimal activation (a = 0.01). The match is won by the muscle that has moved the block to its side of the arena at the end of the simulation. Your challenge is to specify the muscle–tendon parameters and excitation time history necessary to overcome the other challengers.

References

1. Anderson, F.C. and Pandy, M.G. (1999). A dynamic optimization solution for vertical jumping in three dimensions. Computer Methods in Biomechanics and Biomedical Engineering, 2(3):201–231.
2. Hatze, H. (1976). The complete optimization of a human motion. Mathematical Biosciences, 28(1–2):99–135.
3. McMahon, T.A. (1984). Muscles, Reflexes, and Locomotion. Princeton University Press, Princeton, New Jersey.
4. Millard, M., Uchida, T., Seth, A., Delp, S.L. (2013). Flexing computational muscle: modeling and simulation of musculotendon dynamics. ASME Journal of Biomechanical Engineering, 135(2):021005.
5. Schutte, L.M. (1993). Using Musculoskeletal Models to Explore Strategies for Improving Performance in Electrical Stimulation-Induced Leg Cycle Ergometry. PhD Dissertation, Mechanical Engineering Department, Stanford University.
6. Thelen, D.G. (2003). Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. ASME Journal of Biomechanical Engineering, 125(1):70–77.
7. Winters, J.M. (1990). Hill-based muscle models: a systems engineering perspective, in Multiple Muscle Systems: Biomechanics and Movement Organization, edited by Winters, J.M. and Woo, S.L., Springer-Verlag, New York.
8. Zajac, F.E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Enginering, 17(4):359–411.