Muscle co-contraction modulates damping and joint stability in a three-link biomechanical limb.

Computational models of neuromotor control require forward models of limb movement that can replicate the natural relationships between muscle activation and joint dynamics without the burdens of excessive anatomical detail. We present a model of a three-link biomechanical limb that emphasizes the dynamics of limb movement within a simplified two-dimensional framework. Muscle co-contraction effects were incorporated into the model by flanking each joint with a pair of antagonist muscles that may be activated independently. Muscle co-contraction is known to alter the damping and stiffness of limb joints without altering net joint torque. Idealized muscle actuators were implemented using the Voigt muscle model which incorporates the parallel elasticity of muscle and tendon but omits series elasticity. The natural force-length-velocity relationships of contractile muscle tissue were incorporated into the actuators using ideal mathematical forms. Numerical stability analysis confirmed that co-contraction of these simplified actuators increased damping in the biomechanical limb consistent with observations of human motor control. Dynamic changes in joint stiffness were excluded by the omission of series elasticity. The analysis also revealed the unexpected finding that distinct stable (bistable) equilibrium positions can co-exist under identical levels of muscle co-contraction. We map the conditions under which bistability arises and prove analytically that monostability (equifinality) is guaranteed when the antagonist muscles are identical. Lastly we verify these analytic findings in the full biomechanical limb model.

Introduction

Forward models of musculoskeletal dynamics replicate the natural relationship between muscle contraction and limb movement. Forward models are necessary for numerical optimization of motor control strategies (Todorov, 2004) and for exercising theoretical models of neuromotor control (e.g., Conforto et al., 2009; Harischandra et al., 2010). At least one commercial package (LifeMOD) is currently available for constructing anatomically precise forward models of specific body parts and such models are routinely applied to the inverse problem of estimating isolated muscle and joint forces from observed limb movements (Erdemir et al., 2007). However this level of anatomical detail is excessive when exploring basic theoretical principles of neuromotor control. In such cases the dynamical character of limb movement is of primary importance and the use of simplified anatomical models is justified. We present a planar model of an idealized biomechanical limb (Figure 1) that evokes naturalistic limb movements in response to muscle contractions and co-contractions yet remains amenable to customization by individual researchers.

Figure 1

Schematic of the simulated biomechanical limb where each joint is flanked by an opposing pair of muscle actuators (arrows). Each limb segment was modeled as a long, thin, rigid body of length L and mass m.

Co-contraction (the simultaneous contraction of antagonist muscles) is known to stabilize limb movements (Milner and Cloutier, 1995; Milner, 2002; Zakotnik et al., 2006; Lametti et al., 2007) and is regarded as a distinct component of the motor command in many theoretical models of motor control (Feldman and Levin, 1995; Bhushan and Shadmehr, 1999; Gribble and Ostry, 1999; Todorov, 2000; Gribble, 2003; Neilson and Neilson, 2005). However traditional planar limb models typically lump opposing muscles into unitary joint actuators that do not accommodate co-contraction. The present model overcomes this limitation by actuating each joint with an antagonistic pair of muscle actuators that may be activated independently.

Co-contraction alters the biomechanical operating ranges of muscle and tendon by increasing both muscle damping (viscosity) and musculotendon stiffness (elasticity). Winters and Stark (1985, 1988) show that dynamic limb impedance (resistance to perturbation) is accurately predicted by an eighth-order model of antagonist Hill-based muscles. Such models (after Hill, 1938) include a series elastic (SE) element that represents the passive elasticity of muscle tissue and tendon (see Winters and Stark, 1987; Zajac and Winters, 1990; Pandy, 2001; Winter, 2005; Erdemir et al., 2007, for reviews). In these models the series elastic stiffness increases (becomes less compliant) non-linearly with stretch which results in increased joint stiffness under co-contraction (Winters et al., 1988).

Despite this, musculotendon series elasticity is typically an order of magnitude stiffer than contractile muscle tissue and we treat it as negligible in the present model. This simplification permits insights into the dynamics through a formal stability analysis, although it comes at the loss of co-contraction mediated changes in musculotendon stiffness. Nonetheless, the model retains sufficient kinematic realism to make it a useful platform for exercising neuromotor models of low to moderate speed movements where series elasticity has little impact.

