Millard 2012 Muscle Models

Introduction

This page describes the Hill-type muscle models implemented in OpenSim as OpenSim::Millard2012EquilibriumMuscle by Dr. Matthew Millard and described in detail in Millard et al. (2013). Benchmarking data used for testing this model can be found at https://simtk.org/home/opensim_muscle.

Equilibrium Musculotendon 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).

Damped Equilibrium Musculotendon Model

The singularities in Eq. (2) arise because Eq. (1) is formulated in such a way that prevents the muscle from satisfying the equilibrium equation when it is deactivated (i.e., $//<![CDATA[ \begin{array}{l}a = 0\end{array} //]]>$) or when a nonzero tendon force is applied to a maximally-pennated muscle (i.e., $//<![CDATA[ \begin{array}{l}f^{\,\mathrm{T}} > 0\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}\alpha = 90^{\circ}\end{array} //]]>$). We address these two problems by limiting the maximum pennation angle and introducing a damper (with damping coefficient $//<![CDATA[ \begin{array}{l}\beta\end{array} //]]>$) in parallel with the contractile element, which results in the damped equilibrium musculotendon model:

$$//<![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) + \beta \tilde{v}^{\,\mathrm{M}} \right) \cos \alpha - f_o^{\,\mathrm{M}} \mathbf{f}^{\mathrm{T}} \big( \tilde{\ell}^{\,\mathrm{T}} \big) = 0\qquad(3)\end{array} //]]>$$

Since muscle length $//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{M}}\end{array} //]]>$ is a state, Eq. (3) can be readily and uniquely solved for $//<![CDATA[ \begin{array}{l}\tilde{v}^{\,\mathrm{M}}\end{array} //]]>$ using a derivative-based root-finding algorithm such as Newton's method, provided all the force development curves are C₂-continuous. The damped equilibrium musculotendon model of Eq. (3) should generate force profiles that are similar to those generated by the equilibrium model described by Eq. (1), but in a fraction of the simulation time because the damped equilibrium model is free of numerical singularities.

Rigid-Tendon Musculotendon Model

Some tendons are so stiff that they can be treated as inextensible, effectively replacing the tendon spring with an inextensible cable. The tendon inextensibility assumption is appropriate only when the tendon does not stretch sufficiently to affect the normalized length of the contractile element. This modeling simplification makes it possible to determine the muscle length ($//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{M}}\end{array} //]]>$) and velocity ($//<![CDATA[ \begin{array}{l}v^{\mathrm{M}}\end{array} //]]>$) from the musculotendon length ($//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{MT}}\end{array} //]]>$) and velocity ($//<![CDATA[ \begin{array}{l}v^{\,\mathrm{MT}}\end{array} //]]>$) using a kinematic model of the musculotendon actuator:

$$//<![CDATA[ \begin{array}{l}\displaystyle \ell^{\,\mathrm{MT}} = \ell^{\mathrm{T}} + \ell^{\,\mathrm{M}} \cos \alpha\qquad(4)\end{array} //]]>$$

Differentiating Eq. (4) with respect to time yields a relation between the muscle, tendon, and musculotendon actuator velocities:

$$//<![CDATA[ \begin{array}{l}\displaystyle v^{\,\mathrm{MT}} = v^{\mathrm{T}} + v^{\,\mathrm{M}} \cos \alpha - \ell^{\,\mathrm{M}} \dot{\alpha} \sin \alpha\qquad(5)\end{array} //]]>$$

where $//<![CDATA[ \begin{array}{l}v^{\mathrm{T}} = 0\end{array} //]]>$ if the tendon is rigid. The length of the muscle ($//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{M}}\end{array} //]]>$) and its orientation ($//<![CDATA[ \begin{array}{l}\alpha\end{array} //]]>$) are coupled by a fixed-height-parallelogram pennation model:

$$//<![CDATA[ \begin{array}{l}\displaystyle \ell^{\,\mathrm{M}} \sin \alpha = h\qquad(6)\end{array} //]]>$$

The constant height of the parallelogram ($//<![CDATA[ \begin{array}{l}h\end{array} //]]>$) is computed using the optimal muscle length and pennation angle:

$$//<![CDATA[ \begin{array}{l}\displaystyle h = \ell_o^{\,\mathrm{M}} \sin \alpha_o\end{array} //]]>$$

Differentiating Eq. (6) with respect to time yields an expression that can be used to calculate the pennation angular velocity:

$$//<![CDATA[ \begin{array}{l}\displaystyle \dot{\alpha} = - \frac{ v^{\,\mathrm{M}} \sin \alpha }{ \ell^{\,\mathrm{M}} \cos \alpha }\qquad(7)\end{array} //]]>$$

Since the tendon length and velocity are known (i.e., $//<![CDATA[ \begin{array}{l}\ell^{\mathrm{T}} = \ell_{\mathrm{s}}^{\mathrm{T}}\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}v^{\mathrm{T}} = 0\end{array} //]]>$), we use Eqs. (4) and (6) to solve for muscle length ($//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{M}}\end{array} //]]>$) given musculotendon length ($//<![CDATA[ \begin{array}{l}\ell^{\,\mathrm{MT}}\end{array} //]]>$), and Eqs. (5) and (7) to solve for muscle velocity ($//<![CDATA[ \begin{array}{l}v^{\,\mathrm{M}}\end{array} //]]>$) given musculotendon velocity ($//<![CDATA[ \begin{array}{l}v^{\,\mathrm{MT}}\end{array} //]]>$). As with the elastic-tendon models, we use a unilateral constraint to enforce a lower bound of $//<![CDATA[ \begin{array}{l}\ell_{\mathrm{min}}^{\,\mathrm{M}}\end{array} //]]>$ on muscle length. We compute the force generated by the muscle directly:

$$//<![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) + \beta \tilde{v}^{\,\mathrm{M}} \right) \cos \alpha\end{array} //]]>$$

Since a muscle can generate only tensile force, we constrain this equation to remain positive:

$$//<![CDATA[ \begin{align} f^{\,\mathrm{M}} = \begin{cases}f^{\,\mathrm{M*}} &: f^{\,\mathrm{M*}} > 0\\0 &: \text{otherwise}\end{cases} \end{align} //]]>$$

References

1. 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.
2. 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.
3. Winters, J.M. (1995) An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models. Annals of Biomedical Engineering, 23(4):359–374.
4. Zajac, F.E. (1989) Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering, 17(4):359–411.