The Hill Muscle Model

This is a walkthrough of the Hill-type muscle model based on Geyer & Herr (2010) and Wang et al. (2012). I wrote this because I couldn't find the full model described in a single source while implementing it.

Elements

The model consists of three elements:

A contractile element (CE) that represents the muscle fibers that contract based on the muscle activation state.
A serial elastic element (SEE) that represents the tendons that connect the muscle to the bones.
A parallel elastic element (PEE) that represents the passive elastic material surrounding muscle fibers.

Forces

$F_{max}$: maximum isometric force. Optimized constant.
$F_{CE}$: contractile element force. Eq. (1).
$F_{SEE}$: serial elastic element force. Eq. (2).
$F_{PEE}$: parallel elastic element force. Eq. (3).
$F_M$: total muscle force. Eq. (4).

Force relationships

$$F_M = F_{CE} + F_{PEE} = F_{SEE} \tag{4}$$

Contractile element

From Wang et al. 2012 and confirmed by author:

$$L_{CE}^{norm} = \frac{L_{CE}}{L_{CE}^{opt}}$$

$$V_{CE}^{norm} = \frac{V_{CE}}{L_{CE}^{opt}}$$

Force-length relationship [Geyer et al. 2003]

$$f_L(L_{CE}^{norm}) = \exp \left[ c \left\lvert\frac{L_{CE}^{norm} - 1}{w}\right\rvert^3 \right]$$

$w$: width of the bell-shaped $f_L$ curve.
$c$: $\ln(0.05)$ fulfilling $f_L(1 \pm w) = 0.05$

Force-velocity relationship [Geyer et al. 2003]

$$f_V(V_{CE}^{norm}) = \begin{cases} \dfrac{V_{max} - V_{CE}^{norm}}{V_{max} + K V_{CE}^{norm}} & \text{if } V_{CE}^{norm} < 0 \ N + (N - 1) \dfrac{V_{max} + V_{CE}^{norm}}{7.56 K V_{CE}^{norm} - V_{max}} & \text{if } V_{CE}^{norm} \geq 0 \end{cases}$$

$V_{CE} < 0$ follows the Hill equation (Hill 1938) for muscle shortening where $V_{max} < 0$ is the maximum contraction velocity.
$V_{CE} \geq 0$ muscle lengthening is characterized by an equation based on Aubert (1956).
$V_{max} = -12 , l_{opt} s^{-1}$: maximum shortening velocity.
$K$: curvature constant.
$N$: dimensionless force $\frac{F_M}{F_{max}}$ reached at lengthening velocity $V_{CE} = -V_{max}$.