In the present paper we derive the full biomechanical limb model followed by a numerical stability analysis of an isolated pair of antagonistic muscles using the method of numerical continuation. The stability analysis reveals the dynamic character of opponent muscles and illuminates the effects of co-contraction from a dynamical systems viewpoint. In particular, it reveals that strongly co-contracting muscles can exhibit multiple distinct stable (multistable) equilibrium positions depending on the biomechanical properties of the muscles involved. In biomechanical parlance, monostability satisfies equifinality whereby a limb always returns to same equilibrium position following a perturbation (Bizzi et al., 1978; Kelso and Holt, 1980; Schmidt et al., 1986; Feldman and Latash, 2005). We analytically derive the nullclines of the opponent muscle system and formally prove that monostability is guaranteed for the special case of opponent muscles with identical properties. Two numerical experiments are presented which validate these findings in the full biomechanical limb. Experiment 1 demonstrates the effects of co-contraction on limb damping and Experiment 2 demonstrates the emergence of multiple equilibrium postures with strong co-contraction.

Methods

Skeletal and musculotendon kinematics were modeled as a series of independent transforms following Zajac (1989) and Pandy (2001). The forward dynamics were implemented using the Newton-Euler method following Otten (2003). The muscle dynamics were implemented using the Voigt muscle model (Figure 2A) where the force-length and force-velocity properties of muscle tissue were approximated by simple mathematical forms (Figures 2B,C) that captured the general shape of curves reported in the physiological literature.

Figure 2

(A) Voigt model of muscle elasticity. Contractile element (CE) represents the lumped contraction forces of the muscle sarcomeres. Parallel elastic element (PE) represents the lumped elastic forces of muscle tissue. Series elasticity of muscle and tendon is assumed to be negligible and muscle length is considered a linear function of joint angle. (B) Force-length relationships for CE (blue line) with θ = 60°, θ0 = 180°, θ = 300°, k = − (θ − θ0)2/ln(0.1), and PE (green line) with k = 0.1. Vertical axis is normalized with respect to F. (C) Force-velocity relationship for CE with k = 1. Maximal muscle force asymptotes at 1.313 F.

Contemporary models of forward dynamics are usually derived from the Lagrangian formulation of the overall energy in the system implemented using kinematic chains of rigid links (see Craig, 1989; Spong et al., 2006). However the equations of motion are typically non-trivial and are often derived using a symbolic algebra solver in practice (see Westervelt et al., 2007, for a modern tutorial on this approach). The Newton-Euler method, on the other hand, is the more straightforward approach for modeling small systems such as ours (Otten, 2003). Although it does suffer from the problem of numerical drift whereby rounding errors accumulate over time resulting in a slow dislocation of the limb joints. We therefore extended the method to incorporate spring-like binding forces within the limb joints to constrain numerical drift and maintain the geometric integrity of the limb over the long term (see Appendix B). It should be noted that contemporary methods using kinematic chains do not suffer from numerical drift so our solution does not apply to those methods.

Forward dynamics of skeletal mechanics

Each limb segment was modeled as a long thin rigid body of mass m and length L 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 angular orientation θ, translational velocity and angular velocity ω = dθ/dt.

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

governed the translational acceleration of each limb segment (treated as a point mass) in response to internal joint forces and , gravitational force and an external damping force Likewise Euler’s law,

governed the angular acceleration α = dω/dt of each limb segment (treated as a rigid body with moment of inertia ) in response to internal wrenching torques τ and τ as well as an external torque τ and a damping torque τ = − kτω.

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.

Table 1

Description of all parameters related to the forward dynamics.

Parameter	Description
Li	Length of ith limb segment (m)
mi	Mass of ith limb segment (kg)
I_i=m_iLi2∕12	Moment of inertia (kg m2)
s→_i={s_ix,s_iy,0}	Position of ith segment (m)
v→_i={v_ix,v_iy,0}	Velocity (m s−1)
a→_i={a_ix,a_iy,0}	Acceleration (m s−2)
θi	Orientation of ith segment (rad)
ωi = dθi/dt	Angular velocity (rad s−1)
αi = dωi/dt	Angular acceleration (rad s−2)
r→_i={r_ix,r_iy,0}	Radial arm of ith segment (m)
rix = Li cos(θi)/2	Radial arm x-component (m)
riy = Li sin(θi)/2	Radial arm y-component (m)
ɸia, ɸib, ɸic, ɸid	Angles of muscle insertions (rad)
momia = Li sin (ɸia)/2	Moment arm of muscle a (m)
F→_Pi={F_Pix,F_Piy,0}	Joint force at proximal tip (N)
F→_Qi={F_Qix,F_Qiy,0}	Joint force at distal tip (N)
F→_Di=-k_fv→_i	Translational damping force (N)
G→_i={G_ix,G_iy,0}	Gravitational force (N)
J→_ij={J_ijx,J_ijy,0}	Joint-spring force (N)
τPi = FPimomi	Proximal joint torque (N rad)
τQi = − FQimomi	Distal joint torque (N rad)
τEi	External torque (N rad)
τDi = − kτωi	Damping torque (N rad)
kf	Translational damping constant
kτ	Angular damping constant
ks	Joint-spring stiffness constant
kd	Joint-spring damping constant

