4 Mechanical Models II

4.5 Point-to-point muscles, tendons and ligaments

Point-to-point muscles are a simple type of component in biomechanical models that provide muscle-activated forces acting along a line between two points. ArtiSynth provides this through Muscle, which is a subclass of AxialSpring that generates an active muscle force in response to its excitation property. The excitation property can be set and queried using the methods

   setExcitation (double excitation)
   double getExcitation()

As with AxialSprings, Muscle components use subclasses of AxialMaterial to compute the applied force f(l,\dot{l},a) in response to the muscle’s length l, length velocity \dot{l}, and excitation signal a, which is assumed to lie in the interval a\in[0,1]. Special muscle subclasses of AxialMaterial exist that compute forces that vary in response to the excitation. As with other axial materials, this is done by the material’s computeF() method, first described in Section 3.1.4:

  double computeF (l, ldot, l0, excitation)

Usually the force is the sum of a passive component plus an active component that arises in response to the excitation signal.

Once a muscle material is created, it can be assigned to a muscle or queried using Muscle methods

  setMaterial (AxialMaterial mat)
  AxialMaterial getMaterial()

Muscle materials can also be assigned to axial springs, although the resulting force will always be computed with 0 excitation.

4.5.1 Simple muscle materials

A number of simple muscle materials are described below. All compute force as a function of a and l, with an optional damping force that is proportional to \dot{l}. More complex Hill-type muscles, with force velocity curves, pennation angle, and a tendon component in series, are also available and described in Section 4.5.3.

For historical reasons, the materials ConstantAxialMuscle, LinearAxialMuscle and PeckAxialMuscle, described below, contain a property called forceScaling that uniformly scales the computed force. This property is now deprecated, and should therefore have a value of 1. However, creating these materials with constructors not indicated in the documentation below will cause forceScaling to be set to a default value of 1000, thus requiring that the maximum force and damping values be correspondingly reduced.

4.5.1.1 SimpleAxialMuscle

SimpleAxialMuscle is the default AxialMaterial for Muscle, and is essentially an activated version of LinearAxialMaterial. It computes a simple force according to

f(l,\dot{l})=k(l-l_{0})+d\dot{l}+f_{max}\,a, (4.1)

where l_{0} is the muscle rest length, k and d are stiffness and damping terms, and f_{max} is the maximum excitation force. l_{0} is specified by the restLength property of the muscle component, a by the excitation property of the Muscle, and k, d and f_{max} by the following properties of SimpleAxialMuscle:

Property Description Default
k stiffness stiffness term 0
d damping damping term 0
f_{max} maxForce maximum activation force that can be imparted by a 1

SimpleAxialMuscles can be created with the following constructors:

SimpleAxialMuscle()
SimpleAxialMuscle (stiffness, damping, fmax)

where stiffness, damping, and fmax specify k, d and f_{max}, and properties are otherwise set to their defaults.

4.5.1.2 ConstantAxialMuscle

ConstantAxialMuscle is a simple muscle material that has a contractile force proportional to its activation, a constant passive tension, and linear damping. The resulting force is described by:

f(\dot{l})=f_{max}(a+f_{p})+d\dot{l}. (4.2)

The parameters f_{max}, f_{p} and d are specified as properties of ConstantAxialMuscle:

Property Description Default
f_{max} maxForce maximum contractile force 1
f_{p} passiveFraction proportion of f_{max} to apply as passive tension 0
d damping damping parameter 0

ConstantAxialMuscles can be created with the following factory methods and constructors:

ConstantAxialMuscle.create()
ConstantAxialMuscle.create (fmax)
ConstantAxialMuscle (fmax, pfrac)
ConstantAxialMuscle (fmax, pfrac, damping)

where fmax, pfrac, and damping specify f_{max}, f_{p} and d, and properties are otherwise set to their defaults.

Creating a ConstantAxialMuscle with the no-args constructor, or another constructor not listed above, will cause its forceScaling property (Section 4.5.1) to be set to 1000 instead of 1, thus requiring that fmax and damping be corresponding reduced.

4.5.1.3 LinearAxialMuscle

LinearAxialMuscle is a simple muscle material that has a linear relationship between length and tension, as well as linear damping. Given a normalized length described by

\hat{l}=\frac{l-l_{opt}}{l_{max}-l_{opt}},\quad\mathrm{with}~{}0\leq\hat{l}%
\leq 1~{}\mathrm{enforced}, (4.3)

the force generated by this material is:

f(l,\dot{l})=f_{max}\left(a\hat{l}+f_{p}\hat{l}\,\right)+d\dot{l}. (4.4)

The parameters are specified as properties of LinearAxialMaterial:

Property Description Default
f_{max} maxForce maximum contractile force 1
f_{p} passiveFraction proportion of f_{max} that forms the maximum passive force 0
l_{max} maxLength length beyond which maximum passive and active forces are generated 1
l_{opt} optLength length below which zero active force is generated 0
d damping damping parameter 0

LinearAxialMuscles can be created with the following factory methods and constructors:

LinearAxialMuscle.create()
LinearAxialMuscle.create (fmax, lrest)
LinearAxialMuscle (fmax, lopt, lmax, pfrac)
LinearAxialMuscle (fmax, lopt, lmax, pfrac, damping)

