Introduction

This page is a description of a Hill-type muscle model implemented in OpenSim that is based on the work of Thelen (2003). The model was modified by Matt Millard, Ajay Seth, and Peter Loan, and it is implemented in OpenSim as OpenSim::Thelen2003Muscle. In previous versions of OpenSim (2.4 and earlier), there were numerical errors in the implementation.

Background

The force-producing properties of muscle are complex and nonlinear (see McMahon (1984) for review) (Fig. 1). For simplicity, lumped-parameter, dimensionless muscle models, capable of representing a range of muscles with different architectures, are most commonly used in the dynamic simulation of movement (Zajac, 1989). In a complex musculoskeletal model, the model can be actuated by 50 or more muscle–tendon units, each of which is represented as a Hill-type contractile element in series with a tendon.

The Thelen2003Muscle model uses a standard equilibrium muscle model based on the Hill model. The muscle–tendon complex consists of three components: a contractile element (CE), a parallel element (PE), and a series element (SE). The muscle force generated is a function of three factors: the activation value (a), the normalized length of the muscle unit, and the normalized velocity of the muscle unit. The functions describing the force generated by a muscle as its length varies are called the active length curve (AL) for the contractile element and the passive length curve (PL) for the parallel element.

The parameters used to characterize each muscle are maximum isometric force, optimal muscle fiber length, tendon slack length, maximum contraction velocity, and pennation angle. For a table of muscle and tendon parameters, consult Anderson and Pandy (1999). During a forward dynamic simulation, the muscle force is calculated using two states: the activation value and the muscle fiber length.

Figure 1 from Thelen (2003).

Implementation

Using Newton's third law (under the assumption that the muscle and tendon units are massless), the differential equation of the muscle–tendon element is

\[ f_{iso} \left(a\left(t\right)f_{AL}\left(l^M\,\right)f_{v}\big(\dot{l}^M\,\big) + f_{PL}\left(l^M\,\right)\right) cos\alpha - f_{iso}\,f_{SE}\left(l^T\,\right) = 0\]

 

Rearranging the above equation and solving for \( f_{v}\big(\dot{l}^{M}\,\big)\) yields

 

\[ f_{v}\big(\dot{l}^{M}\,\big) = \frac{\frac{f_{SE}\left(l^{T}\right)}{cos\alpha} - f_{PL}\left(l^{M}\right)}{a\left(t\right)f_{AL}\left(l^{M}\,\right)}\]

 

The force–velocity curve is usually inverted to compute the fiber velocity:

 

\[ \dot{l}^{M} = f_{v}^{-1} \left\{ \frac{\frac{f_{SE}\left(l^T\right)}{cos\alpha} - f_{PL}\left(l^{M}\right)}{a\left(t\right)f_{AL}\left(l^{M}\,\right)} \right\}\]

 

which is then integrated to simulate the musculotendon dynamics. In general, the previous equation has 4 singularity conditions:
  1. \( a\left(t\right) \rightarrow 0\)

     

  2. \( f_{AL}\left(l^{M}\,\right) \rightarrow 0\)

     

  3. \( \alpha \rightarrow 90°\)

     

  4. \( f_{v}\big(\dot{l}^{M}\,\big) \leq 0\)  or  \( f_{v}\big(\dot{l}^{M}\,\big) \geq F^{m}_{len}\)

     

This implementation has been modified from the model presented in Thelen (2003) to avoid some of the numerical singularities:

  1. \( a\left(t\right) \rightarrow a_{min} \gt 0\)     A modified activation dynamic equation is used (MuscleFirstOrderActivationDynamicModel) that smoothly approaches a minimum value that is greater than zero.

  2. \( f_{AL}\left(l^{M}\,\right) \rightarrow 0\)    The active-force–length curve of the Thelen muscle is a Gaussian, which is always greater than 0.

  3. \( \alpha \rightarrow 90°\)   This singularity cannot be removed without changing the first equation, and still exists in the present Thelen2003Muscle implementation.

  4. \( f_{v}\big(\dot{l}^{M}\,\big) \leq 0\) or \( f_{v}\big(\dot{l}^{M}\,\big) \geq F^{m}_{len}\)  Equation 6 in Thelen (2003) has been modified so that \( V^{M}\) is linearly extrapolated when \( F^{M} \lt 0\) (during a concentric contraction) and when \( F^{M} \gt 0.95 F^{M}_{len}\) (during an eccentric contraction). These two modifications make the force–velocity curve invertible. The original force–velocity curve published by Thelen was not invertible.

  5. A unilateral constraint has been implemented to prevent the fiber from approaching a fiber length that is smaller than 0.01*optimal fiber length, or a fiber length that creates a pennation angle greater than the maximum pennation angle specified by the pennation model. Note that this unilateral constraint does not prevent the muscle fiber from becoming shorter than is physiologically possible (i.e., shorter than approximately half the normalized fiber length).

References

  1. McMahon, T.A. (1984) Muscles, Reflexes, and Locomotion. Princeton University Press, Princeton, New Jersey.

  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. Download PDF here.
  3. 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.