Muscle contraction dynamics

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, Fmax 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 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 and specified the limits of joint movement and the constant k specified the slope. The CE force-velocity relationship (Figure 2C) was approximated by an exponential hyperbolic tangent,

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

Table 2

Description of parameters related to muscle dynamics.

Parameter	Description
θ^	Joint angle (rad)
θ^_min	Lower limit of joint angle (rad)
θ^_max	Upper limit of joint angle (rad)
θ^_0	Resting length of CE (rad)
ω^	Joint opening velocity (rad s−1)
Fm	Instantaneous muscle force (N)
Fce	Force imparted by CE (N)
Fpe	Force imparted by PE (N)
Fmax	Maximal muscle force (N)
a(t)	Instantaneous muscle activation
f_l(θ^)	CE force-length relation
f_v(ω^)	CE force-velocity relation
f_pe(θ^)	PE force-length relation
kl	Slope of CE force-length relation
kv	Slope of CE force-velocity relation
kpe	Slope of PE force-length relation

Stability analysis of opposing muscles

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 (θ versus ω) were computed for conditions of nil muscle activation (a = a = 0), 25% muscle activation (a = a = 0.25), and 50% muscle activation (a = a = 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 θ0 had been shifted away from the midpoint of the muscle’s range of extension by the arbitrary amount of 1 rad under various conditions of balanced co-contraction (a = a) and imbalanced co-contraction (a varied while a 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).

Numerical experiment 1: Effect of co-contraction on damping

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 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 = 1, L = 1, k = 0, k = 0, k = 1, k = 1, G = 0, G = − 9.81, Fmax,1 = 4000, Fmax,2 = 2000, Fmax,3 = 1000, k = 0.1, k = − π2/ln(0.1), k = 0.2, ɸ = − 0.2, ɸ = 0.2, ɸ = 0.2, ɸ = − 0.2, θmin = 0.3π, θmax = 1.6π, θ0 = 0.95π).

Table 3

Muscle activation levels for the conditions of .

Condition	Muscle activation levels
Medium	a→_A(t)={0.3,0.1,0.3},a→_B(t)={0.1,0.3,0.1}
Strong	a→_A(t)={0.6,0.4,0.6},a→_B(t)={0.4,0.6,0.4}
Extreme	a→_A(t)={0.9,0.7,0.9},a→_B(t)={0.7,0.9,0.7}

The net muscle activation .

Numerical experiment 2: Effect of asymmetric muscles

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 All other parameters were the same as those in Experiment 1.

Results

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 3

Phase portraits of the motion of an isolated limb segment actuated by an opposing pair of symmetric muscles under conditions of . The latter approximates the biological limit of muscle co-contraction in nature. All other musculoskeletal parameters are identical in all three conditions (m = 5, L = 1, ɸ = 0.2, Fmax = 1, k = −π2/In(0.1), k = 2, k = 0.1, , , ). Horizontal axis in each panel represents the angular position of the limb segment (θ) and is synonymous with both joint angle and muscle length. Vertical axis represents the angular velocity of the limb segment (ω) and is synonymous with muscle lengthening velocity. Faint green lines indicate randomly selected motion trajectories. Heavy blue line in each panel indicates the nullcline of ω. Its zero-crossings denote the equilibrium positions of the limb segment (θ = 180 degrees in all panels). The curvature of the nullcline qualitatively characterizes the degree of damping in the vector field. Damping is also quantified by the eigenvalues at the equilibrium point (see text) which confirm that damping is absent under nil co-activation and increases with higher levels of co-activation.

Figure 4A reveals the existence of bistable equilibrium positions in the same joint when the muscles were configured with asymmetric CE force-length properties (θ = θ = π − 1) and subjected to 50% co-contraction. Since the nullcline of is always 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 = 0, …, 1). It reveals a supercritical pitchfork bifurcation in the equilibrium positions as co-contraction exceeds the critical value of a = a = 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.

Figure 4