where fmax, lopt, lmax, pfrac and damping specify f_{max}, l_{opt}, l_{max}, f_{p} and d, lrest specifies l_{opt} and l_{max} via l_{opt}=\mathrm{\tt lrest} and l_{max}=3/2\;\mathrm{\tt lrest}, and other properties are set to their defaults.

Creating a LinearAxialMuscle with the no-args constructor, or another constructor not listed above, will cause its forceScaling property (Section 4.5.1) to be set to 1000 instead of 1, thus requiring that fmax and damping be corresponding reduced.

4.5.1.4 PeckAxialMuscle

The PeckAxialMuscle material generates a force as described in [19]. It has a typical Hill-type active force-length relationship (modeled as a cosine), but the passive force-length properties are linear. This muscle model was empirically verified for jaw muscles during wide jaw opening.

Figure 4.5: Left: active (red) and passive force length curve (green) for a PeckAxialMuscle with f_{max}=1, f_{p}=1, T_{r}=0, l_{opt}=1, l_{max}=2, and a=1. Note that the passive force curve is linear for l\in[l_{opt},l_{max}] and saturates for l>l_{max}. Right: combined active and passive force (blue).

Given a normalized fibre length described by

\hat{l}_{f}=\frac{l-l_{opt}T_{r}}{l_{opt}(1-T_{r})},\quad\mathrm{with}~{}1/2%
\leq\hat{l}\leq 3/2~{}\mathrm{enforced}, (4.5)

and a normalized muscle length

\hat{l}_{m}=\frac{l-l_{opt}}{l_{max}-l_{opt}},\quad\mathrm{with}~{}0\leq\hat{l%
}\leq 1~{}\mathrm{enforced}, (4.6)

the force generated by this material is:

f(l,\dot{l})=f_{max}\left(a\frac{1+\cos(2\pi\hat{l}_{f})}{2}+f_{p}\hat{l}_{m}%
\right)+d\dot{l}. (4.7)

The parameters are specified as properties of PeckAxialMuscle:

Property Description Default
f_{max} maxForce maximum contractile force 1
f_{p} passiveFraction proportion of f_{max} that forms the maximum passive force 0
T_{r} tendonRatio tendon to fibre length ratio 0
l_{max} maxLength length at which maximum passive force is generated 1
l_{opt} optLength length at which maximum active force is generated 0
d damping damping parameter 0

Figure 4.5 illustrates the force length relationship for a PeckAxialMuscle.

PeckAxialMuscles can be created with the following factory methods:

PeckAxialMuscle.create()
PeckAxialMuscle.create (fmax, lopt, lmax, tratio, pfrac)
PeckAxialMuscle.create (fmax, lopt, lmax, tratio, pfrac, damping)

where fmax, lopt, lmax, tratio, pfrac and damping specify f_{max}, l_{opt}, l_{max}, T_{r}, f_{p} and d, and other properties are set to their defaults.

Creating a PeckAxialMuscle with the no-args constructor, or another constructor not listed above, will cause its forceScaling property (Section 4.5.1) to be set to 1000 instead of 1, thus requiring that fmax and damping be corresponding reduced.

4.5.1.5 BlemkerAxialMuscle

The BlemkerAxialMuscle material generates a force as described in [5]. It is the axial muscle equivalent to the constitutive equation along the muscle fiber direction specified in the BlemkerMuscle FEM material.

Figure 4.6: Left: active force curve f_{L}(\hat{l}) (red) and passive force curve f_{P}(\hat{l}) (green) for a BlemkerAxialMuscle with l_{opt}=1, l_{max}=1.4, P_{1}=0.05, and P_{2}=6.6. Right: total force length curve (blue) with f_{max}=1 and a=1.

The force produced is a combination of active and passive terms, plus a damping term, given by

f(l,\dot{l})=f_{max}\left(af_{L}(\hat{l})+f_{P}(\hat{l})\right)+d\dot{l}, (4.8)

where f_{max} is the maximum force, \hat{l}\equiv l/l_{opt} is the normalized muscle length, and f_{L}(\hat{l}) and f_{P}(\hat{l}) are the active and passive force length curves, given by

f_{L}(\hat{l})=\begin{cases}9\,(\hat{l}-0.4),&\hat{l}\in[0.4,0.6]\\
1-4(1-\hat{l}),&\hat{l}\in[0.6,1.4]\\
9\,(\hat{l}-1.6),&\hat{l}\in[1.4,1.6]\\
0,&\mathrm{otherwise},\\
\end{cases} (4.9)

and

f_{P}(\hat{l})=\begin{cases}0,&\hat{l}<1\\
P_{1}(e^{P_{2}(\hat{l}-1)}-1),&\hat{l}\in[l_{opt},l_{max}/l_{opt}]\\
P_{3}\hat{l}+P_{4},&\hat{l}>l_{max}/l_{opt}.\\
\end{cases} (4.10)

For the passive force length curve, P_{3} and P_{4} are computed to provide linear extrapolation for \hat{l}>l_{max}/l_{opt}. The other parameters are specified by properties of BlemkerAxialMuscle:

Property Description Default
f_{max} maxForce the maximum contractile force 3\times 10^{5}
P_{1} expStressCoeff exponential stress coefficient 0.05
P_{2} uncrimpingFactor fibre uncrimping factor 6.6
l_{max} maxLength length at which passive force becomes linear 1.4
l_{opt} optLength length at which maximum active force is generated 1
d damping damping parameter 0

Figure 4.6 illustrates the force length relationship for a BlemkerAxialMuscle.

BlemkerAxialMuscles can be created with the following constructors,

BlemkerAxialMuscle()
BlemkerAxialMuscle (lmax, lopt, fmax, ecoef, uncrimp)

where lmax, lopt, fmax, ecoef, and uncrimp specify l_{max}, l_{opt}, f_{max}, P_{1}, P_{2} and d, and properties are otherwise set to their defaults.

4.5.2 Example: muscle attached to a rigid body

Figure 4.7: SimpleMuscle model loaded into ArtiSynth.

A simple model showing a single muscle connected to a rigid body is defined in

  artisynth.demos.tutorial.SimpleMuscle

This model is identical to RigidBodySpring described in Section 3.2.2, except that the code to create the spring is replaced with code to create a muscle with a SimpleAxialMuscle material:

      // create the muscle:
      muscle = new Muscle ("mus", /*restLength=*/0);
      muscle.setPoints (p1, mkr);
      muscle.setMaterial (
         new SimpleAxialMuscle (/*stiffness=*/20, /*damping=*/10, /*fmax=*/10));

Also, so that the muscle renders differently, the rendering style for lines is set to SPINDLE using the convenience method

      RenderProps.setSpindleLines (muscle, 0.02, Color.RED);

To run this example in ArtiSynth, select All demos > tutorial > SimpleMuscle from the Models menu. The model should load and initially appear as in Figure 4.7. Running the model (Section 1.5.3) will cause the box to fall and sway under gravity. To see the effect of the excitation property, select the muscle in the viewer and then choose Edit properties ... from the right-click context menu. This will open an editing panel that allows the muscle’s properties to be adjusted interactively. Adjusting the excitation property using the adjacent slider will cause the muscle force to vary.

4.5.3 Equilibrium muscles

ArtiSynth supplies several muscle materials that simulate a pennated muscle in series with a tendon. Both the muscle and tendon generate forces which are functions of their respective lengths l_{m} and l_{t}, and because these components are in series, their respective forces must be equal when the system is in equilibrium. Given an overall muscle-tendon length l\equiv l_{m}+l_{t}, ArtiSynth solves for l_{m} at each time step to ensure that this equilibrium condition is met.

Figure 4.8: Schematic illustration of an pennated muscle in series with a tendon.

A general muscle-tendon system is illustrated by Figure 4.8, where l_{m} and l_{t} are the muscle and tendon lengths. These two components generate forces given by f_{m}(l_{m},v_{m}) and f_{t}(l_{t}), where v_{m}\equiv\dot{l}_{m}, and the fact that these components are in series implies that their forces must be in equilibrium:

f_{m}(l_{m},v_{m})=f_{t}(l_{t}). (4.11)

The muscle force is in turn produced by a fibre force f_{f} acting at an pennation angle \phi with respect to the principal muscle direction, such that

f_{m}=\cos(\phi)f_{f}(l_{f},v_{f}),

where l_{f} is the fibre length, which satisfies l_{m}=\cos(\phi)l_{f}, and v_{f}\equiv\dot{l}_{f}.

The fibre force is usually computed from the normalized fibre length \hat{l}_{f} and normalized fibre velocity \hat{v}_{f}, defined by

\hat{l}_{f}=\frac{l_{f}}{l_{o}},\quad\hat{v}_{f}=\frac{v_{f}}{l_{o}V_{m}}, (4.12)

where l_{o} is the optimal fibre length and V_{m} is the maximum contraction velocity. It is composed of three components in parallel, such that

f_{f}(\hat{l}_{f},\hat{v}_{f})=F_{o}\left(af_{L}(\hat{l}_{f})f_{V}(\hat{v}_{f}%
)+f_{P}(\hat{l}_{f})+\beta\hat{v}_{f}\right), (4.13)

where F_{o} is the maximum isometric force, af_{L}(\hat{l}_{f})f_{V}(\hat{v}_{f}) is the active force term induced by an activation level a and modulated by the active force length curve f_{L}(\hat{l}_{f}) and the force velocity curve f_{V}(\hat{v}_{f}), f_{P}(\hat{l}_{f}) is the passive force length curve, and \beta\hat{v}_{f} is an optional damping term induced by a fibre damping parameter \beta.

The tendon force f_{t}(\hat{l}_{t}) is computed from

f_{t}(l_{t})=F_{o}f_{T}(\hat{l}_{t}),

where F_{o} is (again) the maximum isometric force, f_{T}(\hat{l}_{t}) is the tendon force length curve and \hat{l}_{t} is the normalized tendon length, defined by dividing l_{t} by the tendon slack length T:

\hat{l}_{t}=\frac{l_{t}}{T}.

As the muscle moves, it is assumed that the height H from the fibre origin to the main muscle line of action remains constant. This height is defined by

