Each limb segment was modeled as a long thin rigid body of mass m_i and length L_i where the subscript (i = 1, 2, 3) identifies the individual segments. The motion of the center of mass of each body was described by position s→i=(xi,yi,0) angular orientation θ_i, translational velocity v→i=ds→i∕dt, and angular velocity ω_i = dθ_i/dt.

The motions of all bodies were solved simultaneously according to the Newton-Euler method where Newton’s law,

governed the translational acceleration a→i=dv→i∕dt of each limb segment (treated as a point mass) in response to internal joint forces F→Pi and F→Qi, gravitational force G→i and an external damping force F→Di=-kfv→i Likewise Euler’s law,

governed the angular acceleration α_i = dω_i/dt of each limb segment (treated as a rigid body with moment of inertia Ii=miLi2∕12) in response to internal wrenching torques τ_Pi and τ_Qi as well as an external torque τ_Ei and a damping torque τ_Di = − k_τω_i.

Equations (1) and (2) were rearranged as a set of first-order ordinary differential equations and integrated numerically in Matlab. Full details of the forward model are provided in Appendix A. Our extensions to the conventional method are presented in Appendix B. All parameters of the forward model are listed in Table 1.

Muscle was modeled by an active contractile element (CE) in parallel to a passive elastic element (PE) according to the Voigt model (Figure 2A). The net force imparted by the muscle was thus

where a(t) ∈ [0, 1] denotes the instantaneous level of muscle activation, F_max denotes maximal muscle force, θ^ denotes joint angle and ω^ denotes joint opening velocity. Muscle length was treated as a linear function of joint angle, as justified by cadaver studies (Grieve et al., 1978). Notice we use hat notation to distinguish joint angle θ^ and joint opening velocity ω^ from limb segment orientation θ and limb segment angular velocity ω.

The CE force-length relationship (blue line in Figure 2B) was modeled by a Gaussian curve,

centered on resting muscle length θ^0 as is conventional (see Winters and Stark, 1987; Zajac, 1989). The PE force-length relationship (green line in Figure 2B) was approximated by a hyperbolic curve,

where θ^min and θ^max specified the limits of joint movement and the constant k_pe specified the slope. The CE force-velocity relationship (Figure 2C) was approximated by an exponential hyperbolic tangent,

with a single slope parameter, k_v, rather than piecewise hyperbolics after Hill (1938). See Table 2 for a list of all parameters related to muscle dynamics.

We analyzed the angular motion of an isolated limb segment actuated by a pair of opposing muscles to gain insights into the effects of co-contraction. The muscles imposed opposing torques,

and

on the limb segment where F denotes contractile muscle force (3) and mom denotes the moment arm of the muscle. These torques were substituted into Euler’s law (2) and integrated numerically as a set of first-order ordinary differential equations.

The damping effects of co-contraction were investigated by comparing the motions of the isolated limb under differing conditions of muscle co-activation but identical net joint torques. Phase portraits of the motion (θ_a versus ω_a) were computed for conditions of nil muscle activation (a_a = a_b = 0), 25% muscle activation (a_a = a_b = 0.25), and 50% muscle activation (a_a = a_b = 0.5) where the latter represents the upper limit of co-contraction in physiological conditions. Each phase portrait portrays all possible motions of the limb segment for a given set of muscle activations. The flow of the vector field in the phase portrait can be characterized by the nullclines of the dynamic variables which corresponds to those cross sections through phase space where growth of one or more variables is zero. Equilibrium points occur where the nullclines intersect and the growth of all variables is zero. Linearizing the flow around equilibrium points allows their stability to be quantified by their eigenvalues. The real part of the eigenvalue quantifies the rate of convergence to the equilibrium point (damping) whereas the imaginary part of the eigenvalue quantifies the rotational component (oscillation frequency) of the vector field around the equilibrium point. When a parameter change causes the real part of the eigenvalue to cross zero from below, stability of the equilibrium point is lost as fluctuations grow exponentially. This is known as a local bifurcation of a continuous dynamical system and is classified as either a saddle-node, transcritical, pitchfork, or Hopf bifurcation according to the nature of the ensuing flow (see Strogatz, 2000; Breakspear and Jirsa, 2007, for reviews).

