- Email: [email protected]
Stability analysis of the elbow with a load
ARTICLE IN PRESS
Journal of Theoretical Biology 228 (2004) 115–125
Stability analysis of the elbow with a load Peter Giesla,*, Dorothea Meisela, Jurgen . Scheurlea, Heiko Wagnerb b
a Zentrum Mathematik, TU Munchen, Boltzmannstr. 3, 85747 Garching bei Munchen, Germany . . Department Science of Motion, Institute of Sports Science, Friedrich-Schiller-University, Seidelstr. 20, Jena, Jena 07749, Germany
Received 21 May 2003; received in revised form 3 December 2003; accepted 10 December 2003
Abstract We present a model of the human elbow and study the problem of existence and stability of equilibrium states. Our main goal is to demonstrate that stable equilibrium states exist just on grounds of the mechanical properties of the muscles and the skeleton. In particular, additional control mechanisms such as reflexes are not necessary to obtain stability. We assume that the activation of flexor and extensor muscles is constant and such that the right angle is an equilibrium state. We give a complete bifurcation diagram of all equilibrium states in terms of the elbow angle, the activation of the muscles and the mass of a load. Moreover, we define a dimensionless model parameter which allows to determine whether or not there are stable equilibria at an angle of ninety degrees. It turns out that the dependency of the muscle forces on the length of the muscles is the crucial factor for the stability of such an equilibrium. r 2003 Elsevier Ltd. All rights reserved. Keywords: Biomechanics; Stability; Bifurcation
1. Introduction Stability is an important aspect of human locomotion. In this paper, we will address the following question: assuming a certain position of the elbow, how will the arm react to small perturbations. Of course, reflexes play an important role in the control of a certain movement. However, also the intrinsic components of the muscles and skeleton alone are able to provide stability (cf. Wagner (1999), Wagner and Blickhan (1999) and the literature cited there). This means, that given a small perturbation, the elbow will return to its original position without intervention of the brain or any reflexes. Thus, a position can be held in a stable way without the activation of reflexes. This property can be studied rigorously by means of mathematical analysis, as will be done in this paper. In practice, stability will be obtained by a combination of such self-stabilizing properties of the muscles and—especially for greater perturbations—by additional control mechanisms. We consider the model of a human elbow joint including the muscle properties proposed in Wagner et al. (2003) and refine it concerning the following points: we model the *Corresponding author. E-mail address: [email protected] (P. Giesl). 0022-5193/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2003.12.006
dependency of the muscle forces on the muscle length by a function following Murray et al. (2000) and we propose a new formula for the excentric branch of the Hill function in order to obtain a twice continuously differentiable Hill function at 0. In our model, we assume that the humerus is fixed at a position parallel to the ground and the ulna can only move in the plane which is orthogonal to the ground and includes the humerus. Thus, the position of the ulna is totally described by the angle between humerus and ulna. Also, we assume that the ulna is loaded by a weight of mass ml at its free end. This is a very simple model for acrobatics, where a first person lying on the ground holds a second person standing on the hands of the first person. This second person is assumed to represent the aforementioned load. Assume that the flexor and extensor muscles are innervated at a constant activation level in such a way that the elbow joint attains an equilibrium state at the right angle. If both muscles are inactivated, then the right angle is an equilibrium, but it is unstable due to the gravitational force. One of the questions which will be addressed in this paper is whether this equilibrium position can be stabilized by activating both muscles strongly enough. It turns out that the answer is yes for positive values of a certain model parameter and for a
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
116
sufficiently small mass ml of the load. This is the expected behavior in a realistic situation. For positive values of the model parameter, we find that there is a critical mass mcrit up to which we obtain stability for sufficiently strong activation. Beyond this critical mass one cannot hold the weight at an angle of 90 in a stable way on grounds of the present mathematical model. We also analyse the occurrence and stability of other equilibrium states at other angles and draw a bifurcation diagram. An important task of mathematical analysis is not only to make qualitative and quantitative predictions but also to clarify which aspects of the model are important for the stability of the system. It turns out that the dependency of the force on the length of the flexor muscles plays a crucial role regarding stability. In fact, neglecting that dependency would result in a negative value of the model parameter and would imply that the angle of 90 is always unstable. The paper is organized as follows: in Section 2, we introduce the basic equation of motion underlying our mathematical analysis. In Section 3, we state our results, which are discussed further in Section 4. In three appendices we give a detailed derivation of the mathematical description of the muscle forces and the equation of motion considered, the proofs of the results of Section 3, and data of a real person which are used to illustrate the results.
2. The equation of motion Mechanically, we think of the ulna as a solid rod of uniform density and mass mu which is free to rotate about the fixed point where ulna and humerus are joined (cf. Fig. 1). A load of mass ml is fixed at the free end of
the ulna. The forces acting on the ulna and the load are gravitational and muscle forces. We distinguish between the two gravitational forces Fu ¼ g mu acting on the ulna and Fl ¼ g ml acting on the load. Concerning the muscle forces, we only consider the three most important muscles of the elbow, namely the two flexor muscles M. biceps brachii and M. brachioradialis, and the extensor M. triceps brachii. All corresponding muscle forces FBiz ; FBrach and Fext depend on the nerval activation level, on the length of the muscle and on the contraction velocity of the muscle. Hence, we model each of these forces as the product of three functions, namely the activation level E which is assumed to be constant and in the interval ½0; 1; the function Fl depending on the length of the muscle, and the so-called Hill-function H depending on the contraction velocity of the muscle. F ¼ E Fl H:
ð1Þ
The position of the ulna within its plane of motion is described by the angle between humerus and ulna which will be denoted by b; its angular velocity will be denoted ’ The equation of motion for the ulna is then by o :¼ b: given by L’ ¼ T; where Tðb; oÞ is the torque and L :¼ Jo is the angular momentum. Here J denotes the total moment of inertia which is given as the sum of the moment of inertia of the rod representing the ulna with mass mu and length l; and the moment of inertia of the load with mass ml : J :¼ ð13 mu þ ml Þl 2 :
ð2Þ
We write the second-order equation L’ ¼ J o ’ ¼ J b. ¼ T as the following system of first-order ordinary differential equations: 8 < b’ ¼ o; 1 :o ’ ¼ Tðb; oÞ: J
l
α
ð3Þ
I Br
h
ac
uln h rac
hBiz
a
kU B
γ
h ac
r
kU
β
hB
IBiz
Bi z
humerus kH Brach kH Biz
Fig. 1. Schematic illustration of the elbow joint with the flexor muscles M. biceps brachii, M. brachioradialis and their effective moment arms hBiz ; hBrach ; respectively. The extensor muscle M. triceps brachii is considered in the model but not shown in this figure.
The total torque Tðb; oÞ is the sum of the torques of three muscle and two gravitational forces, each of these torques is given as the product of the force and the corresponding so-called effective moment arm. The activation levels of both flexor muscles are assumed to be equal and are denoted by Eflex : Following Murray et al. (2000) we set Flext ðbÞ 1; i.e. for the extensor muscle we neglect the dependency of the force on the length (cf. Appendix A.2). Using Eq. (1) for the three muscle forces we obtain the following expression for the
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
total torque: Tðb; oÞ ¼ Eext Hext ðhext ðbÞ oÞhext ðbÞ þ Eflex FlBiz ðbÞHBiz ðhBiz ðbÞ oÞhBiz ðbÞ þ Eflex FlBrach ðbÞHBrach ðhBrach ðbÞ oÞhBrach ðbÞ m u þ ml : ð4Þ
cos b g l 2 Here h ¼ hðbÞ denotes the so-called effective moment arm and we think of any factor Fl as being a function of b and of any factor H as being a function of the contraction velocity of the muscle h o: For explicit expressions of these functions, cf. Appendix A. For abbreviation, we finally set 1 Fðb; oÞ :¼ Tðb; oÞ: ð5Þ J Thus, the basic equation of motion (3) turns out to be of the form ( b’ ¼ o; ð6Þ o ’ ¼ Fðb; oÞ: From now on we force the right angle b ¼ p=2; o ¼ 0 to be an equilibrium point of system (6) by setting Eflex ¼ Eflex ðEext Þ for Eext A½0; Emax ; explicitly given by Eflex ðEext Þ ¼
hext ðp=2ÞHext ð0Þ Eext : hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
ð7Þ Lemma 1. (i) Relation (7) is sufficient and necessary for b ¼ p=2; o ¼ 0 to be an equilibrium point of system (6). (ii) Eext ; Eflex A½0; 1 is equivalent to Eext A½0; Emax ; ( Emax :¼
where
1
if Eflex ð1Þp1;
1 Eflex ð1Þ
otherwise:
ð8Þ
ð9Þ
Proof. In order to obtain an equilibrium, the right-hand side of the ODE system (6) has to be zero. Thus, o ¼ 0 and Fðb; 0Þ ¼ 0: This implies (i).
k :¼
117
positive. Thus, for the maximal activation either Eflex or Eext is equal to 1. We consider two cases: If the corresponding value Eflex for Eext ¼ 1 satisfies Eflex ð1Þp1; then 1 is the maximal value
1 for Eext : Otherwise, Eflex (1) is the maximal value for Eext : & Note that for Eflex ¼ Eflex ðEext Þ with relation (7), ðp=2; 0Þ is an equilibrium state for any mass ml of the load. This relies on the fact that Fðp=2; 0Þ does not depend on ml since cosðp=2Þ ¼ 0:
3. Results In this section, we state a number of theorems. The proofs are given in Appendix B. In Theorem 1, we determine the stability of the equilibrium ðp=2; 0Þ in terms of ml and Eext : It turns out that this equilibrium is asymptotically stable provided that the mass of the load ml is sufficiently small and—under suitable assumptions—the activation Eext is sufficiently large. In Theorem 2, we give a formula for a critical mass mcrit of the load which has the following meaning: for masses ml omcrit the equilibrium state ðp=2; 0Þ is asymptotically stable for sufficiently large activation levels and for masses ml > mcrit the equilibrium state ðp=2; 0Þ is unstable for all activation levels. In additional theorems, we study the occurrence of other equilibria at angles bap=2: In Theorem 3, we assume maximal activation, i.e. Eext ¼ Emax : We determine all equilibrium points ðb; 0Þ and their stability in terms of b and ml ; and we draw a corresponding bifurcation diagram. Theorem 3 is a special case of Theorem 4 where we, more generally, let Eext A½0; Emax and determine all equilibrium points ðb; 0Þ and their stability in terms of b; ml and Eext : Finally, we explain in Theorem 5 how unstable equilibria bound the basin of attraction of an asymptotically stable equilibrium in some sense. In an example, we illustrate the results using the data of a person given in Appendix C. Before we state the first theorem, let us define the following quantity k: Here, 0 denotes the derivative.
HBiz ð0Þ½h0Biz ðp=2ÞFlBiz ðp=2Þ þ hBiz ðp=2ÞFl0Biz ðp=2Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ HBrach ð0Þ½h0Brach ðp=2ÞFlBrach ðp=2Þ þ hBrach ðp=2ÞFl0Brach ðp=2Þ h0 ðp=2Þ
ext : þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ hext ðp=2Þ
To prove (ii) we will determine the maximal value of Eext such that both Eext and Eflex attain values in ½0; 1: Note that the factor in Eq. (7) is positive since hBiz ðp=2Þ; hBrach ðp=2Þo0 while all other quantities are strictly
ð10Þ
We now introduce f to denote the rate of change of F with respect to the angle b: We consider f as a function of the extensor activation Eext and the mass of the load ml : More precisely, we
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
118
define p 1 Hext ð0Þ
kE h fðEext ; ml Þ ¼ 1 ext ext 2 ð3 mu þ ml Þl 2 m i u þ ml : ð11Þ þgl 2 Note that fðEext ; ml Þ ¼ Fb ðp=2; 0Þ; where Fb denotes the derivative of F with respect to b and where we fix the activations Eext and Eflex such that the right angle is an equilibrium, i.e. we set Eflex ¼ Eflex ðEext Þ cf. Eq. (7). h
Theorem 1. With the notations in Eqs. (10) and (11) we have: (i) If fðEext ; ml Þ > 0; then the equilibrium ðp=2; 0Þ is a saddle and thus unstable. (ii) If fðEext ; ml Þo0; then the equilibrium ðp=2; 0Þ is asymptotically stable. (iii) ð13 mu þ ml Þl 2 fðEext ; ml Þ is strictly monotonically increasing with respect to ml :
We have limb-p=2 ml ðbÞ ¼ mcrit : Thus, adding the trivial branch of equilibria b ¼ p=2 for all ml X0 we can draw a complete bifurcation diagram (cf. the example at the end of this section and Fig. 2). Moreover, the trivial branch is asymptotically stable for ml A½0; mcrit ) and unstable for ml > mcrit : The bifurcating branch near ðmcrit ; p=2Þ is unstable for momcrit and asymptotically stable for ml > mcrit and again changes its stability at angles b0 with ðml Þb ðb0 Þ ¼ 0: There are no periodic orbits. In the following theorem, we allow Eext to attain any fixed value in the interval ½0; Emax : Theorem 4. Assume k > 0 and let Eext A½0; Emax : Then all equilibria ðb; 0Þ of system (6) with bap=2 are given by the function
ml ðb; Eext Þ :¼
p hBiz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBrach ðbÞFlBrach ðbÞHBrach ð0Þ Eext Hext ð0Þ mu hext ðbÞ hext
: cos b g l 2 hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ 2 (iv) fðEext ; ml Þ is strictly monotonically decreasing with respect to Eext ; if and only if k > 0: Theorem 2. Assume that k > 0: Let Eext Hext ð0Þhext ðp=2Þ mu k ; mcrit ðEext Þ :¼ gl 2 mcrit :¼ mcrit ðEmax Þ ¼
Emax Hext ð0Þhext ðp=2Þ mu k : gl 2
ð12Þ ð13Þ
Then we have (i) If ml omcrit ðEext Þ; then ðp=2; 0Þ is asymptotically stable. (ii) If ml > mcrit ðEext Þ; then ðp=2; 0Þ is unstable. (iii) For loads with mass ml > mcrit the equilibrium ðp=2; 0Þ is unstable for each activation Eext A½0; Emax : (iv) For loads with mass ml omcrit there is an activation E* ext A½0; Emax such that the equilibrium ðp=2; 0Þ is asymptotically stable for all Eext AðE* ext ; Emax : In particular, the equilibrium is asymptotically stable for Eext ¼ Emax : Theorem 3. Assume k > 0 and set Eext ¼ Emax : Then all equilibria ðb; 0Þ of system (6) with bap=2 are given by the function ml ðbÞ :¼ p Emax Hext ð0Þ hext ðbÞ hext cos b g l 2
hBiz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBrach ðbÞFlBrach ðbÞHBrach ð0Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
mu : 2
ð14Þ
ð15Þ
We have limb-p=2 ml ðb; Eext Þ ¼ mcrit ðEext Þ: Thus, adding the plane of trivial equilibria b ¼ p=2 for ml X0 and Eext A½0; Emax we can draw a complete bifurcation diagram. Moreover, the trivial branch is asymptotically stable for ml A½0; mcrit ðEext ÞÞ and unstable for ml > mcrit ðEext ). The bifurcating branch near ðmcrit ðEext Þ; p=2Þ is unstable for momcrit ðEext Þ and asymptotically stable for ml > mcrit ðEext Þ and again changes its stability at angles b0 with ðml Þb ðb0 ; Eext Þ ¼ 0: There are no periodic orbits. For an example of the situation of the following Theorem 5 cf. the example at the end of this section and Fig. 3. Theorem 5. Assume that ðbs ; 0Þ is an asymptotically stable equilibrium of system (6). Assume additionally that bu oð>Þ bs is an equilibrium such that Fðb; 0Þo0 ð> 0Þ for all bobu ðb > bu Þ and bA½0; p (i.e. in particular that ðbu ; 0Þ is the only equilibrium for boð>Þ bs and bA½0; p; and bu is unstable). Then the set fðb; oÞ j boð>Þ bu ; opðXÞ 0g does not belong to the basin of attraction of ðbs ; 0Þ: Example. As an example we use the data of the person given in Appendix C. The value of k (cf. Eq. (10)) is in this case k ¼ 1:14 > 0: We have Emax ¼ 1; because of Eflex ð1Þ ¼ 0:50p1: We calculate mcrit ¼ 15:74 kg using Eq. (13) of Theorem 2. We assume maximal activation
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
2 +
( mcrit , π2− )
1.5
angle β [rad]
(m*, β *) 1
119
b E1:298 is asymptotically stable. For masses ml > m we have no stable equilibrium for activations Eext ¼ Emax ; Eflex ¼ Eflex ðEext Þ: Fig. 3 shows the behavior of solutions as an illustration of Theorem 5 with ml ¼ 26 kg and Eext ¼ Emax : For these parameter values ðbs ; 0Þ with bs ¼ 1:403 is an asymptotically stable equilibrium. It is surrounded by the two unstable equilibria ðbu1 ; 0Þ and ðbu2 ; 0Þ with bu1 ¼ 1:152 and bu2 ¼ p=2: They mark the border of the basin of attraction of ðbs ; 0Þ in the sense of Theorem 5.
+ 0.5
4. Discussion
0
5
10
15
20
25
30
mass of the load ml [kg]
Fig. 2. All equilibria and their stability for Eext ¼ Emax depending on the mass ml of the load. Dashed lines denote unstable, straight lines stable equilibria. The signs denote the sign of Fðb; 0Þ in the corresponding regions of the ðb; ml Þ-plane, cf. the proof of Theorem 3 in Appendix B.2. The equilibrium ðp=2; 0Þ is asymptotically stable for ml omcrit : For mcrit oml om the bifurcating equilibrium is asymptotically stable.
3 Φ (β,0) 2 1 1
1.1
1.2
1.3
1.4
0
1.5
1.6
1.7 β
-1 Ð2 ω
. β=ω > 0
Ð3 β Ð4 . β =ω < 0
Fig. 3. Upper diagram: Fðb; 0Þ with ml ¼ 26 kg and Eext ¼ Emax ; lower diagram: the three equilibria ð1:152; 0Þ (unstable), ð1:403; 0Þ (asymptotically stable) and ðp=2; 0Þ (unstable) and the behavior of solutions of system (6) in the b–o plane.
Eext ¼ Emax and draw the bifurcation diagram in Fig. 2 using Theorem 3. The equilibrium ðp=2; 0Þ is asymptotically stable for masses smaller than mcrit : For masses between mcrit and the mass m E28:33 kg an equilibrium point at a smaller angle than p=2 bounded below by
The sign of the quantity k is crucial for most theorems. If k is negative, then for all positive loads ml the equilibrium ðp=2; 0Þ is unstable for any activation Eext by Theorem 1, since then (cf. Eq. (11)) h p 1 fðEext ; ml Þ ¼ 1
kE h Hext ð0Þ ext ext 2 ð3 mu þ ml Þl 2 m i u þ gl þ ml > 0: 2 Note that if k > 0; increasing the activation has a stabilizing effect, while for ko0 it has a destabilizing effect. We consider the influence of the different terms on the sign of k (cf. the definition of k in Eq. (10)). Since the slope of hext ; hBiz and hBrach is positive at p=2; it is necessary to have force–length relations FlBiz and FlBrach with a large positive slope at p=2 in order to obtain k > 0 (recall that hBiz ðp=2Þ; hBrach ðp=2Þo0Þ: Hence, for a realistic model it is essential to determine the force– length relation properly. In particular, a constant force– 0 length relation for the flexors would imply FBiz ðp=2Þ ¼ 0 FBrach ðp=2Þ ¼ 0 and thus ko0: On grounds of such a model the equilibrium at ðp=2; 0Þ would be unstable for any load. Moreover, regarding Eq. (13), the value of k is responsible for the critical mass mcrit you still can hold in a stable way. Hence, in order to increase the critical mass, besides augmentation of the isometric forces through training it could be important to improve the muscle–length dependency. From the bifurcation diagram we see that in our example loads with masses ml > mcrit can be held in a stable way with constant muscle activations, but at smaller angles than p=2: If the mass is augmented continuously, e.g. by filling a glass in the hand of the person with water, the equilibrium point at p=2 becomes unstable when the mass of the load exceeds mcrit : Then the elbow will—without intervention of the brain or reflexes—leave the unstable equilibrium position at p=2 (dashed line) and move to the stable position at a smaller angle (straight line).
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
120
By the bifurcation diagram, loads with mass ml > mcrit can be held in a stable way up to a certain mass m > mcrit and at a certain angle bap=2 (cf. the bifurcating branch of equilibria in Fig. 2). This shows that the right angle is not the optimal angle to hold large loads in a stable way. The question arises whether one can also hold masses beyond the mass m at an angle of 90 : In fact, as experiments have shown, it is possible to hold much higher masses, but not with constant activation of the muscles. This is in agreement with our analysis which showed that other stabilizing mechanisms are necessary to hold these masses in a stable way. In fact, one holds these heavy loads by switching on and off the muscle activations which can be shown by recording the electrical nerve stimulation of the muscles using electrodes (EMG). This situation of time-dependent activations is not included in our model. Stability is thus obtained through a combination of the self-stabilizing properties of the musculoskeletal system studied in this paper on the one hand, and the reflexes and the control through the brain on the other hand.
Appendix A. Derivation of the torque In this appendix, we give a detailed derivation of the torques corresponding to the gravitational and the muscle forces; the total torque Tðb; oÞ will be the sum of these torques. In general the torque is defined by the product of the length of the vector from the point of joint to the point of action of the applied force and the component of the force orthogonal on the aforementioned vector with appropriate sign. In our case, however, it will be simpler to determine the torque by the equivalent formula
l cosðp bÞ ¼ l cos b (cf. Fig. 1). The effective moment arm of the gravitational force acting on the ulna is hu ¼ 12l cos b due to the integration over the length of the ulna. Hence, m u Tðb; oÞ ¼ cos b g l þ ml 2 þ hBiz FBiz þ hBrach FBrach þ hext Fext : ðA:1Þ In the next sections, we will consider the torques corresponding to the muscle forces. In Section A.1, we determine the effective moment arms hBiz ; hBrach and hext and in Section A.2, we derive the corresponding forces FBiz ; FBrach and Fext : In Section A.3, finally, we will summarize the previous sections giving a formula of the total torque. A.1. The effective moment arms We calculate the effective moment arms h of each muscle. As mentioned before, we restrict ourselves to the most important muscles, namely the two flexor muscles M. biceps brachii and M. brachioradialis, and the extensor M. triceps brachii. We start with the two flexor muscles. Flexors: In Fig. 1 we marked the relevant quantities. The distances of the insertions of the muscles M. biceps brachii and M. brachioradialis to the elbow joint are denoted by kUBiz ; kHBiz ; kUBrach and kHBrach ; respectively. We calculate the effective moment arms hBrach and hBiz : Their sign will be negative. Using the sine-theorem, we get kHBrach sin a ¼ sin b ; ðA:2Þ lBrach sin g ¼ sin b
kHBiz : lBiz
ðA:3Þ
We conclude (note that the sign of hBrach and hBiz is negative)
T ¼ hF ; where F denotes the total strength of the force and h the so-called effective moment arm. It is given in the following way: draw a line which includes the point of the action of the force along the vector of the action of the force. Then there is a unique vector from the point of joint to a point of this line which is orthogonal on this line. h is 7 the length of this vector. The sign is negative, if this vector points in direction of decreasing b and positive otherwise. The total torque is the sum of the torques corresponding to the different forces Tðb; oÞ ¼ hu Fu þ hl Fl þ hBiz FBiz þ hBrach FBrach þ hext Fext : The gravitational forces are the force Fu ¼ g mu acting on each point of the ulna, and the force Fl ¼ g ml acting on the load. The effective moment arm of the gravitational force acting on the load is hl ¼
kHBrach kUBrach ; lBrach ðA:4Þ kHBiz kUBiz ðA:3Þ ¼ sin g kUBiz ¼ sin b : lBiz ðA:2Þ
hBrach ¼ sin a kUBrach ¼ sin b hBiz
The lengths of the muscles are given by the cosinetheorem as follows: lBrach ðbÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ¼ kHBrach þ kUBrach
2kHBrach kUBrach cos b;
ðA:5Þ
lBiz ðbÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 þ kUBiz
2kHBiz kUBiz cos b: ¼ kHBiz
ðA:6Þ
Extensor: In order to calculate the effective moment arm of the extensor we cannot give a geometric argumentation as in the case of the flexors. Instead, we
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
follow Gerbeaux et al. (1996) who approximated the effective moment arm (which they call geometry function) by a polynomial of third degree, namely hext ðbÞ ¼ l½0:0993 þ 0:0371ðp bÞ 0:0524ðp bÞ2 þ 0:0147ðp bÞ3 :
121
160 ð88p=9Þ: Hence, p p zBiz lBiz ¼ 0:45 ) c1 lBiz þ c2 ¼ 0:45; 3 3
ðA:7Þ zBiz
8p 8p lBiz ¼ 1 ) c1 lBiz þ c2 ¼ 1; 9 9
A.2. The muscle forces from which we can conclude Now we consider the muscle forces FBiz ; FBrach and Fext : Each of them depends on the activation of the muscle, the length of the muscle and the velocity of contraction of the muscle. We write the forces as the product of three functions, namely the activations Eflex ; Eext ; respectively, the force–length functions FlBiz ðbÞ and FlBrach ðbÞ; and the velocity–length functions HBiz ðb; oÞ; HBrach ðb; oÞ; Hext ðb; oÞ; respectively, which are also called the Hill functions. Note that we set the force– length function of the extensor Flext ðbÞ :¼ 1: FBiz ðb; oÞ ¼ Eflex FlBiz ðbÞ HBiz ðhBiz ðbÞ oÞ;
zBiz ðxÞ ¼
0:55 lBiz ð8p=9Þ lBiz ðp=3Þ
8p x lBiz þ 1: 9
A similar formula holds for zBrach : Hence, FlBiz ðbÞ ¼
ðA:8Þ
1:1 arctan½80ðzBiz ðlBiz ðbÞÞ 1:14Þ9 p þ 7ðzBiz ðlBiz ðbÞÞ 0:57Þ þ 0:55;
FBrach ðb; oÞ ¼ Eflex FlBrach ðbÞ HBrach ðhBrach ðbÞ oÞ; ðA:9Þ FlBrach ðbÞ ¼ Fext ðb; oÞ ¼ Eext Hext ðhext ðbÞ oÞ:
ðA:10Þ
Activation: We assume that both flexor muscles are activated at the same activation level EBiz ¼ EBrach ¼: Eflex : The activation of the muscle is given as the amount of the maximal activation and will be assumed to be constant in our model. Thus, the activations have to fulfill Eflex ; Eext A½0; 1: Force–length relation: The length-dependency of the (maximal) force which can be generated by the muscle has been studied by Murray et al. (2000) analysing dead bodies. Depending on the numbers of overlapping filaments the force is maximal at a medium length and decreases both if the muscle is longer and if it is shorter than that (cf. Gordon et al., 1966). The curve is the same for all muscles but the region in which the muscle works is different for different muscles. For the extensor, the force–length function is approximately 1 and thus we set Flext ðbÞ ¼ 1: For both flexors (M. biceps and M. brachioradialis), we model the function given in Murray et al. (2000) (at the left-hand region) by Fl ðzÞ :¼
1:1 arctan½80ðz 1:14Þ9 þ 7ðz 0:57Þ p þ 0:55;
ðA:11Þ
where z denotes the normalized muscle length. We assume zðxÞ which maps the total muscle length x to the normalized muscle length to be a linear relation of the form zBiz ðxÞ ¼ c1 x þ c2 : Following Murray et al. (2000) we postulate that the flexors attain 45% of their optimal length at an angle of 60 ð8p=3Þ and 100% at
ðA:12Þ
ðA:13Þ
1:1 arctan½80ðzBrach ðlBrach ðbÞÞ 1:14Þ9 p þ 7ðzBrach ðlBrach ðbÞÞ 0:57Þ þ 0:55: ðA:14Þ
Force–velocity relation: the Hill function: A common model for the dependency of the muscle force on the velocity in the literature is the so-called Hill force– velocity relation. Hill (1938) only studied the concentric part, i.e. movements, where the muscle is shortened. We, however, will need a function which is also defined for the excentric part and is C 2 ; and therefore we suggest the following Hill function: Hðvm Þ 8 c > > < vm þ b a ¼ B C > > þ :A þ vm D ðvm DÞ2
for vm X0; for vm o0;
ðA:15Þ
where a; b; c and A; B; C; D are positive parameters. The Hill functions of both the flexors and the extensor are of this form with different parameters. The contraction velocity vm is given by the derivative of the negative muscle length lm with respect to the time t: vm ¼
dlm : dt
The value Hð0Þ is called the isometric force Fiso :¼ c=b a: For vm - N the muscle can generate a stronger force, i.e. H tends to Fiso cexz ; where cexz > 1: In general, one assumes cexz ¼ 1:5: Together with the requirement that HAC 2 ðR; RÞ; the four parameters
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
122
A; B; C; D are determined, if a; b; c and cexz are given.
Dempster, 1955):
A ¼ cexz Fiso > 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b3 D ¼ b þ b2 þ Fiso ðcexz 1Þ > 0; c c B ¼ ð3b þ 2DÞ Fiso ðcexz 1Þ 2D 2 > 0; b c > 0: C ¼ D Fiso ðcexz 1ÞðD bÞ þ 2D b
mu ¼ 0:022m:
The derivative of the Hill function is strictly negative for all vm AR: For vm X0; this follows from c=ðvm þ bÞ2 o0: For vm o0 we have to show that ðBðvm DÞ þ 2CÞ=ðvm DÞ3 o0: It is sufficient to show 2CoBD: This follows by a straightforward calculation using the above definitions of B; C and D: In order to derive a relation between the muscle velocity vm and the velocity o; we use the following equality: the power of the muscle Fext vmext equals the power Text o: Thus, we have Fext vmext ¼ Fext hext o; vmext ¼ hext o:
ðA:16Þ
For the flexor muscles we have in an analogous way
Sust (1996) developed a method to determine the parameters of the Hill functions for individual persons through maximal voluntary contraction. We used this method to determine the parameters a; b; c of the Hill functions, and we additionally assumed that they are equal for both flexor muscles M. biceps brachii and M. brachioradialis. Appendix B. Proofs B.1. Stability of the equilibrium ðp=2; 0Þ Proof of Theorem 1. We denote the derivatives of F with respect to b and o by Fb and Fo ; respectively. The Jacobian of the right-hand side of the system of ordinary differential equations (6) is ! 0 1 ðB:1Þ Fb ðb; oÞ Fo ðb; oÞ
vmBiz ¼ hBiz o;
ðA:17Þ
and its eigenvalues at the equilibrium ðb; 0Þ are qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l1;2 ¼ 12 ðFo ðb; 0Þ7 Fo ðb; 0Þ2 þ 4Fb ðb; 0ÞÞ:
vmBrach ¼ hBrach o:
ðA:18Þ
Note that
A.3. The total torque We have described the muscle forces by Eqs. (A.8)– (A.10). Thus, the total torque is given by (cf. Eq. (A.1)) Tðb; oÞ ¼ hext Fext þ hBiz FBiz þ hBrach FBrach m u þ ml
cos b g l 2 ¼ Eext hext ðbÞHext ðhext ðbÞ oÞ þ Eflex hBiz ðbÞFlBiz ðbÞHBiz ðhBiz ðbÞ oÞ þ Eflex hBrach ðbÞFlBrach ðbÞHBrach ðhBrach ðbÞ oÞ m u þ ml ;
cos b g l 2 which corresponds to Eq. (4). As an example the data of a specific person are given in the Appendix C. The lengths of the ulna and the humerus are taken from those data, and we assume l ¼ lulna : The insertions of the muscles are determined following the studies of Wank (1996) kHBrach :¼ 0:104 lhumerus ; kHBiz :¼ lhumerus ; kUBrach :¼ lulna ; kUBiz :¼ 0:186 lulna : The mass of the arm is determined by the following formula through the total mass m of the person (cf.
ðB:2Þ
Fo ðb; oÞ 1 0 ðhext ðbÞ oÞh2ext ðbÞ ¼ ½Eext Hext J 0 þ Eflex FlBiz ðbÞHBiz ðhBiz ðbÞ oÞh2Biz ðbÞ 0 ðhBrach ðbÞ oÞh2Brach ðbÞ þ Eflex FlBrach ðbÞHBrach
ðB:3Þ
is strictly negative for all values of b and o; because all components are strictly positive except for the deriva0 0 0 tives of the Hill functions Hext ; HBiz and HBrach which are strictly negative (cf. Appendix A). Thus, there is at least one eigenvalue with negative real part. The second eigenvalue has negative real part if and only if Fb ðb; 0Þo0: Depending on the sign of Fo ðb; 0Þ2 þ 4Fb ðb; 0Þ we have either real or complex eigenvalues which means the equilibrium is either a node or a sink. The stability of ðp=2; 0Þ is determined by the sign of Fb ðp=2; 0Þ which we denote by fðEext ; ml Þ—note that we have used Eq. (7) to derive Eq. (11). This proves (i) and (ii). (iii) and (iv) follow directly from Eq. (11). & Proof of Theorem 2. Defining mcrit ðEext Þ implicitly by fðEext ; mcrit Þ ¼ 0; we have fðEext ; ml Þo0 ð> 0Þ for ml oð>Þ mcrit ðEext Þ: This relies on the specific dependency of f on ml and Eext which was shown in Theorem 1. This proves (i) and (ii). To prove (iii) we have ml > mcrit ¼ mcrit ðEmax Þ Xmcrit ðEext Þ for all Eext A½0; Emax and thus by (ii) the equilibrium ðp=2; 0Þ is unstable. For (iv) note that for
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
ml omcrit there is one and only one value E* ext A½0; Emax such that ml ¼ mcrit ðE* ext Þ since mcrit ð0Þo0; mcrit ðEmax Þ ¼ mcrit > ml and mcrit ðEext Þ is strictly monotonically increasing with respect to Eext : Thus, (iv) follows by (i). &
B.2. The bifurcating branch of equilibria and the basin of attraction Proof (of Theorems 4 and 3). We prove Theorem 4. Theorem 3 follows from Theorem 4 setting Eext :¼ Emax : The formula for the mass is derived from the condition Fðb; 0Þ ¼ 0 (cf. also Eqs. (5) and (7)). The # limit is calculated using l’Hospital: lim ml ðb; Eext Þ Eext Hext ð0Þ ½h0ext ðbÞ ¼ lim b-p=2 sin b g l
b-p=2
hext ðp=2Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
ðh0Biz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBiz ðbÞFl0Biz ðbÞHBiz ð0Þ þ h0Brach ðbÞFlBrach ðbÞHBrach ð0Þ þhBrach ðbÞFl0Brach ðbÞHBrach ð0ÞÞ
mu
2
¼ mcrit ðEext Þ ðcf: Eqs: ð10Þ and ð12ÞÞ: We calculate the trace of the Jacobian, which is Fo ðb; oÞo0 (cf. Eq. (B.1) and Appendix A). By the Bendixson criterion we can conclude that there is no periodic solution. Fix Eext A½0; Emax : To illustrate the argumentation also consider Fig. 2 which shows the bifurcation diagram for a certain person the data of whom are reported in Appendix C. The stability of the trivial branch at b ¼ p=2 is determined by the sign of Fb ðb; 0Þ; which is negative for ml omcrit ðEext Þ and which is positive for ml > mcrit ðEext Þ; cf. Theorem 2. Hence, we can determine the sign of Fðb; 0Þ for fixed ml near the trivial branch. Since we know that we marked all zeros of Fðb; 0Þ; we know the sign of this function for all b; because by the intermediate value theorem it cannot change its sign unless it reaches zero. Now we know the sign of Fðb; 0Þ everywhere and thus conclude the stability of the nontrivial branch. Near ðmcrit ðEext Þ; p=2Þ all different cases result in asymptotic stability for the bifurcating equilibria if ml > mcrit ðEext Þ; and instability if ml omcrit ðEext Þ: Assume that ðb0 ; 0Þ is an equilibrium and that the sign of Fðb; 0Þ changes from to þ at b0 for increasing b: Then ðb0 ; 0Þ is unstable. If Fðb; 0Þ changes its sign from þ to
at b0 for increasing b; then ðb0 ; 0Þ is asymptotically stable. In case of hyperbolic equilibria this is ensured through linearization, but even in the non-hyperbolic
123
case this is true as we show in Proposition 2 in Section B.3. Thus, we have determined the stability of all equilibria except for the critical ones. This proves Theorem 4. & Proof of Theorem 5. For illustration cf. Fig. 3. We only consider the case of an equilibrium bu obs ; a similar argumentation holds for the other case. We know that for op0 we have b’ ¼ op0 and solutions cannot cross the b-axis, because then o ’ ¼ Fðb; 0ÞX0 which contradicts the assumption. Thus, the solution with initial value ðb; oÞ; 0pbobu ; op0 cannot tend to ðbs ; 0Þ: &
B.3. A proposition using a center manifold We prove a general proposition which links the stability properties of equilibria of the first-order system b’ ¼ Fðb; 0Þ to the stability properties of corresponding ’ In equilibria of the second-order system b. ¼ Fðb; bÞ: case of a hyperbolic equilibrium one has to check the signs of the eigenvalues, but the proposition is also valid in the non-hyperbolic case, where we use a center manifold for the proof. Note that we need FAC 2 ðR2 ; RÞ in order to use the center manifold theorem. Proposition 2. Let FAC 2 ðR2 ; RÞ: Consider the two ordinary differential equations ( b’ ¼ o; ðB:4Þ o ’ ¼ Fðb; oÞ; b’ ¼ Fðb; 0Þ;
ðB:5Þ
b0 is an equilibrium of Eq. ðB:5Þ if and only if ðb0 ; 0Þ is an equilibrium of Eq. ðB:4Þ: If b0 is an unstable/asymptotically stable equilibrium of Eq. ðB:5Þ such that Fo ðb0 ; 0Þo0 holds, then the corresponding equilibrium ðb0 ; 0Þ of Eq. ðB:4Þ is unstable/asymptotically stable. Proof. ðb0 ; oÞ is an equilibrium of Eq. (B.4), if and only if o ¼ 0 and Fðb0 ; 0Þ ¼ 0: Now assume that b0 is an equilibrium and that Fo ðb0 ; 0Þo0: The Jacobian of the right-hand side of Eq. (B.4) is ! 0 1 : Fb ðb; oÞ Fo ðb; oÞ The eigenvalues at ðb ; 0Þ are l1;2 ¼12ðFo ðb0 ; 0Þ7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞÞ: Let us first assume that the equilibrium b0 of Eq. (B.5) is hyperbolic. If b0 is unstable with respect to Eq. (B.5), then Fb ðb0 ; 0Þ > 0: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi We can conclude Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0Þ > jFo ðb0 ; 0Þj and hence one eigenvalue is negative and the other positive. Thus, ðb0 ; 0Þ is a saddle and unstable. If b0 is asymptotically stable, then Fb ðb0 ; 0Þo0: Moreover,
ARTICLE IN PRESS 124
P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
either Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0Þo0 in which case we have two complex conjugate eigenvalues with real part 1 Fo ðb0 ; 0Þo0; or Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞX0 and then 2q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞojFo ðb0 ; 0Þj and both eigen values are real and negative. In both cases ðb0 ; 0Þ is asymptotically stable. Now assume that the equilibrium b0 of Eq. (B.5) is not hyperbolic. Then we have Fb ðb0 ; 0Þ ¼ 0: If the equilibrium is asymptotically stable, then Fðb0 þ x; 0Þ is negative for small x > 0 and positive for small xo0: If the equilibrium is unstable, then Fðb0 þ x; 0Þ is positive for small x > 0 or negative for small xo0: In the following we assume in the unstable case that Fðb0 þ x; 0Þ is positive for small x > 0: The other case can be dealt with similarly. We show more generally that the behavior of the one-dimensional system (B.5) is the same as the behavior on the center manifold of the twodimensional system (B.4). Since Fb ðb0 ; 0Þ ¼ 0 the eigenvalues of Eq. (B.4) satisfy l1 ¼ 0 and l2 o0: Thus, there is a one-dimensional stable and a one-dimensional center manifold. The stability of ðb0 ; 0Þ with respect to Eq. (B.4) is determined by the behavior on the center manifold. To study the dynamics we first determine the equation of the center manifold, then we study the behavior on the center manifold and finally we conclude the stability of the two-dimensional equilibrium. To determine the equation of the center manifold, we first apply a coordinate transformation in order to transform the system of differential equations into a suitable form. In fact, we choose the eigenvectors to be the basis of the new coordinate system. We set c :¼ Fo ðb0 ; 0Þo0: The Jacobian in this case ðFb ðb0 ; 0Þ ¼ 0Þ reads ! 0 1 : 0 c The eigenvalues are 0 with corresponding eigenvector ð1; 0ÞT ; and co0 with corresponding eigenvector ð1; cÞT : Hence, we define the following new coordinates: ! ! ! b b0 x 1 1c ¼ : y 0 1c o The system of differential equations in the new coordinates reads 8 1 > < x’ ¼ cy Fðb0 þ x þ y; cyÞ; c ðB:6Þ > : y’ ¼ 1 Fðb þ x þ y; cyÞ: 0 c ð0; 0Þ is the corresponding equilibrium of Eq. (B.6). By the theorem of the existence of a center manifold (cf. Carr, 1981, Theorem 1), we can describe a center manifold by a C 2 -function hðxÞ; such that locally near ð0; 0Þ a point ðx; yÞ lies on the center manifold if and only
if y ¼ hðxÞ: Moreover, hð0Þ ¼ h0 ð0Þ ¼ 0: To determine h we use the invariance of the center manifold, which implies y’ ¼ h0 ðxÞx’ and hence 1 Fðb0 þ x þ hðxÞ; chðxÞÞ c
1 ðB:7Þ ¼ h0 ðxÞ chðxÞ Fðb0 þ x þ hðxÞ; chðxÞÞ : c By Taylor’s formula there is a yA½0; 1 such that Fðb0 þ x þ hðxÞ; chðxÞÞ ¼ Fðb0 þ x; 0Þ þ Fb ðb0 þ x þ yhðxÞ; ychðxÞÞ hðxÞ þ Fo ðb0 þ x þ yhðxÞ; ychðxÞÞ chðxÞ:
ðB:8Þ
Plugging Eq. (B.8) into Eq. (B.7) we get ½Fðb0 þ x; 0Þ þ Fb ðb0 þ x þ yhðxÞ; ychðxÞÞ hðxÞ 1 þ Fo ðb0 þ x þ yhðxÞ; ychðxÞÞ chðxÞ ð1 þ h0 ðxÞÞ c ¼ ch0 ðxÞhðxÞ: If hðxÞ ¼ 0; then the above formula shows that Fðb0 þ x; 0Þ ¼ 0: But if b0 is asymptotically stable, then we have hðxÞa0 for all x in a neighborhood of zero. If on the other hand b0 is unstable, then we have hðxÞa0 for all xAð0; eÞ with an e > 0; by the assumption Fðb0 þ x; 0Þ > 0 for those x: Thus, hðxÞa0 holds in a neighborhood of zero, if b0 is asymptotically stable, and for all xAð0; eÞ; if b0 is unstable. Hence, we are allowed to divide both sides by hðxÞ and then let x-0; if b0 is asymptotically stable, and xr0; if b0 is unstable. Since hð0Þ ¼ h0 ð0Þ ¼ 0; we get with Fb ðb0 ; 0Þ ¼ 0 Fðb0 þ x; 0Þ ¼ c2 o0: ðB:9Þ lim x-0 hðxÞ Hence, Fðb0 þ x; 0Þ and hðxÞ have opposite signs for small jxj: The dynamics on the center manifold are governed locally near 0 by (cf. Carr, 1981, Theorem 1) 1 u’ ¼ chðuÞ Fðb0 þ u þ hðuÞ; chðuÞÞ: ðB:10Þ c We can expand the right-hand side using Eq. (B.8) again. There is a yA½0; 1 such that the right-hand side of Eq. (B.10) reads 1 chðuÞ ½Fðb0 þ u; 0Þ þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ hðuÞ c þ Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ chðuÞ 1 Fðb0 þ u; 0Þ ¼ c
þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ c hðuÞ
þ c Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ hðuÞ: Since by Eq. (B.9) Fðb0 þ u; 0Þ þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ lim u-0 hðuÞ þ c Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ ¼ c2 þ 0 þ c2 ¼ 0;
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
the sign of the right-hand side of Eq. (B.10) is determined by chðuÞ for small juj: Since co0; it has the opposite sign as compared to hðuÞ near u ¼ 0 and thus the same sign as Fðb0 þ u; 0Þ: Since the right-hand side of Eq. (B.10) has the same sign as Fðb0 þ u; 0Þ for small u as we just have shown, the stability properties of b0 with respect to Eq. (B.5) and u ¼ 0 with respect to Eq. (B.10) are the same. By Carr (1981, Theorem 2), the asymptotic stability/ instability on the center manifold (B.10) implies the asymptotic stability/instability of the two-dimensional system (B.4) and thus the theorem is proved. &
Appendix C. Data
Constant Description
Value
Unit
lulna lhumerus m aext
0.26 0.31 55 403.7752
m m kg N
bext cext cexzext aflex bflex
Length of the ulna Length of the humerus Mass of the person Hill function parameter (extensor) Hill function parameter (extensor) Hill function parameter (extensor) Excentric force increase factor (extensor) Hill function parameter (both flexors) Hill function parameter (both flexors)
0.2440 m/s 502.3032 W 1.5 21.4177 N 0.0233 m/s
cflex cexzflex
Hill function parameter (both flexors) Excentric force increase factor (both flexors)
125
30.7825 W 1.5
References Carr, J., 1981. Applications of Centre Manifold Theory. Springer, Berlin. Dempster, W.T., 1955. Space requirements of the seated operator (WADC-TR-55-159). Wright-Patterson Air Force Base, OH. Gerbeaux, M., Turpin, E., Lensel-Corbeil, G., 1996. MusculoArticular modelling of the triceps brachii. J. Biomech. 29, 171–180. Gordon, A.M., Huxley, A.F., Julian, F.J., 1966. The variation in isometric tension with sarcomere length in vertebrate muscle fibres. J. Physiol. 184, 170–192. Hill, A.V., 1938. The heat of shortening and the dynamic constants of muscle. Proc. Roy. Soc. (B) Lond. 126, 136–195. Murray, W., Buchanan, T., Scott, D., 2000. The isometric functional capacity of muscles that cross the elbow. J. Biomech. 33, 943–952. Sust, M., 1996. Modular aufgebaute deterministische Modelle zur Beschreibung menschlicher Bewegungen. In: Ballreich, R., Baumann, W. (Eds.), Grundlagen der Biomechanik des Sports. Ferdinand-Enke-Verlag, Stuttgart. Wagner, H., 1999. Zur Stabilit.at der menschlichen Bewegung. Dissertation, J.W. Goethe-Universit.at, Frankfurt a.M. Wagner, H., Blickhan, R., 1999. Stabilizing function of skeletal muscles: an analytical investigation. J. Theor. Biol. 199, 163–179. Wagner, H., Machalett, M., Blickhan, R., 2003. Bewegungslernen und . Attraktorstabilit.at—vorl.aufige Ergebnisse und mogliche Anwendungen. In: Riele, H. (Ed.), Biomechanik als Anwendungsforschung—Transfer zwischen Theorie und Praxis, Vol. 132. Czwalina Verlag, Hamburg. Wank, V., 1996. Modellierung und Simulation von Muskelkontrak. tionen. Sport und Buch StrauX GmbH, Koln.
Journal of Theoretical Biology 228 (2004) 115–125
Stability analysis of the elbow with a load Peter Giesla,*, Dorothea Meisela, Jurgen . Scheurlea, Heiko Wagnerb b
a Zentrum Mathematik, TU Munchen, Boltzmannstr. 3, 85747 Garching bei Munchen, Germany . . Department Science of Motion, Institute of Sports Science, Friedrich-Schiller-University, Seidelstr. 20, Jena, Jena 07749, Germany
Received 21 May 2003; received in revised form 3 December 2003; accepted 10 December 2003
Abstract We present a model of the human elbow and study the problem of existence and stability of equilibrium states. Our main goal is to demonstrate that stable equilibrium states exist just on grounds of the mechanical properties of the muscles and the skeleton. In particular, additional control mechanisms such as reflexes are not necessary to obtain stability. We assume that the activation of flexor and extensor muscles is constant and such that the right angle is an equilibrium state. We give a complete bifurcation diagram of all equilibrium states in terms of the elbow angle, the activation of the muscles and the mass of a load. Moreover, we define a dimensionless model parameter which allows to determine whether or not there are stable equilibria at an angle of ninety degrees. It turns out that the dependency of the muscle forces on the length of the muscles is the crucial factor for the stability of such an equilibrium. r 2003 Elsevier Ltd. All rights reserved. Keywords: Biomechanics; Stability; Bifurcation
1. Introduction Stability is an important aspect of human locomotion. In this paper, we will address the following question: assuming a certain position of the elbow, how will the arm react to small perturbations. Of course, reflexes play an important role in the control of a certain movement. However, also the intrinsic components of the muscles and skeleton alone are able to provide stability (cf. Wagner (1999), Wagner and Blickhan (1999) and the literature cited there). This means, that given a small perturbation, the elbow will return to its original position without intervention of the brain or any reflexes. Thus, a position can be held in a stable way without the activation of reflexes. This property can be studied rigorously by means of mathematical analysis, as will be done in this paper. In practice, stability will be obtained by a combination of such self-stabilizing properties of the muscles and—especially for greater perturbations—by additional control mechanisms. We consider the model of a human elbow joint including the muscle properties proposed in Wagner et al. (2003) and refine it concerning the following points: we model the *Corresponding author. E-mail address: [email protected] (P. Giesl). 0022-5193/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2003.12.006
dependency of the muscle forces on the muscle length by a function following Murray et al. (2000) and we propose a new formula for the excentric branch of the Hill function in order to obtain a twice continuously differentiable Hill function at 0. In our model, we assume that the humerus is fixed at a position parallel to the ground and the ulna can only move in the plane which is orthogonal to the ground and includes the humerus. Thus, the position of the ulna is totally described by the angle between humerus and ulna. Also, we assume that the ulna is loaded by a weight of mass ml at its free end. This is a very simple model for acrobatics, where a first person lying on the ground holds a second person standing on the hands of the first person. This second person is assumed to represent the aforementioned load. Assume that the flexor and extensor muscles are innervated at a constant activation level in such a way that the elbow joint attains an equilibrium state at the right angle. If both muscles are inactivated, then the right angle is an equilibrium, but it is unstable due to the gravitational force. One of the questions which will be addressed in this paper is whether this equilibrium position can be stabilized by activating both muscles strongly enough. It turns out that the answer is yes for positive values of a certain model parameter and for a
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
116
sufficiently small mass ml of the load. This is the expected behavior in a realistic situation. For positive values of the model parameter, we find that there is a critical mass mcrit up to which we obtain stability for sufficiently strong activation. Beyond this critical mass one cannot hold the weight at an angle of 90 in a stable way on grounds of the present mathematical model. We also analyse the occurrence and stability of other equilibrium states at other angles and draw a bifurcation diagram. An important task of mathematical analysis is not only to make qualitative and quantitative predictions but also to clarify which aspects of the model are important for the stability of the system. It turns out that the dependency of the force on the length of the flexor muscles plays a crucial role regarding stability. In fact, neglecting that dependency would result in a negative value of the model parameter and would imply that the angle of 90 is always unstable. The paper is organized as follows: in Section 2, we introduce the basic equation of motion underlying our mathematical analysis. In Section 3, we state our results, which are discussed further in Section 4. In three appendices we give a detailed derivation of the mathematical description of the muscle forces and the equation of motion considered, the proofs of the results of Section 3, and data of a real person which are used to illustrate the results.
2. The equation of motion Mechanically, we think of the ulna as a solid rod of uniform density and mass mu which is free to rotate about the fixed point where ulna and humerus are joined (cf. Fig. 1). A load of mass ml is fixed at the free end of
the ulna. The forces acting on the ulna and the load are gravitational and muscle forces. We distinguish between the two gravitational forces Fu ¼ g mu acting on the ulna and Fl ¼ g ml acting on the load. Concerning the muscle forces, we only consider the three most important muscles of the elbow, namely the two flexor muscles M. biceps brachii and M. brachioradialis, and the extensor M. triceps brachii. All corresponding muscle forces FBiz ; FBrach and Fext depend on the nerval activation level, on the length of the muscle and on the contraction velocity of the muscle. Hence, we model each of these forces as the product of three functions, namely the activation level E which is assumed to be constant and in the interval ½0; 1; the function Fl depending on the length of the muscle, and the so-called Hill-function H depending on the contraction velocity of the muscle. F ¼ E Fl H:
ð1Þ
The position of the ulna within its plane of motion is described by the angle between humerus and ulna which will be denoted by b; its angular velocity will be denoted ’ The equation of motion for the ulna is then by o :¼ b: given by L’ ¼ T; where Tðb; oÞ is the torque and L :¼ Jo is the angular momentum. Here J denotes the total moment of inertia which is given as the sum of the moment of inertia of the rod representing the ulna with mass mu and length l; and the moment of inertia of the load with mass ml : J :¼ ð13 mu þ ml Þl 2 :
ð2Þ
We write the second-order equation L’ ¼ J o ’ ¼ J b. ¼ T as the following system of first-order ordinary differential equations: 8 < b’ ¼ o; 1 :o ’ ¼ Tðb; oÞ: J
l
α
ð3Þ
I Br
h
ac
uln h rac
hBiz
a
kU B
γ
h ac
r
kU
β
hB
IBiz
Bi z
humerus kH Brach kH Biz
Fig. 1. Schematic illustration of the elbow joint with the flexor muscles M. biceps brachii, M. brachioradialis and their effective moment arms hBiz ; hBrach ; respectively. The extensor muscle M. triceps brachii is considered in the model but not shown in this figure.
The total torque Tðb; oÞ is the sum of the torques of three muscle and two gravitational forces, each of these torques is given as the product of the force and the corresponding so-called effective moment arm. The activation levels of both flexor muscles are assumed to be equal and are denoted by Eflex : Following Murray et al. (2000) we set Flext ðbÞ 1; i.e. for the extensor muscle we neglect the dependency of the force on the length (cf. Appendix A.2). Using Eq. (1) for the three muscle forces we obtain the following expression for the
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
total torque: Tðb; oÞ ¼ Eext Hext ðhext ðbÞ oÞhext ðbÞ þ Eflex FlBiz ðbÞHBiz ðhBiz ðbÞ oÞhBiz ðbÞ þ Eflex FlBrach ðbÞHBrach ðhBrach ðbÞ oÞhBrach ðbÞ m u þ ml : ð4Þ
cos b g l 2 Here h ¼ hðbÞ denotes the so-called effective moment arm and we think of any factor Fl as being a function of b and of any factor H as being a function of the contraction velocity of the muscle h o: For explicit expressions of these functions, cf. Appendix A. For abbreviation, we finally set 1 Fðb; oÞ :¼ Tðb; oÞ: ð5Þ J Thus, the basic equation of motion (3) turns out to be of the form ( b’ ¼ o; ð6Þ o ’ ¼ Fðb; oÞ: From now on we force the right angle b ¼ p=2; o ¼ 0 to be an equilibrium point of system (6) by setting Eflex ¼ Eflex ðEext Þ for Eext A½0; Emax ; explicitly given by Eflex ðEext Þ ¼
hext ðp=2ÞHext ð0Þ Eext : hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
ð7Þ Lemma 1. (i) Relation (7) is sufficient and necessary for b ¼ p=2; o ¼ 0 to be an equilibrium point of system (6). (ii) Eext ; Eflex A½0; 1 is equivalent to Eext A½0; Emax ; ( Emax :¼
where
1
if Eflex ð1Þp1;
1 Eflex ð1Þ
otherwise:
ð8Þ
ð9Þ
Proof. In order to obtain an equilibrium, the right-hand side of the ODE system (6) has to be zero. Thus, o ¼ 0 and Fðb; 0Þ ¼ 0: This implies (i).
k :¼
117
positive. Thus, for the maximal activation either Eflex or Eext is equal to 1. We consider two cases: If the corresponding value Eflex for Eext ¼ 1 satisfies Eflex ð1Þp1; then 1 is the maximal value
1 for Eext : Otherwise, Eflex (1) is the maximal value for Eext : & Note that for Eflex ¼ Eflex ðEext Þ with relation (7), ðp=2; 0Þ is an equilibrium state for any mass ml of the load. This relies on the fact that Fðp=2; 0Þ does not depend on ml since cosðp=2Þ ¼ 0:
3. Results In this section, we state a number of theorems. The proofs are given in Appendix B. In Theorem 1, we determine the stability of the equilibrium ðp=2; 0Þ in terms of ml and Eext : It turns out that this equilibrium is asymptotically stable provided that the mass of the load ml is sufficiently small and—under suitable assumptions—the activation Eext is sufficiently large. In Theorem 2, we give a formula for a critical mass mcrit of the load which has the following meaning: for masses ml omcrit the equilibrium state ðp=2; 0Þ is asymptotically stable for sufficiently large activation levels and for masses ml > mcrit the equilibrium state ðp=2; 0Þ is unstable for all activation levels. In additional theorems, we study the occurrence of other equilibria at angles bap=2: In Theorem 3, we assume maximal activation, i.e. Eext ¼ Emax : We determine all equilibrium points ðb; 0Þ and their stability in terms of b and ml ; and we draw a corresponding bifurcation diagram. Theorem 3 is a special case of Theorem 4 where we, more generally, let Eext A½0; Emax and determine all equilibrium points ðb; 0Þ and their stability in terms of b; ml and Eext : Finally, we explain in Theorem 5 how unstable equilibria bound the basin of attraction of an asymptotically stable equilibrium in some sense. In an example, we illustrate the results using the data of a person given in Appendix C. Before we state the first theorem, let us define the following quantity k: Here, 0 denotes the derivative.
HBiz ð0Þ½h0Biz ðp=2ÞFlBiz ðp=2Þ þ hBiz ðp=2ÞFl0Biz ðp=2Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ HBrach ð0Þ½h0Brach ðp=2ÞFlBrach ðp=2Þ þ hBrach ðp=2ÞFl0Brach ðp=2Þ h0 ðp=2Þ
ext : þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ hext ðp=2Þ
To prove (ii) we will determine the maximal value of Eext such that both Eext and Eflex attain values in ½0; 1: Note that the factor in Eq. (7) is positive since hBiz ðp=2Þ; hBrach ðp=2Þo0 while all other quantities are strictly
ð10Þ
We now introduce f to denote the rate of change of F with respect to the angle b: We consider f as a function of the extensor activation Eext and the mass of the load ml : More precisely, we
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
118
define p 1 Hext ð0Þ
kE h fðEext ; ml Þ ¼ 1 ext ext 2 ð3 mu þ ml Þl 2 m i u þ ml : ð11Þ þgl 2 Note that fðEext ; ml Þ ¼ Fb ðp=2; 0Þ; where Fb denotes the derivative of F with respect to b and where we fix the activations Eext and Eflex such that the right angle is an equilibrium, i.e. we set Eflex ¼ Eflex ðEext Þ cf. Eq. (7). h
Theorem 1. With the notations in Eqs. (10) and (11) we have: (i) If fðEext ; ml Þ > 0; then the equilibrium ðp=2; 0Þ is a saddle and thus unstable. (ii) If fðEext ; ml Þo0; then the equilibrium ðp=2; 0Þ is asymptotically stable. (iii) ð13 mu þ ml Þl 2 fðEext ; ml Þ is strictly monotonically increasing with respect to ml :
We have limb-p=2 ml ðbÞ ¼ mcrit : Thus, adding the trivial branch of equilibria b ¼ p=2 for all ml X0 we can draw a complete bifurcation diagram (cf. the example at the end of this section and Fig. 2). Moreover, the trivial branch is asymptotically stable for ml A½0; mcrit ) and unstable for ml > mcrit : The bifurcating branch near ðmcrit ; p=2Þ is unstable for momcrit and asymptotically stable for ml > mcrit and again changes its stability at angles b0 with ðml Þb ðb0 Þ ¼ 0: There are no periodic orbits. In the following theorem, we allow Eext to attain any fixed value in the interval ½0; Emax : Theorem 4. Assume k > 0 and let Eext A½0; Emax : Then all equilibria ðb; 0Þ of system (6) with bap=2 are given by the function
ml ðb; Eext Þ :¼
p hBiz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBrach ðbÞFlBrach ðbÞHBrach ð0Þ Eext Hext ð0Þ mu hext ðbÞ hext
: cos b g l 2 hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ 2 (iv) fðEext ; ml Þ is strictly monotonically decreasing with respect to Eext ; if and only if k > 0: Theorem 2. Assume that k > 0: Let Eext Hext ð0Þhext ðp=2Þ mu k ; mcrit ðEext Þ :¼ gl 2 mcrit :¼ mcrit ðEmax Þ ¼
Emax Hext ð0Þhext ðp=2Þ mu k : gl 2
ð12Þ ð13Þ
Then we have (i) If ml omcrit ðEext Þ; then ðp=2; 0Þ is asymptotically stable. (ii) If ml > mcrit ðEext Þ; then ðp=2; 0Þ is unstable. (iii) For loads with mass ml > mcrit the equilibrium ðp=2; 0Þ is unstable for each activation Eext A½0; Emax : (iv) For loads with mass ml omcrit there is an activation E* ext A½0; Emax such that the equilibrium ðp=2; 0Þ is asymptotically stable for all Eext AðE* ext ; Emax : In particular, the equilibrium is asymptotically stable for Eext ¼ Emax : Theorem 3. Assume k > 0 and set Eext ¼ Emax : Then all equilibria ðb; 0Þ of system (6) with bap=2 are given by the function ml ðbÞ :¼ p Emax Hext ð0Þ hext ðbÞ hext cos b g l 2
hBiz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBrach ðbÞFlBrach ðbÞHBrach ð0Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
mu : 2
ð14Þ
ð15Þ
We have limb-p=2 ml ðb; Eext Þ ¼ mcrit ðEext Þ: Thus, adding the plane of trivial equilibria b ¼ p=2 for ml X0 and Eext A½0; Emax we can draw a complete bifurcation diagram. Moreover, the trivial branch is asymptotically stable for ml A½0; mcrit ðEext ÞÞ and unstable for ml > mcrit ðEext ). The bifurcating branch near ðmcrit ðEext Þ; p=2Þ is unstable for momcrit ðEext Þ and asymptotically stable for ml > mcrit ðEext Þ and again changes its stability at angles b0 with ðml Þb ðb0 ; Eext Þ ¼ 0: There are no periodic orbits. For an example of the situation of the following Theorem 5 cf. the example at the end of this section and Fig. 3. Theorem 5. Assume that ðbs ; 0Þ is an asymptotically stable equilibrium of system (6). Assume additionally that bu oð>Þ bs is an equilibrium such that Fðb; 0Þo0 ð> 0Þ for all bobu ðb > bu Þ and bA½0; p (i.e. in particular that ðbu ; 0Þ is the only equilibrium for boð>Þ bs and bA½0; p; and bu is unstable). Then the set fðb; oÞ j boð>Þ bu ; opðXÞ 0g does not belong to the basin of attraction of ðbs ; 0Þ: Example. As an example we use the data of the person given in Appendix C. The value of k (cf. Eq. (10)) is in this case k ¼ 1:14 > 0: We have Emax ¼ 1; because of Eflex ð1Þ ¼ 0:50p1: We calculate mcrit ¼ 15:74 kg using Eq. (13) of Theorem 2. We assume maximal activation
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
2 +
( mcrit , π2− )
1.5
angle β [rad]
(m*, β *) 1
119
b E1:298 is asymptotically stable. For masses ml > m we have no stable equilibrium for activations Eext ¼ Emax ; Eflex ¼ Eflex ðEext Þ: Fig. 3 shows the behavior of solutions as an illustration of Theorem 5 with ml ¼ 26 kg and Eext ¼ Emax : For these parameter values ðbs ; 0Þ with bs ¼ 1:403 is an asymptotically stable equilibrium. It is surrounded by the two unstable equilibria ðbu1 ; 0Þ and ðbu2 ; 0Þ with bu1 ¼ 1:152 and bu2 ¼ p=2: They mark the border of the basin of attraction of ðbs ; 0Þ in the sense of Theorem 5.
+ 0.5
4. Discussion
0
5
10
15
20
25
30
mass of the load ml [kg]
Fig. 2. All equilibria and their stability for Eext ¼ Emax depending on the mass ml of the load. Dashed lines denote unstable, straight lines stable equilibria. The signs denote the sign of Fðb; 0Þ in the corresponding regions of the ðb; ml Þ-plane, cf. the proof of Theorem 3 in Appendix B.2. The equilibrium ðp=2; 0Þ is asymptotically stable for ml omcrit : For mcrit oml om the bifurcating equilibrium is asymptotically stable.
3 Φ (β,0) 2 1 1
1.1
1.2
1.3
1.4
0
1.5
1.6
1.7 β
-1 Ð2 ω
. β=ω > 0
Ð3 β Ð4 . β =ω < 0
Fig. 3. Upper diagram: Fðb; 0Þ with ml ¼ 26 kg and Eext ¼ Emax ; lower diagram: the three equilibria ð1:152; 0Þ (unstable), ð1:403; 0Þ (asymptotically stable) and ðp=2; 0Þ (unstable) and the behavior of solutions of system (6) in the b–o plane.
Eext ¼ Emax and draw the bifurcation diagram in Fig. 2 using Theorem 3. The equilibrium ðp=2; 0Þ is asymptotically stable for masses smaller than mcrit : For masses between mcrit and the mass m E28:33 kg an equilibrium point at a smaller angle than p=2 bounded below by
The sign of the quantity k is crucial for most theorems. If k is negative, then for all positive loads ml the equilibrium ðp=2; 0Þ is unstable for any activation Eext by Theorem 1, since then (cf. Eq. (11)) h p 1 fðEext ; ml Þ ¼ 1
kE h Hext ð0Þ ext ext 2 ð3 mu þ ml Þl 2 m i u þ gl þ ml > 0: 2 Note that if k > 0; increasing the activation has a stabilizing effect, while for ko0 it has a destabilizing effect. We consider the influence of the different terms on the sign of k (cf. the definition of k in Eq. (10)). Since the slope of hext ; hBiz and hBrach is positive at p=2; it is necessary to have force–length relations FlBiz and FlBrach with a large positive slope at p=2 in order to obtain k > 0 (recall that hBiz ðp=2Þ; hBrach ðp=2Þo0Þ: Hence, for a realistic model it is essential to determine the force– length relation properly. In particular, a constant force– 0 length relation for the flexors would imply FBiz ðp=2Þ ¼ 0 FBrach ðp=2Þ ¼ 0 and thus ko0: On grounds of such a model the equilibrium at ðp=2; 0Þ would be unstable for any load. Moreover, regarding Eq. (13), the value of k is responsible for the critical mass mcrit you still can hold in a stable way. Hence, in order to increase the critical mass, besides augmentation of the isometric forces through training it could be important to improve the muscle–length dependency. From the bifurcation diagram we see that in our example loads with masses ml > mcrit can be held in a stable way with constant muscle activations, but at smaller angles than p=2: If the mass is augmented continuously, e.g. by filling a glass in the hand of the person with water, the equilibrium point at p=2 becomes unstable when the mass of the load exceeds mcrit : Then the elbow will—without intervention of the brain or reflexes—leave the unstable equilibrium position at p=2 (dashed line) and move to the stable position at a smaller angle (straight line).
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
120
By the bifurcation diagram, loads with mass ml > mcrit can be held in a stable way up to a certain mass m > mcrit and at a certain angle bap=2 (cf. the bifurcating branch of equilibria in Fig. 2). This shows that the right angle is not the optimal angle to hold large loads in a stable way. The question arises whether one can also hold masses beyond the mass m at an angle of 90 : In fact, as experiments have shown, it is possible to hold much higher masses, but not with constant activation of the muscles. This is in agreement with our analysis which showed that other stabilizing mechanisms are necessary to hold these masses in a stable way. In fact, one holds these heavy loads by switching on and off the muscle activations which can be shown by recording the electrical nerve stimulation of the muscles using electrodes (EMG). This situation of time-dependent activations is not included in our model. Stability is thus obtained through a combination of the self-stabilizing properties of the musculoskeletal system studied in this paper on the one hand, and the reflexes and the control through the brain on the other hand.
Appendix A. Derivation of the torque In this appendix, we give a detailed derivation of the torques corresponding to the gravitational and the muscle forces; the total torque Tðb; oÞ will be the sum of these torques. In general the torque is defined by the product of the length of the vector from the point of joint to the point of action of the applied force and the component of the force orthogonal on the aforementioned vector with appropriate sign. In our case, however, it will be simpler to determine the torque by the equivalent formula
l cosðp bÞ ¼ l cos b (cf. Fig. 1). The effective moment arm of the gravitational force acting on the ulna is hu ¼ 12l cos b due to the integration over the length of the ulna. Hence, m u Tðb; oÞ ¼ cos b g l þ ml 2 þ hBiz FBiz þ hBrach FBrach þ hext Fext : ðA:1Þ In the next sections, we will consider the torques corresponding to the muscle forces. In Section A.1, we determine the effective moment arms hBiz ; hBrach and hext and in Section A.2, we derive the corresponding forces FBiz ; FBrach and Fext : In Section A.3, finally, we will summarize the previous sections giving a formula of the total torque. A.1. The effective moment arms We calculate the effective moment arms h of each muscle. As mentioned before, we restrict ourselves to the most important muscles, namely the two flexor muscles M. biceps brachii and M. brachioradialis, and the extensor M. triceps brachii. We start with the two flexor muscles. Flexors: In Fig. 1 we marked the relevant quantities. The distances of the insertions of the muscles M. biceps brachii and M. brachioradialis to the elbow joint are denoted by kUBiz ; kHBiz ; kUBrach and kHBrach ; respectively. We calculate the effective moment arms hBrach and hBiz : Their sign will be negative. Using the sine-theorem, we get kHBrach sin a ¼ sin b ; ðA:2Þ lBrach sin g ¼ sin b
kHBiz : lBiz
ðA:3Þ
We conclude (note that the sign of hBrach and hBiz is negative)
T ¼ hF ; where F denotes the total strength of the force and h the so-called effective moment arm. It is given in the following way: draw a line which includes the point of the action of the force along the vector of the action of the force. Then there is a unique vector from the point of joint to a point of this line which is orthogonal on this line. h is 7 the length of this vector. The sign is negative, if this vector points in direction of decreasing b and positive otherwise. The total torque is the sum of the torques corresponding to the different forces Tðb; oÞ ¼ hu Fu þ hl Fl þ hBiz FBiz þ hBrach FBrach þ hext Fext : The gravitational forces are the force Fu ¼ g mu acting on each point of the ulna, and the force Fl ¼ g ml acting on the load. The effective moment arm of the gravitational force acting on the load is hl ¼
kHBrach kUBrach ; lBrach ðA:4Þ kHBiz kUBiz ðA:3Þ ¼ sin g kUBiz ¼ sin b : lBiz ðA:2Þ
hBrach ¼ sin a kUBrach ¼ sin b hBiz
The lengths of the muscles are given by the cosinetheorem as follows: lBrach ðbÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ¼ kHBrach þ kUBrach
2kHBrach kUBrach cos b;
ðA:5Þ
lBiz ðbÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 þ kUBiz
2kHBiz kUBiz cos b: ¼ kHBiz
ðA:6Þ
Extensor: In order to calculate the effective moment arm of the extensor we cannot give a geometric argumentation as in the case of the flexors. Instead, we
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
follow Gerbeaux et al. (1996) who approximated the effective moment arm (which they call geometry function) by a polynomial of third degree, namely hext ðbÞ ¼ l½0:0993 þ 0:0371ðp bÞ 0:0524ðp bÞ2 þ 0:0147ðp bÞ3 :
121
160 ð88p=9Þ: Hence, p p zBiz lBiz ¼ 0:45 ) c1 lBiz þ c2 ¼ 0:45; 3 3
ðA:7Þ zBiz
8p 8p lBiz ¼ 1 ) c1 lBiz þ c2 ¼ 1; 9 9
A.2. The muscle forces from which we can conclude Now we consider the muscle forces FBiz ; FBrach and Fext : Each of them depends on the activation of the muscle, the length of the muscle and the velocity of contraction of the muscle. We write the forces as the product of three functions, namely the activations Eflex ; Eext ; respectively, the force–length functions FlBiz ðbÞ and FlBrach ðbÞ; and the velocity–length functions HBiz ðb; oÞ; HBrach ðb; oÞ; Hext ðb; oÞ; respectively, which are also called the Hill functions. Note that we set the force– length function of the extensor Flext ðbÞ :¼ 1: FBiz ðb; oÞ ¼ Eflex FlBiz ðbÞ HBiz ðhBiz ðbÞ oÞ;
zBiz ðxÞ ¼
0:55 lBiz ð8p=9Þ lBiz ðp=3Þ
8p x lBiz þ 1: 9
A similar formula holds for zBrach : Hence, FlBiz ðbÞ ¼
ðA:8Þ
1:1 arctan½80ðzBiz ðlBiz ðbÞÞ 1:14Þ9 p þ 7ðzBiz ðlBiz ðbÞÞ 0:57Þ þ 0:55;
FBrach ðb; oÞ ¼ Eflex FlBrach ðbÞ HBrach ðhBrach ðbÞ oÞ; ðA:9Þ FlBrach ðbÞ ¼ Fext ðb; oÞ ¼ Eext Hext ðhext ðbÞ oÞ:
ðA:10Þ
Activation: We assume that both flexor muscles are activated at the same activation level EBiz ¼ EBrach ¼: Eflex : The activation of the muscle is given as the amount of the maximal activation and will be assumed to be constant in our model. Thus, the activations have to fulfill Eflex ; Eext A½0; 1: Force–length relation: The length-dependency of the (maximal) force which can be generated by the muscle has been studied by Murray et al. (2000) analysing dead bodies. Depending on the numbers of overlapping filaments the force is maximal at a medium length and decreases both if the muscle is longer and if it is shorter than that (cf. Gordon et al., 1966). The curve is the same for all muscles but the region in which the muscle works is different for different muscles. For the extensor, the force–length function is approximately 1 and thus we set Flext ðbÞ ¼ 1: For both flexors (M. biceps and M. brachioradialis), we model the function given in Murray et al. (2000) (at the left-hand region) by Fl ðzÞ :¼
1:1 arctan½80ðz 1:14Þ9 þ 7ðz 0:57Þ p þ 0:55;
ðA:11Þ
where z denotes the normalized muscle length. We assume zðxÞ which maps the total muscle length x to the normalized muscle length to be a linear relation of the form zBiz ðxÞ ¼ c1 x þ c2 : Following Murray et al. (2000) we postulate that the flexors attain 45% of their optimal length at an angle of 60 ð8p=3Þ and 100% at
ðA:12Þ
ðA:13Þ
1:1 arctan½80ðzBrach ðlBrach ðbÞÞ 1:14Þ9 p þ 7ðzBrach ðlBrach ðbÞÞ 0:57Þ þ 0:55: ðA:14Þ
Force–velocity relation: the Hill function: A common model for the dependency of the muscle force on the velocity in the literature is the so-called Hill force– velocity relation. Hill (1938) only studied the concentric part, i.e. movements, where the muscle is shortened. We, however, will need a function which is also defined for the excentric part and is C 2 ; and therefore we suggest the following Hill function: Hðvm Þ 8 c > > < vm þ b a ¼ B C > > þ :A þ vm D ðvm DÞ2
for vm X0; for vm o0;
ðA:15Þ
where a; b; c and A; B; C; D are positive parameters. The Hill functions of both the flexors and the extensor are of this form with different parameters. The contraction velocity vm is given by the derivative of the negative muscle length lm with respect to the time t: vm ¼
dlm : dt
The value Hð0Þ is called the isometric force Fiso :¼ c=b a: For vm - N the muscle can generate a stronger force, i.e. H tends to Fiso cexz ; where cexz > 1: In general, one assumes cexz ¼ 1:5: Together with the requirement that HAC 2 ðR; RÞ; the four parameters
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
122
A; B; C; D are determined, if a; b; c and cexz are given.
Dempster, 1955):
A ¼ cexz Fiso > 0; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b3 D ¼ b þ b2 þ Fiso ðcexz 1Þ > 0; c c B ¼ ð3b þ 2DÞ Fiso ðcexz 1Þ 2D 2 > 0; b c > 0: C ¼ D Fiso ðcexz 1ÞðD bÞ þ 2D b
mu ¼ 0:022m:
The derivative of the Hill function is strictly negative for all vm AR: For vm X0; this follows from c=ðvm þ bÞ2 o0: For vm o0 we have to show that ðBðvm DÞ þ 2CÞ=ðvm DÞ3 o0: It is sufficient to show 2CoBD: This follows by a straightforward calculation using the above definitions of B; C and D: In order to derive a relation between the muscle velocity vm and the velocity o; we use the following equality: the power of the muscle Fext vmext equals the power Text o: Thus, we have Fext vmext ¼ Fext hext o; vmext ¼ hext o:
ðA:16Þ
For the flexor muscles we have in an analogous way
Sust (1996) developed a method to determine the parameters of the Hill functions for individual persons through maximal voluntary contraction. We used this method to determine the parameters a; b; c of the Hill functions, and we additionally assumed that they are equal for both flexor muscles M. biceps brachii and M. brachioradialis. Appendix B. Proofs B.1. Stability of the equilibrium ðp=2; 0Þ Proof of Theorem 1. We denote the derivatives of F with respect to b and o by Fb and Fo ; respectively. The Jacobian of the right-hand side of the system of ordinary differential equations (6) is ! 0 1 ðB:1Þ Fb ðb; oÞ Fo ðb; oÞ
vmBiz ¼ hBiz o;
ðA:17Þ
and its eigenvalues at the equilibrium ðb; 0Þ are qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l1;2 ¼ 12 ðFo ðb; 0Þ7 Fo ðb; 0Þ2 þ 4Fb ðb; 0ÞÞ:
vmBrach ¼ hBrach o:
ðA:18Þ
Note that
A.3. The total torque We have described the muscle forces by Eqs. (A.8)– (A.10). Thus, the total torque is given by (cf. Eq. (A.1)) Tðb; oÞ ¼ hext Fext þ hBiz FBiz þ hBrach FBrach m u þ ml
cos b g l 2 ¼ Eext hext ðbÞHext ðhext ðbÞ oÞ þ Eflex hBiz ðbÞFlBiz ðbÞHBiz ðhBiz ðbÞ oÞ þ Eflex hBrach ðbÞFlBrach ðbÞHBrach ðhBrach ðbÞ oÞ m u þ ml ;
cos b g l 2 which corresponds to Eq. (4). As an example the data of a specific person are given in the Appendix C. The lengths of the ulna and the humerus are taken from those data, and we assume l ¼ lulna : The insertions of the muscles are determined following the studies of Wank (1996) kHBrach :¼ 0:104 lhumerus ; kHBiz :¼ lhumerus ; kUBrach :¼ lulna ; kUBiz :¼ 0:186 lulna : The mass of the arm is determined by the following formula through the total mass m of the person (cf.
ðB:2Þ
Fo ðb; oÞ 1 0 ðhext ðbÞ oÞh2ext ðbÞ ¼ ½Eext Hext J 0 þ Eflex FlBiz ðbÞHBiz ðhBiz ðbÞ oÞh2Biz ðbÞ 0 ðhBrach ðbÞ oÞh2Brach ðbÞ þ Eflex FlBrach ðbÞHBrach
ðB:3Þ
is strictly negative for all values of b and o; because all components are strictly positive except for the deriva0 0 0 tives of the Hill functions Hext ; HBiz and HBrach which are strictly negative (cf. Appendix A). Thus, there is at least one eigenvalue with negative real part. The second eigenvalue has negative real part if and only if Fb ðb; 0Þo0: Depending on the sign of Fo ðb; 0Þ2 þ 4Fb ðb; 0Þ we have either real or complex eigenvalues which means the equilibrium is either a node or a sink. The stability of ðp=2; 0Þ is determined by the sign of Fb ðp=2; 0Þ which we denote by fðEext ; ml Þ—note that we have used Eq. (7) to derive Eq. (11). This proves (i) and (ii). (iii) and (iv) follow directly from Eq. (11). & Proof of Theorem 2. Defining mcrit ðEext Þ implicitly by fðEext ; mcrit Þ ¼ 0; we have fðEext ; ml Þo0 ð> 0Þ for ml oð>Þ mcrit ðEext Þ: This relies on the specific dependency of f on ml and Eext which was shown in Theorem 1. This proves (i) and (ii). To prove (iii) we have ml > mcrit ¼ mcrit ðEmax Þ Xmcrit ðEext Þ for all Eext A½0; Emax and thus by (ii) the equilibrium ðp=2; 0Þ is unstable. For (iv) note that for
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
ml omcrit there is one and only one value E* ext A½0; Emax such that ml ¼ mcrit ðE* ext Þ since mcrit ð0Þo0; mcrit ðEmax Þ ¼ mcrit > ml and mcrit ðEext Þ is strictly monotonically increasing with respect to Eext : Thus, (iv) follows by (i). &
B.2. The bifurcating branch of equilibria and the basin of attraction Proof (of Theorems 4 and 3). We prove Theorem 4. Theorem 3 follows from Theorem 4 setting Eext :¼ Emax : The formula for the mass is derived from the condition Fðb; 0Þ ¼ 0 (cf. also Eqs. (5) and (7)). The # limit is calculated using l’Hospital: lim ml ðb; Eext Þ Eext Hext ð0Þ ½h0ext ðbÞ ¼ lim b-p=2 sin b g l
b-p=2
hext ðp=2Þ hBiz ðp=2ÞFlBiz ðp=2ÞHBiz ð0Þ þ hBrach ðp=2ÞFlBrach ðp=2ÞHBrach ð0Þ
ðh0Biz ðbÞFlBiz ðbÞHBiz ð0Þ þ hBiz ðbÞFl0Biz ðbÞHBiz ð0Þ þ h0Brach ðbÞFlBrach ðbÞHBrach ð0Þ þhBrach ðbÞFl0Brach ðbÞHBrach ð0ÞÞ
mu
2
¼ mcrit ðEext Þ ðcf: Eqs: ð10Þ and ð12ÞÞ: We calculate the trace of the Jacobian, which is Fo ðb; oÞo0 (cf. Eq. (B.1) and Appendix A). By the Bendixson criterion we can conclude that there is no periodic solution. Fix Eext A½0; Emax : To illustrate the argumentation also consider Fig. 2 which shows the bifurcation diagram for a certain person the data of whom are reported in Appendix C. The stability of the trivial branch at b ¼ p=2 is determined by the sign of Fb ðb; 0Þ; which is negative for ml omcrit ðEext Þ and which is positive for ml > mcrit ðEext Þ; cf. Theorem 2. Hence, we can determine the sign of Fðb; 0Þ for fixed ml near the trivial branch. Since we know that we marked all zeros of Fðb; 0Þ; we know the sign of this function for all b; because by the intermediate value theorem it cannot change its sign unless it reaches zero. Now we know the sign of Fðb; 0Þ everywhere and thus conclude the stability of the nontrivial branch. Near ðmcrit ðEext Þ; p=2Þ all different cases result in asymptotic stability for the bifurcating equilibria if ml > mcrit ðEext Þ; and instability if ml omcrit ðEext Þ: Assume that ðb0 ; 0Þ is an equilibrium and that the sign of Fðb; 0Þ changes from to þ at b0 for increasing b: Then ðb0 ; 0Þ is unstable. If Fðb; 0Þ changes its sign from þ to
at b0 for increasing b; then ðb0 ; 0Þ is asymptotically stable. In case of hyperbolic equilibria this is ensured through linearization, but even in the non-hyperbolic
123
case this is true as we show in Proposition 2 in Section B.3. Thus, we have determined the stability of all equilibria except for the critical ones. This proves Theorem 4. & Proof of Theorem 5. For illustration cf. Fig. 3. We only consider the case of an equilibrium bu obs ; a similar argumentation holds for the other case. We know that for op0 we have b’ ¼ op0 and solutions cannot cross the b-axis, because then o ’ ¼ Fðb; 0ÞX0 which contradicts the assumption. Thus, the solution with initial value ðb; oÞ; 0pbobu ; op0 cannot tend to ðbs ; 0Þ: &
B.3. A proposition using a center manifold We prove a general proposition which links the stability properties of equilibria of the first-order system b’ ¼ Fðb; 0Þ to the stability properties of corresponding ’ In equilibria of the second-order system b. ¼ Fðb; bÞ: case of a hyperbolic equilibrium one has to check the signs of the eigenvalues, but the proposition is also valid in the non-hyperbolic case, where we use a center manifold for the proof. Note that we need FAC 2 ðR2 ; RÞ in order to use the center manifold theorem. Proposition 2. Let FAC 2 ðR2 ; RÞ: Consider the two ordinary differential equations ( b’ ¼ o; ðB:4Þ o ’ ¼ Fðb; oÞ; b’ ¼ Fðb; 0Þ;
ðB:5Þ
b0 is an equilibrium of Eq. ðB:5Þ if and only if ðb0 ; 0Þ is an equilibrium of Eq. ðB:4Þ: If b0 is an unstable/asymptotically stable equilibrium of Eq. ðB:5Þ such that Fo ðb0 ; 0Þo0 holds, then the corresponding equilibrium ðb0 ; 0Þ of Eq. ðB:4Þ is unstable/asymptotically stable. Proof. ðb0 ; oÞ is an equilibrium of Eq. (B.4), if and only if o ¼ 0 and Fðb0 ; 0Þ ¼ 0: Now assume that b0 is an equilibrium and that Fo ðb0 ; 0Þo0: The Jacobian of the right-hand side of Eq. (B.4) is ! 0 1 : Fb ðb; oÞ Fo ðb; oÞ The eigenvalues at ðb ; 0Þ are l1;2 ¼12ðFo ðb0 ; 0Þ7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞÞ: Let us first assume that the equilibrium b0 of Eq. (B.5) is hyperbolic. If b0 is unstable with respect to Eq. (B.5), then Fb ðb0 ; 0Þ > 0: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi We can conclude Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0Þ > jFo ðb0 ; 0Þj and hence one eigenvalue is negative and the other positive. Thus, ðb0 ; 0Þ is a saddle and unstable. If b0 is asymptotically stable, then Fb ðb0 ; 0Þo0: Moreover,
ARTICLE IN PRESS 124
P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
either Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0Þo0 in which case we have two complex conjugate eigenvalues with real part 1 Fo ðb0 ; 0Þo0; or Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞX0 and then 2q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fo ðb0 ; 0Þ2 þ 4Fb ðb0 ; 0ÞojFo ðb0 ; 0Þj and both eigen values are real and negative. In both cases ðb0 ; 0Þ is asymptotically stable. Now assume that the equilibrium b0 of Eq. (B.5) is not hyperbolic. Then we have Fb ðb0 ; 0Þ ¼ 0: If the equilibrium is asymptotically stable, then Fðb0 þ x; 0Þ is negative for small x > 0 and positive for small xo0: If the equilibrium is unstable, then Fðb0 þ x; 0Þ is positive for small x > 0 or negative for small xo0: In the following we assume in the unstable case that Fðb0 þ x; 0Þ is positive for small x > 0: The other case can be dealt with similarly. We show more generally that the behavior of the one-dimensional system (B.5) is the same as the behavior on the center manifold of the twodimensional system (B.4). Since Fb ðb0 ; 0Þ ¼ 0 the eigenvalues of Eq. (B.4) satisfy l1 ¼ 0 and l2 o0: Thus, there is a one-dimensional stable and a one-dimensional center manifold. The stability of ðb0 ; 0Þ with respect to Eq. (B.4) is determined by the behavior on the center manifold. To study the dynamics we first determine the equation of the center manifold, then we study the behavior on the center manifold and finally we conclude the stability of the two-dimensional equilibrium. To determine the equation of the center manifold, we first apply a coordinate transformation in order to transform the system of differential equations into a suitable form. In fact, we choose the eigenvectors to be the basis of the new coordinate system. We set c :¼ Fo ðb0 ; 0Þo0: The Jacobian in this case ðFb ðb0 ; 0Þ ¼ 0Þ reads ! 0 1 : 0 c The eigenvalues are 0 with corresponding eigenvector ð1; 0ÞT ; and co0 with corresponding eigenvector ð1; cÞT : Hence, we define the following new coordinates: ! ! ! b b0 x 1 1c ¼ : y 0 1c o The system of differential equations in the new coordinates reads 8 1 > < x’ ¼ cy Fðb0 þ x þ y; cyÞ; c ðB:6Þ > : y’ ¼ 1 Fðb þ x þ y; cyÞ: 0 c ð0; 0Þ is the corresponding equilibrium of Eq. (B.6). By the theorem of the existence of a center manifold (cf. Carr, 1981, Theorem 1), we can describe a center manifold by a C 2 -function hðxÞ; such that locally near ð0; 0Þ a point ðx; yÞ lies on the center manifold if and only
if y ¼ hðxÞ: Moreover, hð0Þ ¼ h0 ð0Þ ¼ 0: To determine h we use the invariance of the center manifold, which implies y’ ¼ h0 ðxÞx’ and hence 1 Fðb0 þ x þ hðxÞ; chðxÞÞ c
1 ðB:7Þ ¼ h0 ðxÞ chðxÞ Fðb0 þ x þ hðxÞ; chðxÞÞ : c By Taylor’s formula there is a yA½0; 1 such that Fðb0 þ x þ hðxÞ; chðxÞÞ ¼ Fðb0 þ x; 0Þ þ Fb ðb0 þ x þ yhðxÞ; ychðxÞÞ hðxÞ þ Fo ðb0 þ x þ yhðxÞ; ychðxÞÞ chðxÞ:
ðB:8Þ
Plugging Eq. (B.8) into Eq. (B.7) we get ½Fðb0 þ x; 0Þ þ Fb ðb0 þ x þ yhðxÞ; ychðxÞÞ hðxÞ 1 þ Fo ðb0 þ x þ yhðxÞ; ychðxÞÞ chðxÞ ð1 þ h0 ðxÞÞ c ¼ ch0 ðxÞhðxÞ: If hðxÞ ¼ 0; then the above formula shows that Fðb0 þ x; 0Þ ¼ 0: But if b0 is asymptotically stable, then we have hðxÞa0 for all x in a neighborhood of zero. If on the other hand b0 is unstable, then we have hðxÞa0 for all xAð0; eÞ with an e > 0; by the assumption Fðb0 þ x; 0Þ > 0 for those x: Thus, hðxÞa0 holds in a neighborhood of zero, if b0 is asymptotically stable, and for all xAð0; eÞ; if b0 is unstable. Hence, we are allowed to divide both sides by hðxÞ and then let x-0; if b0 is asymptotically stable, and xr0; if b0 is unstable. Since hð0Þ ¼ h0 ð0Þ ¼ 0; we get with Fb ðb0 ; 0Þ ¼ 0 Fðb0 þ x; 0Þ ¼ c2 o0: ðB:9Þ lim x-0 hðxÞ Hence, Fðb0 þ x; 0Þ and hðxÞ have opposite signs for small jxj: The dynamics on the center manifold are governed locally near 0 by (cf. Carr, 1981, Theorem 1) 1 u’ ¼ chðuÞ Fðb0 þ u þ hðuÞ; chðuÞÞ: ðB:10Þ c We can expand the right-hand side using Eq. (B.8) again. There is a yA½0; 1 such that the right-hand side of Eq. (B.10) reads 1 chðuÞ ½Fðb0 þ u; 0Þ þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ hðuÞ c þ Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ chðuÞ 1 Fðb0 þ u; 0Þ ¼ c
þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ c hðuÞ
þ c Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ hðuÞ: Since by Eq. (B.9) Fðb0 þ u; 0Þ þ Fb ðb0 þ u þ yhðuÞ; ychðuÞÞ lim u-0 hðuÞ þ c Fo ðb0 þ u þ yhðuÞ; ychðuÞÞ ¼ c2 þ 0 þ c2 ¼ 0;
ARTICLE IN PRESS P. Giesl et al. / Journal of Theoretical Biology 228 (2004) 115–125
the sign of the right-hand side of Eq. (B.10) is determined by chðuÞ for small juj: Since co0; it has the opposite sign as compared to hðuÞ near u ¼ 0 and thus the same sign as Fðb0 þ u; 0Þ: Since the right-hand side of Eq. (B.10) has the same sign as Fðb0 þ u; 0Þ for small u as we just have shown, the stability properties of b0 with respect to Eq. (B.5) and u ¼ 0 with respect to Eq. (B.10) are the same. By Carr (1981, Theorem 2), the asymptotic stability/ instability on the center manifold (B.10) implies the asymptotic stability/instability of the two-dimensional system (B.4) and thus the theorem is proved. &
Appendix C. Data
Constant Description
Value
Unit
lulna lhumerus m aext
0.26 0.31 55 403.7752
m m kg N
bext cext cexzext aflex bflex
Length of the ulna Length of the humerus Mass of the person Hill function parameter (extensor) Hill function parameter (extensor) Hill function parameter (extensor) Excentric force increase factor (extensor) Hill function parameter (both flexors) Hill function parameter (both flexors)
0.2440 m/s 502.3032 W 1.5 21.4177 N 0.0233 m/s
cflex cexzflex
Hill function parameter (both flexors) Excentric force increase factor (both flexors)
125
30.7825 W 1.5
References Carr, J., 1981. Applications of Centre Manifold Theory. Springer, Berlin. Dempster, W.T., 1955. Space requirements of the seated operator (WADC-TR-55-159). Wright-Patterson Air Force Base, OH. Gerbeaux, M., Turpin, E., Lensel-Corbeil, G., 1996. MusculoArticular modelling of the triceps brachii. J. Biomech. 29, 171–180. Gordon, A.M., Huxley, A.F., Julian, F.J., 1966. The variation in isometric tension with sarcomere length in vertebrate muscle fibres. J. Physiol. 184, 170–192. Hill, A.V., 1938. The heat of shortening and the dynamic constants of muscle. Proc. Roy. Soc. (B) Lond. 126, 136–195. Murray, W., Buchanan, T., Scott, D., 2000. The isometric functional capacity of muscles that cross the elbow. J. Biomech. 33, 943–952. Sust, M., 1996. Modular aufgebaute deterministische Modelle zur Beschreibung menschlicher Bewegungen. In: Ballreich, R., Baumann, W. (Eds.), Grundlagen der Biomechanik des Sports. Ferdinand-Enke-Verlag, Stuttgart. Wagner, H., 1999. Zur Stabilit.at der menschlichen Bewegung. Dissertation, J.W. Goethe-Universit.at, Frankfurt a.M. Wagner, H., Blickhan, R., 1999. Stabilizing function of skeletal muscles: an analytical investigation. J. Theor. Biol. 199, 163–179. Wagner, H., Machalett, M., Blickhan, R., 2003. Bewegungslernen und . Attraktorstabilit.at—vorl.aufige Ergebnisse und mogliche Anwendungen. In: Riele, H. (Ed.), Biomechanik als Anwendungsforschung—Transfer zwischen Theorie und Praxis, Vol. 132. Czwalina Verlag, Hamburg. Wank, V., 1996. Modellierung und Simulation von Muskelkontrak. tionen. Sport und Buch StrauX GmbH, Koln.