H=l_{o}\sin\phi_{o},

where \phi_{o} is the optimal pennation angle at the optimal fibre length when l_{f}=l_{o}.

4.5.4 Equilibrium muscle materials

The equilibrium muscle materials supplied by ArtiSynth include Thelen2003AxialMuscle and Millard2012AxialMuscle. These are all controlled by properties which specify the parameters presented in Section 4.5.3:

Property Description Default
F_{o} maxIsoForce maximum isometric force 1000
l_{o} optFibreLength optimal fibre length 0.1
T tendonSlackLength length beyond which tendon exerts force 0.2
\phi_{o} optPennationAngle pennation angle at optimal fibre length (radians) 0
V_{m} maxContractionVelocity maximum contraction velocity 10
\beta fibreDamping damping parameter for normalized fibre velocity 0

The materials differ from each other with respect to their active, passive and tendon force length curves (f_{L}(\hat{l}_{f}), f_{P}(\hat{l}_{f}), and f_{T}(\hat{l}_{t})) and their force velocity curves (f_{V}(\hat{v}_{f})).

For a given muscle instance, it is typically only necessary to specify F_{o}, l_{o}, T and \phi_{o}, where F_{o}, l_{o} and T are given in the model’s basic force and length units. V_{m} is given in units of l_{o} per second and has a default value of 10; changing this value will stretch or contract the domain of the force velocity curve f_{V}(\hat{v}_{f}). The damping parameter \beta imparts a damping force proportional to \hat{v}_{f} and has a default value of 0 (i.e., no damping). It is often not necessary to add fibre damping (since damping can be readily applied to the model in other ways), but \beta is supplied for formulations that do specify damping. If non-zero, \beta is usually set to a value less than or equal to 0.1.

In addition to the above parameters, equilibrium muscle materials also export the following properties to adjust their behavior:

Property Description Default
rigidTendon forces the tendon to be rigid false
ignoreForceVelocity ignore the force velocity curve f_{V}(\hat{v}_{f}) false

If rigidTendon is true, the tendon will be assumed to be rigid with a length given by the tendon slack length parameter T. This simplifies the muscle computations since l_{m} is then given by l_{m}=l-T and there is no need to compute an equilibrium position. If ignoreForceVelocity is true, then force velocity effects are ignored by replacing the force velocity curve f_{V}(\hat{v}_{f}) with 1.

The different equilibrium muscle materials are now summarized:

4.5.4.1 Millard2012AxialMuscle

Millard2012AxialMuscle implements the default version of the Millard2012EquilibriumMuscle model supplied by OpenSim [7] and described in [15].

The active, passive and tendon force length curves (f_{L}(\hat{l}_{f}), f_{P}(\hat{l}_{f}), and f_{T}(\hat{l}_{t})) and force velocity curve (f_{V}(\hat{v}_{f})) are implemented using cubic Hermite spline curves to conform closely to the default curve values provided by OpenSim’s Millard2012EquilibriumMuscle. Plots of these curves are shown in Figures 4.9 and 4.10. Both the passive and tendon force curves are linearly extrapolated (and hence exhibit constant stiffness) past \hat{l}_{f}=1.7 and \hat{l}_{t}=1.049, respectively, with stiffness values of 2.857 and 29.06.

Figure 4.9: Default active force length curve f_{L}(\hat{l}_{f}) (left) and passive force length curve f_{P}(\hat{l}_{f}) (right) for the Millard 2012 muscle material.
Figure 4.10: Default force velocity curve f_{V}(\hat{v}_{f}) (left) and tendon force length curve f_{T}(\hat{l}_{t}) (right) for the Millard 2012 muscle material. Note that the tendon force curve is about 10 times stiffer than the passive force curve, as seen by that fact the tendon curve is shown with a horizontal range of [0.9,1.1] vs. [0,2] for the passive curve.

OpenSim requires that the Millard2012EquilibriumMuscle always exhibits a small bias activation, even when the activation should be zero, in order to avoid a singularity in the computation of the muscle length. ArtiSynth computes the muscle length in a manner that makes this unnecessary.

Millard2012AxialMuscles can be created using the constructors

Millard2012AxialMuscle()
Millard2012AxialMuscle(fmax, lopt, tslack, optPenAng)

where fmax, lopt, tslack, and optPenAng specify F_{o}, l_{o}, T, and \phi_{o}, and other properties are set to their defaults.

4.5.4.2 Thelen2003AxialMuscle

Thelen2003AxialMuscle implements the Thelen2003Muscle model supplied by OpenSim [7] and introduced in [26].

The active and passive force length curves are described by

f_{L}(\hat{l}_{f})=e^{-(\hat{l}_{f}-1)^{2}/\gamma} (4.14)

and