Preliminary analysis of these nullclines also suggested that gross asymmetries in the force-length properties of the opponent muscles may cause the equilibrium position of the joint to bifurcate into a co-existing pair of non-unique equilibrium positions under high levels of co-contraction. Phase portraits were thus computed for opponent muscles in which the resting length of the CE force-length property θ₀ had been shifted away from the midpoint of the muscle’s range of extension by the arbitrary amount of 1 rad (θ^0a=θ0b=π-1) under various conditions of balanced co-contraction (a_a = a_b) and imbalanced co-contraction (a_a varied while a_b held fixed). Parameter space was explored by numerical continuation of the observed bifurcation points was performed using the CL_MATCONT numerical toolkit (Dhooge et al., 2003).

The effect of co-contraction mediated muscle damping was verified in the full biomechanical limb by comparing the motion of identical limbs under conditions of medium (near 20%), strong (near 50%), and extreme (near 80%) co-contraction (Table 3) with identical net muscle activations (a→A-a→B) across conditions. Extreme co-contraction represents the theoretical limit of co-contraction in the mathematical sense and does not occur in nature. Muscle activations were held constant for the duration of the simulation (10 s) with the limb initially hanging at rest from the base pivot. All other parameters were held fixed (m_i = 1, L_i = 1, k_f = 0, k_τ = 0, k_s = 1, k_d = 1, G_ix = 0, G_iy = − 9.81, F_max,1 = 4000, F_max,2 = 2000, F_max,3 = 1000, k_pe = 0.1, k_l = − π²/ln(0.1), k_v = 0.2, ɸ_a = − 0.2, ɸ_b = 0.2, ɸ_c = 0.2, ɸ_d = − 0.2, θ_min = 0.3π, θ_max = 1.6π, θ₀ = 0.95π).

The effect of co-contraction mediated multistability in asymmetric muscles was verified in the full biomechanical limb model by comparing the final postures adopted (at t = 30 s) by n = 200 simulation runs of identical limbs undergoing medium, strong, and extreme co-contraction from random initial postures (always at rest). The same set of initial postures were used in both co-conditions and muscle asymmetries were imposed by shifting the CE resting length of opposing muscles to θ^0=π-1. All other parameters were the same as those in Experiment 1.

Analysis of the isolated limb segment revealed co-contraction modulated muscle damping effects emerged natively from the muscle model. Furthermore, a wide range of muscle activations were found to yield bistable equilibria when the opponent muscles were configured with asymmetric CE force-length properties.

Figure 3 reveals the differing levels of damping exhibited by the isolated limb segment under conditions of 0, 25, and 50% co-contraction. Damping was entirely absent under 0% co-contraction (Figure 3A) as seen by the closed orbits around the equilibrium position. The lack of damping was confirmed quantitatively by noticing that the real part of the eigenvalues of the equilibrium position were zero (λ_1,2 = 0 ± 0.174i). By comparison, progressively more damping was observed in the 25% co-contraction (Figure 3B) and 50% co-contraction (Figure 3C) conditions, as quantified by progressively larger negative values in the real parts of the eigenvalues (λ_1,2 = − 0.0657 ± 0.161i and λ_1,2 = − 0.1314 ± 0.114i respectively).

Figure 4A reveals the existence of bistable equilibrium positions in the same joint when the muscles were configured with asymmetric CE force-length properties (θ_a0 = θ_b0 = π − 1) and subjected to 50% co-contraction. Since the nullcline of θ^ is always ω^=0 in the present model, we need only consider the nullcline of ω^ to understand the motion of the joint. Figure 4B shows the nullclines for the asymmetrically actuated limb under differing levels of balanced co-contraction (a_a = a_b = 0, …, 1). It reveals a supercritical pitchfork bifurcation in the equilibrium positions as co-contraction exceeds the critical value of a_a = a_b = 0.17 (indicated by the branch point BP). The pitchfork bifurcation occurs when the slope of the nullcline at the central stable equilibrium point changes sign resulting in a nullcline with three distinct zero crossings. In the supercritical case, these three zero crossings correspond to a central unstable equilibrium point flanked by a pair of stable equilibrium points.

Similarly, Figures 4D–F describe the motion of the asymmetrically actuated limb when a_b = 0.5 is held fixed while a_b is manipulated. Here we observe a saddle-node bifurcation that yields bistable equilibrium positions when muscle activations are within the critical range 0.338 < a_a < 0.998. The saddle-node (or fold) bifurcation occurs when the turning points in the nullcline fold back far enough to support multiple zero crossings, specifically, one unstable equilibrium flanked by two stable equilibria.

