Numerical instability of Hill-type muscle models

Abstract

Hill-type muscle models are highly preferred as phenomenological models for musculoskeletal simulation studies despite their introduction almost a century ago. The use of simple Hill-type models in simulations, instead of more recent cross-bridge models, is well justified since computationally 'light-weight'-although less accurate-Hill-type models have great value for large-scale simulations. However, this article aims to invite discussion on numerical instability issues of Hill-type muscle models in simulation studies, which can lead to computational failures and, therefore, cannot be simply dismissed as an inevitable but acceptable consequence of simplification. We will first revisit the basic premises and assumptions on the force-length and force-velocity relationships that Hill-type models are based upon, and their often overlooked but major theoretical limitations. We will then use several simple conceptual simulation studies to discuss how these numerical instability issues can manifest as practical computational problems. Lastly, we will review how such numerical instability issues are dealt with, mostly in an ad hoc fashion, in two main areas of application: musculoskeletal biomechanics and computer animation.

Keywords: Hill-type muscle model; biomechanical models; computer graphics; muscle mechanics; musculoskeletal simulation; numerical instability.

Figure 1.

Force–length–velocity relationship assumed in HMMs and corresponding FL (red) and FV (green) curves. Given the length and velocity of a fully activated muscle, the surface predicts how much force it can produce. Based on FL and FV, four different contractile scenarios can be defined. The grey-shaded area is the only area on the surface predicted by Hill's equation.

Figure 2.

A typical FV and FL curve of HMMs, based on Millard et al. [21]. (a) FV combines concentric (negative velocity) and eccentric (positive velocity) regions, where the concentric part is based on Hill's equation. (b) Total FL (red) is a combination of active (orange) and passive (blue) FL curves. Names of different regions of the curve are based on the shape: ascending (with shoulder) and descending limbs are divided by the plateau formed around L⁰, the optimal length. The passive force that begins to develop after L⁰ cancels out the negative slope of the active FL, but often leaves a dip region where the negative slope still persists.

Figure 3.

Simulations of one hundred serially connected HMM elements. In the first simulation (first row), the initial lengths of HMM elements were all set on the ascending limb (0.85L⁰) with less than 0.1% random non-uniformity across them. The first three columns are snapshots of the states of each element on the FL curve, taken at 0.0, 2.5 and 5.0 s, respectively, where filled red circles indicate lengths of HMM elements in these time points—note that this circle represents one hundred circles superimposed on each other. The second simulation (second row) had the same initial conditions as the first one, with the exception that the initial lengths were set on the descending limb (1.15L⁰). The snapshots and the length profiles clearly indicate a bifurcation of the HMM element lengths due to numerical instability. The rightmost column plots the simulated time–length profiles of all HMM elements for both conditions. Note that the top and bottom plots are at very different length scales. The visualizations at the bottom row are the corresponding snapshots of the volumetric deformation of the muscle based on the length profile simulated in the second simulation scenario. Each HMM element was modelled as an isovolumetric cylinder becoming thicker when shortened and thinner when lengthened.

Figure 4.

Prevalence of numerical instability in an agonist–antagonist system. The effective stiffness of the system was tested for different activation levels (X and Y axes) and total combined lengths (from 1.8L⁰ to 2.8L⁰). In this three-dimensional space, a group of points with negative effective stiffness is visualized as a grey volume, and its projection to the X–Y plane is drawn as a grey polygon. All the other points outside of the grey volume had a positive equilibrium.

Figure 5.

An example simulation of a two-muscle system with negative effective stiffness. Two muscles were activated at 91 and 98% respectively, and their FL relationships (normalized FL curve scaled by activation level) are plotted as blue and red curves in the left-hand panel. Initially, lengths of M₁ and M₂ muscles were set to 1.15L⁰ and 0.95L⁰, respectively, to form a force equilibrium at 0.894$F_{\max }^0$ (red and blue circles). When a small perturbation of 0.1% L⁰ was applied, the system deviated away from the initial condition and reached a stable equilibrium at 0.794$F_{\max }^0$ (red and blue squares).

Figure 6.

Schematic illustration of the residual force enhancement on the descending limb (the dip region) of the FL curve. The colour scheme of the FL curve is identical to the one previously used in figure 2b. A muscle is activated at L_a, a length slightly beyond the plateau and therefore belongs to the dip region where negative stiffness is predicted by the HMM and stretched to a longer length L_b. Whereas the HMM predicts a reduction of active muscle force (purple arrow) along the descending limb, the actual observed active force increases (brown arrow), indicating that the FL relationship does not stay on the FL curve anymore. The enhancement of active muscle force after the stretch (green arrow on the left) compared to the amount predicted by the FL curve is called residual force enhancement. In the actual experiments, the same amount of enhancement (green arrow on the right) will be observed in total force, under the assumption that passive force always follows the passive FL relationship (blue curve).

Figure 7.

Comparison between measured and predicted dynamic force enhancement patterns. Muscle was activated and stretched with a constant velocity at seven different initial and final lengths, which covers the ascending limb, plateau, and descending limb (left) of the FL relationship. Experimentally recorded force profiles are plotted as red curves (middle) and predicted force profiles by a standardized HMM are plotted as blue curves (right).