f_{P}(\hat{l}_{f})=\frac{e^{k^{PE}(\hat{l}_{f}-1)/\epsilon_{0}^{M}}-1}{e^{k^{%
PE}-1}},` (4.15)

where \gamma, k^{PE}, and \epsilon_{0}^{M} are parameters, described in [26], that control the curve shapes. These are exposed as properties of Thelen2003AxialMuscle and are described in Table 4.3.

The tendon force length curve is described by

f_{T}(\hat{l}_{t})=\begin{cases}0,&\hat{l}_{t}<1\\
\dfrac{F_{toe}}{e^{k_{toe}}-1}(e^{k_{toe}(\hat{l}_{t}-1)/\epsilon^{T}_{toe}}-1%
),&\hat{l}_{t}\leq 1+\epsilon^{T}_{toe}\\
k_{lin}(\hat{l}_{t}-1-\epsilon^{T}_{toe})+F_{toe},&\hat{l}_{t}>1+\epsilon^{T}_%
{toe},\end{cases} (4.16)

where F_{toe}, k_{toe}, k_{lin} and \epsilon^{T}_{toe} are parameters, described in [26], that control the curve shape. The values of these are either fixed, or derived from the maximum isometric tendon strain \epsilon_{0}^{T} (controlled by the fmaxTendonStrain property, Table 4.3), according to

F_{toe}=0.33,\quad k_{toe}=3.0,\quad\epsilon^{T}_{toe}=\frac{99\epsilon_{0}^{T%
}e^{3}}{166e^{3}-67},\quad k_{lin}=\frac{0.67}{\epsilon_{0}^{T}-\epsilon^{T}_{%
toe}}. (4.17)
Property Description Default
\gamma kShapeActive shape factor for active force length curve 0.45
k^{PE} kShapePassive shape factor for active force length curve 5.0
\epsilon_{0}^{M} fmaxMuscleStrain passive muscle strain at maximum isometric muscle force 0.6
\epsilon_{0}^{T} fmaxTendonStrain tendon strain at maximum isometric muscle force 0.04
A_{f} af force velocity shape factor 0.25
\bar{F}^{M}_{len} flen maximum normalized lengthening force 1.4
fvLinearExtrapThreshold \hat{v}_{f} beyond which f_{V}(\hat{v}_{f}) is linearly extrapolated 0.95
Table 4.3: Properties of Thelen2003AxialMuscle that control the shapes of the force curves.

The force velocity curve is determined from equation (6) in [26]. In the notation used by that equation and the accompanying equation (7), V^{M}/V^{M}_{\text{max}}=\hat{v}_{f}, f_{l}=f_{L}(\hat{l}_{f}), and \bar{F}^{M}=af_{L}(\hat{l}_{f})f_{V}(\hat{v}_{f}). Inverting equation (6) yields the force velocity curve:

f_{v}=\begin{cases}\dfrac{\alpha+\bar{v}^{M}}{\alpha-\bar{v}^{M}/A_{f}},&\bar{%
v}^{M}\leq 0\\
\dfrac{\beta+\bar{v}^{M}\bar{F}^{M}_{len}}{\beta+\bar{v}^{M}},&\bar{v}^{M}>0,%
\end{cases}

where

\alpha\equiv 0.25+0.75a,\qquad\beta\equiv\alpha\frac{\bar{F}^{M}_{len}-1}{2+2/%
A_{f}},

and a is the activation level. Parameters include the maximum normalized lengthening force \bar{F}^{M}_{len}, and force velocity shape factor A_{f}, described in Table 4.3.

Plots of the curves resulting from the above equations are shown, for default values, in Figures 4.11 and 4.12.

Figure 4.11: Default active force length curve f_{L}(\hat{l}_{f}) (left) and passive force length curve f_{P}(\hat{l}_{f}) (right) for the Thelen 2003 muscle material. Note that the passive curve is exponential and does not transition to a constant slope for high values of \hat{l}_{f}.
Figure 4.12: Default force velocity curve f_{V}(\hat{v}_{f}) (left) and tendon force length curve f_{T}(\hat{l}_{t}) (right) for the Thelen 2003 muscle model. Since the force velocity curve also depends on the activation, its plot shows the curve for two activation values: 1.0 (dark magenta), and 0.5 (light magenta).

Thelen2003AxialMuscles can be created using the constructors

Thelen2003AxialMuscle()
Thelen2003AxialMuscle(fmax, lopt, tslack, optPenAng)

where fmax, lopt, tslack, and optPenAng specify F_{o}, l_{o}, T, and \phi_{o}, and other properties are set to their defaults.

4.5.5 Tendons and ligaments

Special point-to-point spring materials are also available to model tendons and ligaments. Because these are passive materials, they would normally be assigned to AxialSpring instead of Muscle components, although they will also work correctly in the latter. The materials include:

4.5.5.1 Millard2012AxialTendon

Millard2012AxialTendon implements the default tendon material of the Millard2012AxialMuscle (Section 4.5.4.1), with a force length relationship given by

f(l)=F_{o}f_{T}(\hat{l}_{t}),\quad\hat{l}_{t}\equiv\frac{l}{T}, (4.18)

where F_{o} is the maximum isometric force, f_{T}() is the tendon force length curve (shown in Figure 4.10, right), and T is the tendon slack length. F_{o} and T are specified by the following properties of Millard2012AxialTendon:

Property Description Default
F_{o} maxIsoForce maximum isometric force 1000
T tendonSlackLength length beyond which tendon exerts force 0.2

Millard2012AxialTendons can be created with the constructors

Millard2012AxialTendon()
Millard2012AxialTendon (fmax, tslack)

where fmax and tslack specify F_{o} and T, and other properties are set to their defaults.

4.5.5.2 Thelen2003AxialTendon

Thelen2003AxialTendon implements the default tendon material of the Thelen2003AxialMuscle (Section 4.5.4.2), with a force length relationship given by

f(l)=F_{o}f_{T}(\hat{l}_{t}),\quad\hat{l}_{t}\equiv\frac{l}{T}, (4.19)

where F_{o} is the maximum isometric force, f_{T}() is the tendon force length curve described by equation (4.16), and T is the tendon slack length. F_{o}, T, and the maximum isometric tendon strain \epsilon_{0}^{T} (used to determine the parameters of (4.16), according to (4.17)), are specified by the following properties of Thelen2003AxialTendon:

Property Description Default
F_{o} maxIsoForce maximum isometric force 1000
T tendonSlackLength length beyond which tendon exerts force 0.2
\epsilon_{0}^{T} fmaxTendonStrain tendon strain at maximum isometric muscle force 0.04

Thelen2003AxialTendons can be created with the constructors

Thelen2003AxialTendon()
Thelen2003AxialTendon (fmax, tslack)

where fmax and tslack specify F_{o} and T, and other properties are set to their defaults.

4.5.5.3 Blankevoort1991AxialLigament

Blankevoort1991AxialLigament implements the Blankevoort1991Ligament model supplied by OpenSim [7] and described in [4, 23].

With the ligament strain and its derivative defined by

\epsilon\equiv\frac{l-l_{0}}{l_{0}},\quad\dot{\epsilon}\equiv\frac{\dot{l}}{l_%
{0}}, (4.20)

the ligament force is given by

f(l,\dot{l})=f_{e}(\epsilon)+f_{d}(\dot{\epsilon}), (4.21)

where

f_{e}(\epsilon)=\begin{cases}0,&\epsilon<0\\
\dfrac{1}{2\epsilon_{t}}k\epsilon^{2},&0\leq\epsilon\leq\epsilon_{t}\\
k(\epsilon-\dfrac{\epsilon_{t}}{2}),&\epsilon>\epsilon_{t},\end{cases}

and

f_{d}(\dot{\epsilon})=\begin{cases}d\dot{\epsilon},&\epsilon>0\;\mathrm{and}\;%
\dot{\epsilon}>0\\
0,&\mathrm{otherwise}.\end{cases}

The parameters l_{0}, \epsilon_{t}, k, and d are specified by the following properties of Blankevoort1991Ligament:

Property Description Default
l_{0} slackLength ligament slack length 0
\epsilon_{t} transitionStrain strain at transition from toe to linear 0.06
k linearStiffness maximum stiffness of force-length curve 100
d damping damping parameter 0.003

Blankevoort1991Ligaments can be created with the constructors

Blankevoort1991Ligament()
Blankevoort1991Ligament (stiffness, slackLen, damping)

where stiffness, slackLen and damping specify k, l_{o} and d, and other properties are set to their defaults.

4.5.6 Example: muscles with separate tendons

An alternate way to model the equilibrium muscles of Section 4.5.3 is to use two separate point-to-point muscles, attached in series via a connecting particle with a small mass (and possibly damping), thus implementing the muscle and tendon components separately. This approach allows the muscle equilibrium length to be determined automatically by the physical response of the connecting particle, in a manner similar to that employed by OpenSim’s Millard2012AccelerationMuscle (described in [14]). For the muscle material, one can use any of the tendonless materials of Section 4.5.1, or the materials of Section 4.5.3 with the properties rigidTendon and tendonSlackLength set to true and 0, respectively, while the tendon material can be one described in Section 4.5.5 or some other suitable passive material. The implicit integrators used by ArtiSynth should permit relatively stable simulation of the arrangement.

While it is usually easier to employ the equilibrium muscle materials of Section 4.5.4, especially for models involving wrapping and/or via points (Chapter 9) where it may be difficult to handle the connecting particle correctly, in some situations the separate muscle/tendon method may offer a more versatile solution. It can also be used to validate the equilibrium muscles, as illustrated by the application model defined in

  artisynth.demos.tutorial.EquilibriumMuscleDemo

which provides a direct comparison of the methods. It creates two simple instances of a Millard 2012 muscle, one above the other in the x-z plane, with the top instance implemented using the separate muscle/tendon method and the bottom instance implemented using an equilibrium muscle material (Figure 4.13). The application model uses control and monitoring components introduced in Chapter 5 to drive the simulation and record its results: a controller (Section 5.3) is used to uniformly extend the length of both muscles; output probes (Section 5.4) are used to record the resulting tension forces; and a control panel (Section 5.1) is created to allow the user to interactively adjust some of the muscle parameters.

Figure 4.13: EquilibriumMuscleDemo loaded into ArtiSynth, with different particles and muscle components labeled. The separate muscle/tendon muscle is at the top, the equilibrium muscle at the bottom, and the control panel is overlaid on the viewer.

The code for the EquilibriumMuscleDemo class attributes and build() method is shown below:

1    Particle pr0, pr1; // right end point particles; used by controller
2
3    // default muscle parameter settings
4    private double myOptPennationAng = Math.toRadians(20.0);
5    private double myMaxIsoForce = 10.0;
6    private double myTendonSlackLen = 0.5;
7    private double myOptFibreLen = 0.5;
8
9    // initial total length of the muscles:
10    private double len0 = 0.25 + myTendonSlackLen;
11
12    public void build (String[] args) {
13       // create a mech model with zero gravity
14       MechModel mech = new MechModel ("mech");
15       addModel (mech);
16       mech.setGravity (0, 0, 0);
17
18       // build first muscle, consisting of a tendonless muscle, attached to a
19       // tendon via a connecting particle pc0 with a small mass.
20       Particle pl0 = new Particle ("pl0", 1.0, 0.0, 0, 0); // left end point
21       pl0.setDynamic (false); // point is fixed
22       mech.addParticle (pl0);
23
24       // create connecting particle. x coordinate will be set later.
25       Particle pc0 = new Particle ("pc0", /*mass=*/1e-5, 0, 0, 0);
26       mech.addParticle (pc0);
27
28       pr0 = new Particle ("pr0", 1.0, len0, 0, 0); // right end point
29       pr0.setDynamic (false); // point will be positioned by length controller
30       mech.addParticle (pr0);
31
32       // create muscle and attach it between pl0 and pc0
33       Muscle mus0 = new Muscle("mus0"); // muscle
34       Millard2012AxialMuscle mat0 = new Millard2012AxialMuscle (
35          myMaxIsoForce, myOptFibreLen, myTendonSlackLen, myOptPennationAng);
36       mat0.setRigidTendon (true); // set muscle to rigid tendon with zero length
37       mat0.setTendonSlackLength (0);
38       mus0.setMaterial (mat0);
39       mech.attachAxialSpring (pl0, pc0, mus0);
40
41       // create explicit tendon and attach it between pc0 and pr0
42       AxialSpring tendon = new AxialSpring(); // tendon
43       tendon.setMaterial (
44          new Millard2012AxialTendon (myMaxIsoForce, myTendonSlackLen));
45       mech.attachAxialSpring (pc0, pr0, tendon);
46
47       // build second muscle, using combined muscle/tendom material, and attach
48       // it between pl1 and pr1.
49       Particle pl1 = new Particle (1.0, 0, 0, -0.5); // left end point
50       pl1.setDynamic (false);
51       mech.addParticle (pl1);
52
53       pr1 = new Particle ("pr1", 1.0, len0, 0, -0.5); // right end point
54       pr1.setDynamic (false);
55       mech.addParticle (pr1);
56
57       Muscle mus1 = new Muscle("mus1");
58       Millard2012AxialMuscle mat1 = new Millard2012AxialMuscle (
59          myMaxIsoForce, myOptFibreLen, myTendonSlackLen, myOptPennationAng);
60       mus1.setMaterial (mat1);
61       mech.attachAxialSpring (pl1, pr1, mus1);
62
63       // initialize both muscle excitations to 1, and then adjust the muscle
64       // lengths to the corresponding (zero velocity) equilibrium position
65       mus0.setExcitation (1);
66       mus1.setExcitation (1);
67       // compute equilibrium muscle length with for 0 velocity
68       double lm = mat1.computeLmWithConstantVm (
69          len0, /*vel=*/0, /*excitation=*/1);
70       // set muscle length of mat1 and x coord of pc0 to muscle length:
71       mat1.setMuscleLength (lm);
72       pc0.setPosition (new Point3d (lm, 0, 0));
73
74       // set render properties:
75       // render markers as white spheres, and muscles as red spindles
76       RenderProps.setSphericalPoints (mech, 0.03, Color.WHITE);
77       RenderProps.setSpindleLines (mech, 0.02, Color.RED);
78       // render tendon in blue and the juntion point in cyan
79       RenderProps.setLineColor (tendon, Color.BLUE);
80       RenderProps.setPointColor (pc0, Color.CYAN);
81
82       // create a control panel to adjust material parameters and excitation
83       ControlPanel panel = new ControlPanel();
84       panel.addWidget ("material.optPennationAngle", mus0, mus1);
85       panel.addWidget ("material.fibreDamping", mus0, mus1);
86       panel.addWidget ("material.ignoreForceVelocity", mus0, mus1);
87       panel.addWidget ("excitation", mus0, mus1);
88       addControlPanel (panel);
89
90       // add a controller to extend/contract the muscle end points, and probes
91       // to record both muscle forces
92       addController (new LengthController());
93       addForceProbe ("muscle/tendon force", mus0, 2);
94       addForceProbe ("equilibrium force", mus1, 2);
95    }

Lines 4-10 declare attributes for muscle parameters and the initial length of the combined muscle-tendon. Within the build() method, a MechModel is created with zero gravity (lines 14-16).

Next, the separate muscle/tendon muscle is assembled, consisting of three particles (pl0, pc0, and pr0, lines 20-30), a muscle mus0 connected between pr0 and pc0 (lines 33-39), and an tendon connected between pc0 and pr0 (lines 42-45). Particles pl0 and pr0 are both set to be non-dynamic, since pl0 will be fixed and pr0 will be moved parametrically, while the connecting particle pc0 has a small mass and its position is updated by the simulation and so automatically maintains force equilibrium between the muscle and tendon. The material mat0 used for mus0 is an instance of Millard2012AxialMuscle, with the tendon removed by setting the tendonSlackLength and rigidTendon properties to 0 and true, respectively, while the material for tendon is an instance of Millard2012AxialTendon.

The equilibrium muscle is then assembled, consisting of two particles (pl1 and pr1, lines 49-55) and a muscle mus1 connected between them (lines 57-61). pl1 and pr1 are both set to be non-dynamic, since pl1 will be fixed and pr1 will be moved parametrically, and the material mat1 for mus1 is an instance of Millard2012AxialMuscle in its default equilibrium mode. Excitations for both mus0 and mus1 are initialized to 1, and the zero-velocity equilibrium muscle length is computed using mat1.computeLmWithConstantVm() (lines 65-69). This is used to update the muscle position, for mus0 by setting the x coordinate of pc0 and for mus1 by setting the internal muscle length variable of mat1 (lines 71-72).

Render properties for different components are set at lines 76-80, and then a control panel is created to allow interactive adjustment of the muscle material properties optPennationAngle, fibreDamping, and ignoreForceVelocity, as well the muscle excitations (lines 83-88). As explained in Sections 5.1.1 and 5.1.3, addWidget() methods are used to create widgets that control each of these properties for both mus0 and mus1.

At line 92, a controller is added to the root model to move the right side particles p0r and p1r during the simulation, thus changing the muscles’ lengths and hence their tension forces. Controllers are explained in Section 5.3, and the controller itself is defined by lines 101-127 of the code listing below. At the start of each simulation time step, the controller’s apply() method is called, with t0 and t1 denoting the step’s start and end times. It increases the x positions of p0r and p1r uniformly with a speed of mySpeed until time myRunTime/2, and then reverses direction until time myRunTime. Velocities are also updated since these are needed to determine \dot{l} for the muscles’ computeF() methods.

Each muscle’s tension force is recorded by an output probe connected to its forceNorm property. Probes are explained in Section 5.4; in this example, they are created using the convenience method addForceProbe(), defined by lines 130-140 (below) and called at lines 93-94 in the build() method, with probes for the muscle/tendon and equilibrium muscles named “muscle/tendon force” and “equilibrium force”, respectively.

1    /**
2     * A controller to extend and the contract the muscle length by moving the
3     * rightmost muscle end points.
4     */
5    public class LengthController extends ControllerBase {
6
7       double myRunTime = 1.5; // total extensions/contraction time
8       double mySpeed = 1.0;   // speed of the end point motion
9
10       public LengthController() {
11          // need null args constructor if this model is read from a file
12       }
13
14       public void apply (double t0, double t1) {
15          double xlen = len0; // x position of the end points
16          double xvel = 0; // x velocity of the end points
17          if (t1 <= myRunTime/2) { // extend
18             xlen += mySpeed*t0;
19             xvel = mySpeed;
20          }
21          else if (t1 <= myRunTime) { // contract
22             xlen += mySpeed*(2*myRunTime/2 - t1);
23             xvel = -mySpeed;
24          }
25          // update end point positions and velocities
26          pr0.setPosition (new Point3d (xlen, 0, 0));
27          pr1.setPosition (new Point3d (xlen, 0, -0.5));
28          pr0.setVelocity (new Vector3d (xvel, 0, 0));
29          pr1.setVelocity (new Vector3d (xvel, 0, 0));
30       }
31    }
32
33    // Create and add an output probe to record the tension force of a muscle
34    void addForceProbe (String name, Muscle mus, double stopTime) {
35       try {
36          NumericOutputProbe probe =
37             new NumericOutputProbe (mus, "forceNorm", 0, stopTime, -1);
38          probe.setName (name);
39          addOutputProbe (probe);
40       }
41       catch (Exception e) {
42          e.printStackTrace();
43       }
44    }
Figure 4.14: Tension forces produced by EquilibriumMuscleDemo, displayed by the output probes, as both muscles are extended and then contracted. Left and right images show results with the muscle material property ignoreForceVelocity set to true and false, respectively.

To run this example in ArtiSynth, select All demos > tutorial > EquilibriumMuscleDemo from the Models menu. The model should load and initially appear as in Figure 4.13. As explained above, running the model will cause the muscles to be extended and contracted by moving particles p0r and p1r right and back again. As this happens, the resulting tension forces can be examined by expanding the output probes in the timeline (see section “The Timeline” in the ArtiSynth User Interface Guide). The results are shown in Figure 4.14, with the left and right images showing results with the muscle material property ignoreForceVelocity set to true and false, respectively.

As should be expected, the results for the two muscle implementations are very similar. In both cases, the forces exhibit a local peak near time 0.25 when the muscle lengths are at the maximum of the active force length curve. As time increases, forces then decrease and then rise again as the passive force increases, peaking at time 0.75 when the muscle starts to contract. In the left image, the forces during contraction are symmetrical with those during extension, while in the right image the post-peak forces are more attenuated, because in the latter case the force velocity relationship is enabled as this reduces forces when a muscle is contracting.