(A) Phase portrait of an isolated limb segment actuated by opposing muscles with asymmetric CE force-length properties at 50% co-activation (a = a = 0.5). All other parameters are identical to those of Figure 3. One unstable equilibrium point is observed at and two stable equilibrium points are observed at and respectively. (B) Nullclines for the same limb segment after manipulating muscle co-activation (a = a = a) from 0 to 1 in steps of 0.1. (C) Bifurcation plot showing the onset of bistability in the equilibrium positions through a pitchfork bifurcation as muscle co-activation is increased. Solid lines indicate stable equilibrium points, dashed lines indicate unstable equilibrium points. The branch point (BP) marks the onset of bistability (a = 0.172). (D) Same as (A) except here a = 0., a = 0.5 and only a single stable equilibrium point (þeta = 292°) is observed. (E) Same as (B) except here a = 0.5 is held fixed while a is manipulated from 0 to 1. (F) Bifurcation plot showing the emergence of bistability through a saddle-node bifurcation as muscle activation a is increased while a = 0.5 is held fixed. The limit points (LP) mark the onset (a = 0.338) and offset (a = 0.998) of bistability.

Similarly, Figures 4D–F describe the motion of the asymmetrically actuated limb when a = 0.5 is held fixed while a 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 < 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 θ0 = π − 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) depends upon the CE resting length parameter θ0. Figure 5B plots the migration of the cusp point as θ0 is manipulated. Here we see that the bistability emerges at 100% co-contraction levels when θ0 = 3.00 rad (point A) while it emerges at 50% co-contraction levels when θ0 = 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°).

Figure 5

(A) Bistability map for the isolated limb segment with asymmetric CE force-length properties (θ0 = π − 1). Solid lines indicate the transition boundary between bistable and monostable regions. Dotted line indicates balanced co-contraction (a = a). Cusp point (CP) marks the furthest extent of the bistable region (a = a = 0.171). Points M, S, and X correspond to the conditions of medium, strong, and extreme co-contraction applied in Experiments 1 and 2. (B) Shows the migration of the cusp point along the line of balanced co-contraction when the CE resting length parameter is manipulated between θ0 = θ = θ = 0 (the limit of joint flexion) and θ0 = θ = θ = π (the midpoint of muscle). Point A (θ0 = 3.00, a = 1) marks the emergence of bistability at 100% co-contraction levels. Point B (θ0 = 2.86, a = 0.5) marks the emergence of bistability at 50% co-contraction levels. Point C (θ0 = π − 1, a = 0.171) corresponds to point CP in the previous panel. Point D (θ0 = 1.68, a = 0.154) marks the maximal extent of the bistable region.

In comparison, bistability emerges at 17.1% co-contraction levels when θ0 = π − 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 θ0 = 1.68 rad (point D) and corresponds to shifting θ0 by (π − 1.68) radians (83.7°) from the muscle midpoint – which is an extreme shift.

Numerical experiment 1

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.

Figure 6

Results of Numerical Experiment 1. (A–C) Show the trajectories of the full biomechanical limb under conditions of medium, strong, and extreme co-contraction in the presence of gravity (denoted by G). Here, medium contraction is close to 20% activation capacity. Strong co-contraction is close to 50% capacity and thus approximates the biological limit of co-contraction. Extreme co-contraction is close to 80% capacity and approximates the mathematical limit of co-contraction. Limb position is sampled at 0.1 s intervals in all three panels. (D) Tracks the vertical position of the limb tip for each co-contraction condition. Damping in the full biomechanical limb is observed to increase with higher levels of co-contraction, consistent with the findings for the isolated limb segment. See Movie S1 in Supplementary Material for an animated version of this figure.

Numerical experiment 2

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.

Figure 7

Results of Numerical Experiment 2. (A) Shows the random initial limb postures for all simulation runs (n = 200) of the full biomechanical limb. (B–D) Show the final limb postures adopted by all runs under conditions of medium, strong, and extreme co-contraction respectively. Opposing muscles have asymmetric CE force-length properties (θ = θ = π − 1) otherwise all parameters are the same as Figure 6. All limbs undergoing medium co-contraction converged to same final posture whereas those undergoing strong and extreme co-contraction converged to non-unique final postures. These observations are consistent with the findings of bistability in the isolated limb segment. See Movie S2 in Supplementary Material for an animated version of this figure.

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.

Conditions for monostability

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 for our system we need only consider the zero crossings of the nullcline of Solving Euler’s law of motion for the case where and substituting muscle equations (3–6) allows the nullcline to be expressed as

where

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

where θ + θ = 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 .

Proof. The roots of nullcline (9) are given by equation (10). When the opponent muscles (denoted “a” and “b”) have identical muscuoloskeletal properties (mom = mom, F = F, k = k, k = k, ) then equation (10) reduces to