Figure 5A shows the bistability map for the asymmetric joint with θ₀ = π − 1. It reveals the extent of muscle activations that yield bistable equilibrium positions in this case. The cusp point (CP) corresponds to the branch point in Figure 4C and it marks the furthest extent of the bistable region. Its position (always on the line a_a = a_b) depends upon the CE resting length parameter θ₀. Figure 5B plots the migration of the cusp point as θ₀ is manipulated. Here we see that the bistability emerges at 100% co-contraction levels when θ₀ = 3.00 rad (point A) while it emerges at 50% co-contraction levels when θ₀ = 2.86 rad (point B). Bistability therefore emerges in the biological operating range of co-contraction when the CE resting lengths of the opposing muscles are shifted away from the muscle midpoint by as little as (π − 2.86) radians (16.1°).

In comparison, bistability emerges at 17.1% co-contraction levels when θ₀ = π − 1 (point C) which is the value used in our simulations. This happens to be close to the maximal extent of bistability which emerges at only 15.4% co-contraction levels when θ₀ = 1.68 rad (point D) and corresponds to shifting θ₀ by (π − 1.68) radians (83.7°) from the muscle midpoint – which is an extreme shift.

Co-contraction modulated damping effects were confirmed in the full biomechanical limb. Figures 6A–C show the trajectories of the limb under conditions of medium, strong, and extreme co-contraction. Figure 6D tracks the vertical position of the limb tip for each condition where limb damping is observed to increase with higher levels of co-contraction. The limb converged to the same equilibrium position in all conditions, confirming that multistability did not apply when the opponent muscles were symmetrically matched.

Co-contraction modulated multistability was confirmed in the full biomechanical limb with asymmetric CE force-length properties. Figure 7A shows the random initial limb postures used in all conditions. Figures 7B–D show the final limb postures for medium, strong, and extreme co-contraction respectively. All limbs undergoing medium co-contraction converged to the same posture (Figure 7B) and are indistinguishable. In contrast, limbs undergoing extreme co-contraction (Figure 7D) converged to one of six distinct stable postures which clearly demonstrated that all joints were operating in the bistable regime, whereas limbs undergoing strong co-contraction (Figure 7C) converged to one of only two distinct equilibrium postures. Nonetheless, all joint angles in this condition diverged noticeably from those observed under medium co-contraction suggesting all were operating in the bistable regime even though only some had converged to distinct equilibrium positions.

Not all of the possible combinations of (bi)stable joint solutions yield stable limb postures in Figures 7C,D. This is particularly noticeable in Figure 7D where two of the eight possible postural combinations of bistable joints are rendered unstable by the action of gravity. Critical slowing of the limb is evident in the vicinity of these “missing” postures (see Movie S2 in Supplementary Material) confirming that the system is close to a stable regime even though stability is not achieved in this case.

By dynamical system theory the limb joint is guaranteed to be monostable when the velocity nullcline is monotonically decreasing and thus has only one zero-crossing (as is the case in Figure 3). Here, we derive an analytical expression for the equilibrium position of the isolated limb segment for the special case of identical antagonist muscles. We then use that expression to prove that monostability is guaranteed in the special case. This proof does not preclude the possibility that monostability may also hold when antagonistic muscles are similar but not identical.

The equilibrium position of the limb segment is given by the intersection of the nullclines of ω^ and θ^ but since the nullcline of θ^ is always ω^=0 for our system we need only consider the zero crossings of the nullcline of ω^. Solving Euler’s law of motion for the case where dω^a∕dt=0 and substituting muscle equations (3–6) allows the nullcline to be expressed as

where

are defined for brevity. Observe that ω^a=-ω^b=0 at the zero crossings of (9) and all force-velocity relations have f_v(0) = 1 by definition thus the zero crossings of the nullcline are obtained by solving the expression

where θ_a + θ_b = 2π. We now use this expression to prove the following theorem.

Theorem. The nullcline (9) of an isolated limb segment actuated by an identical pair of opposing Voigt muscles always has exactly one root (ω^=0) for any given pair of fixed muscle activations, a_a and a_b. Thus monostability of the limb segment is guaranteed for all possible muscle activations.

Proof. The roots of nullcline (9) are given by equation (10). When the opponent muscles (denoted “a” and “b”) have identical muscuoloskeletal properties (mom_a = mom_b, F_max,a = F_max,b, k_la = k_lb, k_pea = k_peb, θ^a+θ^b=2π,θ^0a+θ^0b=2π,θ^min,a+θ^max,b=2π,θ^max,a+θ^min,b=2π) then equation (10) reduces to

