4 Mechanical Models II 4.4 Custom rendering 4.6 Distance Grids and Components

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 , length velocity $\dot{l}$ , and excitation signal , 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 and , 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, and 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, by the excitation property of the Muscle, and , and $f_{max}$ by the following properties of SimpleAxialMuscle:

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

SimpleAxialMuscles can be created with the following constructors:

⬇

SimpleAxialMuscle()

SimpleAxialMuscle (stiffness, damping, fmax)

where stiffness, damping, and fmax specify , 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 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
	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 , 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
	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 , 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 [17]. 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

. 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
	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 , 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 [4]. 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

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
	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 , 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 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 :

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

As the muscle moves, it is assumed that the height 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
	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}$ , and $\phi_{o}$ , where $F_{o}$ , $l_{o}$ and 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 .

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

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 [6] and described in [14].

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

vs.

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}$ , , and $\phi_{o}$ , and other properties are set to their defaults.

4.5.4.2 Thelen2003AxialMuscle

Thelen2003AxialMuscle implements the Thelen2003Muscle model supplied by OpenSim [6] and introduced in [24].

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 [24], 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 [24], 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 [24]. 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 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}$ , , 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 is the tendon slack length. $F_{o}$ and are specified by the following properties of Millard2012AxialTendon:

	Property	Description	Default
$F_{o}$	maxIsoForce	maximum isometric force	1000
	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 , 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 is the tendon slack length. $F_{o}$ , , 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
	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 , and other properties are set to their defaults.

4.5.5.3 Blankevoort1991AxialLigament

Blankevoort1991AxialLigament implements the Blankevoort1991Ligament model supplied by OpenSim [6] and described in [3, 21].

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}$ , , and 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
	linearStiffness	maximum stiffness of force-length curve	100
	damping	damping parameter	0.003

Blankevoort1991Ligaments can be created with the constructors

⬇

Blankevoort1991Ligament()

Blankevoort1991Ligament (stiffness, slackLen, damping)

where stiffness, slackLen and damping specify , $l_{o}$ and , 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 [13]). 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 , 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

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;

9 // initial total length of the muscles:

10 private double len0 = 0.25 + myTendonSlackLen;

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);

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);

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);

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);

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);

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);

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);

53 pr1 = new Particle ("pr1", 1.0, len0, 0, -0.5); // right end point

54 pr1.setDynamic (false);

55 mech.addParticle (pr1);

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);

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));

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);

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);

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

7 double myRunTime = 1.5; // total extensions/contraction time

8 double mySpeed = 1.0; // speed of the end point motion

10 public LengthController() {

11 // need null args constructor if this model is read from a file

12 }

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 }

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.

4.4 Custom rendering Bibliography 4.6 Distance Grids and Components

Generated on Tue Mar 26 13:12:56 2024 by LaTeXML [LOGO]