where can be expressed in terms of by substituting as follows,

Similarly, can be expressed in terms of by substituting 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 corresponds to a Gaussian curve centered on with a peak amplitude of (a − a)/k. Notice also that corresponds to the difference of two hyperbolics which asymptotes vertically to +∞ at and to −∞ at with a zero crossing at Functions and are both continuous on the interval

We argue geometrically that functions and always intersect exactly once on the interval for any value of (a − a) 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.

Case 1: aa = ab

Here thus equation (14) reduces to

which simplifies to the unique solution

after rearranging and canceling redundant terms. Thus nullcline (9) has exactly one root when a = a.

Case 2: aa > ab

Here is a continuous, non-negative function that is monotonically increasing on the interval and monotonically decreasing on the interval whereas is a continuous, monotonically decreasing function that is non-negative on the interval and non-positive on the interval

Since and are both non-negative in the interval wherein is monotonically increasing and is monotonically decreasing from its upper bound of +∞ to its lower bound of zero then there must exist exactly one solution in where is satisfied.

On the other hand, is always positive and is always non-positive on the interval Hence there are no solutions in this interval that satisfy Consequently, exactly one solution to equation (14) exists in this case across the entire interval Thus nullcline (9) has exactly one root when a > a.

Case 3: aa < ab

Similar to case 2 except here is a continuous, non-positive function that is monotonically decreasing on the interval and monotonically increasing on the interval whereas is unchanged.

Following the same reasoning as above, no solutions to exist in the interval and exactly one solution exists in the interval Consequently, exactly one solution to equation (11) exists in this case across the entire interval Thus nullcline (9) has exactly one root when when a < a.

Discussion

Our objective was to provide a simplified anatomical forward model of limb movement that replicates the fundamental dynamical relationships between muscle co-contraction and limb movement. Limb anatomy was reduced to a planar three-link rigid body limb where each joint is actuated by an opposing pair of simplified muscle actuators having force-length-velocity properties that approximate those of natural muscle systems. Stability analysis of a pair of antagonist muscle actuators confirmed that co-contraction increases muscle damping as anticipated. The stability analysis also revealed that co-contraction induces bistable equilibrium when the force-length properties of the opponent CEs are sufficiently asymmetric. Both findings were verified in the full biomechanical limb model where overall limb damping increased with co-contraction (Numerical Experiment 1) and multiple co-existing stable equilibrium postures were evoked when co-contraction was high (Numerical Experiment 2).

The effect of co-contraction on muscle damping is best understood by inspecting equation (3) of the muscle model and observing that muscle activation a(t) modulates both the force-length and force-velocity properties of the contractile element. Muscle activation therefore modulates both isometric muscle force and muscle damping to the same extent however the opposing isometric muscle forces cancel to produce nil joint torque whereas the damping forces of both muscles unite against the common muscle movement.

On the other hand, the effect of co-contraction on joint stability is best understood in terms of the velocity nullcline for a single pair of opposing muscles. Increased co-contraction induces turning points in the velocity nullcline when the force-length curves of the opponent muscles are sufficiently asymmetric. Above a critical level of co-contraction these turnings can become large enough for the nullcline to support multiple zero crossings. In such cases, the unique stable equilibrium bifurcates to yield a pair of non-identical stable equilibria. The critical level of co-contraction at which bistability emerges is determined by the degree of asymmetry in the force-length properties of the antagonist muscles. Bistability does not emerge when the antagonist muscles have identical properties.

Implications for motor control theory

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.

Implications for neurorobotics

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.

Limitations

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.

Conclusion

Our analysis of antagonistic Voigt muscles highlights that even simple muscle systems can exhibit unexpected multistable behaviors. Insights into these complex behaviors were gleaned through the application of formal methods from dynamical system theory where knowledge of the nullclines of the dynamic variables allowed us to predict the stability of the equilibrium positions of co-contracting muscles and to map the conditions under which co-contraction induces bistable joint postures. Bistability complicates the problem of achieving unambiguous control over the forward dynamics. Its existence has practical implications for the design of non-linear biomimetic actuators and theoretical implications for accounts of biological motor control that presume antagonist muscle systems are universally monostable. Further numerical analysis using more elaborate antagonist muscle models is required to assess the impact of musculotendon series elasticity on the emergence of bistability. Empirical research is ultimately required to establish whether bistable antagonist muscle dynamics can be observed in nature.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Movies S1 and S2 for this article can be found online at http://www.frontiersin.org/neurorobotics/10.3389/fnbot.2011.00005/abstract