where flb(θ^b) can be expressed in terms of fla(θ^a) by substituting θ^b=2π-θ^a,θ^0b=2π-θ^0a,andklb=kla as follows,

Similarly, fpeb(θ^b) can be expressed in terms of θ^a by substituting θ^b=2π-θ^a,θ^min,b=2π-θ^max,a,θ^max,b=2π-θ^min,a,andkpe=kpea=kpeb as follows,

Substituting (12) and (13) into (11) yields

where

and

are defined for convenience. From this point onward we omit the muscle subscript (a) from the notation for brevity. Notice that G(θ^) corresponds to a Gaussian curve centered on θ^0 with a peak amplitude of (a_a − a_b)/k_pe. Notice also that H(θ^) corresponds to the difference of two hyperbolics which asymptotes vertically to +∞ at θ^=θ^min and to −∞ at θ^=θ^max with a zero crossing at θ^0=(0^min+θ^max)∕2. Functions G(θ^) and H(θ^) are both continuous on the interval (θ^min,θ^max).

We argue geometrically that functions G(θ^) and H(θ^) always intersect exactly once on the interval (θ^min,θ^max) for any value of (a_a − a_b) and therefore equation (11) always has exactly one solution. It follows from the definition of (11) that nullcline (9) always has exactly one root when the opposing muscles are symmetric.

While co-contraction is known to modulate both muscle damping and musculotendon stiffness, it has not previously been implicated with bistable equilibria to the best of our knowledge. Indeed, mainstream theoretical accounts of biological motor control (e.g. Feldman, 1966; Feldman and Levin, 1995) typically assume that co-activations of antagonist muscles implicitly translate into monostable equilibrium positions. However our findings suggest that such an assumption may not always be justified. Existing theoretical accounts may therefore need to accommodate the additional complexities of either controlling or avoiding bistability in the muscle apparatus to achieve unambiguous forward control of the joint.

Whether co-contraction mediated bistability occurs in nature remains an open question. Our analysis suggests that bistable postures can emerge at biologically plausible levels of co-contraction (50% maximal isometric force) when the CE resting length parameters are offset from the muscle midpoint by as little as 16.1°. This critical offset value corresponds to a peak-to-peak discrepancy in the force-length curves of the antagonist muscles of 32.2° (twice the offset). This critical value happens to be exceeded by the peak-to-peak discrepancies in human elbow (40°), knee (40°), and ankle (60°) reported by Winters and Stark (1985). So it is not unreasonable to consider that bistable antagonist muscle systems may indeed exist in nature.

The increasing use of biomimetic actuators in robotic systems (e.g. Ayers et al., 2002; Safak and Adams, 2002) complements the neurorobotic doctrine that the brain cannot be studied separately from the body (Chiel and Beer, 1997). The central nervous system’s ability to use muscle co-contraction to actively control limb damping highlights the tight integration between the functionality of the motor control system and the contractile dynamics of muscle tissue. However not all muscle properties have equivalent functional relevance to the motor control system so only the most relevant muscle properties need be incorporated into biomimetic actuators.

Our biomechanical model demonstrates that simplified actuators with idealized forms of CE force-length-velocity properties are sufficient to permit active control of joint damping through co-contraction. More specifically, it is crucial that the force-velocity property of the actuator be modulated by its activation level for co-contraction modulated damping effects to occur. Musculotendon series elasticity, for example, is not necessary for this purpose. As we have shown, care must be taken with the force-length property of the actuator to avoid gross non-linearities in the dynamics of opposing actuators which can introduce non-unique equilibria into the forward dynamics that merely complicate the problem of joint control. Thankfully, numerical stability analysis allows the monostable operating limits of antagonist actuators to be ascertained from the contractile dynamics of a single biomimetic actuator.

The contractile dynamics of the present biomechanical model were simplified by the use of the Voigt muscle model which lacks a series elastic element. This simplification made it possible to derive an analytical solution to the equilibrium position of antagonist muscles and to prove that monostability is guaranteed for identical muscles. However musculotendon series elasticity is known to have a significant impact on muscle stiffness during rapid movement (Winters and Stark, 1985; Winters et al., 1988). The physiological accuracy of the present model is therefore limited to low and moderate speed movements where dynamic changes in muscle stiffness are negligible.