Relationship (cannot compute directly for $V_{CE}$ because it's recursive):

$$f_V(V_{CE}^{norm}) = \frac{F_{SEE} - F_{PEE}}{a F_{max} f_L(L_{CE}^{norm})}$$

Force-velocity inverse

$$f_V^{-1}(y) = \begin{cases} \dfrac{V_{max} (1 - y)}{1 + y K} & \text{if } y < 1 \ \dfrac{V_{max} (y - 1)}{7.56 K (y - N) - N + 1} & \text{if } y \geq 1 \end{cases}$$

Contractile-element velocity

$V_{CE}$ is derived from Equations (1) and (4), by inverting the force-velocity relation $f_V$:

$$V_{CE} = f_V^{-1} \left[\frac{F_{SEE} - F_{PEE}}{a F_{max} f_L(L_{CE}^{norm})}\right] \tag{6}$$

Because:

$$\frac{F_{CE}}{\frac{F_{CE}}{f_V(V_{CE}^{norm})}} = f_V(V_{CE}^{norm})$$

And the inverse function relationship $x = f^{-1}(f(x)) = f(f^{-1}(x))$, thus:

$$V_{CE} = f_V^{-1}(f_V(V_{CE}))$$

Contractile-element force

$$F_{CE} = a F_{max} f_L(L_{CE}^{norm}) f_V(V_{CE}^{norm}) \tag{1}$$

Contractile element force as a function of normalized length and velocity.

Serial elastic element

The passive forces produced by the elastic elements, $F_{PEE}$ and $F_{SEE}$, are modeled as non-linear springs based on their lengths:

$$F_{SEE} = f_{SEE}(L_M - L_{CE}) = f_{SEE}(L_{SEE}) \tag{2}$$

$f_{SEE}$ is described in Geyer et al. 2003:

$$f_{SEE} = \begin{cases} \left(\dfrac{\epsilon}{\epsilon_{ref}}\right)^2 & \text{if } \epsilon > 0 \ 0 & \text{if } \epsilon \leq 0 \end{cases}$$

$\epsilon$: tendon strain. $\epsilon = \frac{L_{SEE} - L_{SEE}^{rest}}{L_{SEE}^{rest}}$
$L_{SEE}^{rest}$: tendon's resting length
$\epsilon_{ref}$: tendon reference strain with $f_{SEE}(\epsilon_{ref}) = 1$

Parallel elastic element

$$F_{PEE} = f_{PEE}(L_{CE}^{norm}) \tag{3}$$

$f_{PEE}$ isn't described in Geyer et al. 2003 but is described in Geyer et al. 2010:

$$f_{PEE} = F_{max} \left(\frac{L_{CE}^{norm} - L_{CE}^{opt}}{L_{CE}^{opt} \epsilon_{PE}}\right)^2 f_V(V_{CE}^{norm})$$

$\epsilon_{PE} = w$ according to Geyer et al. 2010.

Parallel elastic element force as a function of normalized length and velocity.

Lengths

$L_M$: total muscle length. Eq. (8).
$L_{CE}$: contractile element length. Initialized to $L_{CE}^{opt}$ then integrated from $V_{CE}$.
$L_{SEE}$: serial element length. Derived from relationship.
$L_M^{rest}$: muscle rest length. Derived from relationship.
$L_{SEE}^{slack}$: tendon slack length. Optimized parameter.
$L_{CE}^{opt}$: optimal contractile element length. Optimized parameter.

$$L_M = L_{SEE} + L_{CE}$$

$$L_M^{rest} = L_{SEE}^{slack} + L_{CE}^{opt} \tag{5}$$

$$L_M = \sum_{i=1}^{n-1} |s_i| \tag{8}$$

Muscle segment length:

$$s_i = P_{i+1}^W - P_i^W$$

Contractile element length

The muscle state parameter $L_{CE}$ is updated each time step through integration of the contraction velocity $V_{CE}$.

Be sure to use the velocity, then update it for Euler integration.

$$L_{CE} = \int V_{CE} , dt \tag{6}$$

$$L_{CE} = L_{CE} + V_{CE} , dt$$

Activation

$a$: muscle activation
$u$: neural activation signal
$c_a$: constant activation and deactivation rate. Default: $c_a = 100 s^{-1}$ or $0.01$.

$$\frac{\partial a}{\partial t} = c_a (u - a) \tag{7}$$

Force Application

The total muscle length $L_M$ with activation state $a$ to compute the scalar muscle contraction force $F_M$.

$k$: a joint spanned by the muscle.
$\tau_k$: the torque generated in the direction of the moment arm $r_k$.
$r_k$: moment arm corresponding to the cross-product between the direction of the segment that crosses the joint $s_c$ and a vector from the joint center $j_k$ to any point on the segment $s_c$, such as $P_c^W$.
$s_c$: muscle segment.
$j_k$: joint center.

$$\tau_k = F_M |r_k| \tag{9}$$

$$r_k = (P_{i+1}^W - P_i^W) \times \frac{s_c}{|s_c|}$$

Force relationships

Torque: $T = \frac{\partial L}{\partial t}$
Net torque: $T_{net} = I \cdot a$
Moment of inertia: $I = \frac{L}{w}$
Angular momentum: $L = I \cdot w$
Force: $F = m \cdot a$
Angular acceleration: $a = \frac{T}{I}$

References

Geyer, H. & Herr, H. (2010) — A Muscle-Reflex Model That Encodes Principles of Legged Mechanics Produces Human Walking Dynamics and Muscle Activities
Wang, J.M. et al. (2012) — Optimizing Locomotion Controllers Using Biologically-Based Actuators and Objectives
Geyer, H. et al. (2003) — Positive force feedback in bouncing gaits?
Hill, A.V. (1938) — The Heat of Shortening and the Dynamic Constants of Muscle