A Purely Mechanical Theory of Muscular Contraction WILLIAM 0. WILLIAMS Department of Mathematics,

Carnegie Mellon University, Pittsburgh, Pennsylvania


Received II July 1989; revised 2 March 1990

ABSTRACT A simple mathematic model is constructed for the cross-bridge dynamics that govern muscular contraction, based on ideas introduced by Huxley and Simmons. The model differs from earlier ones in that it is based on distributions over time rather than displacement and uses a sharp disengagement rule. The constitutive functions of the model are set using classical experimental results and its predictions for the length-step experiment of Ford, Huxley, and Simmons are examined.



It has been observed for years that the contraction of muscles is effected by the relative sliding motions of the actin and myocin filaments that are the basic building blocks of the muscle fibers. It also is almost universally accepted that the relative motion of these molecular filaments is brought about by the cross-bridges: club-shaped parts of the myocin molecule that protrude from the myocin-constructed thick filaments and cyclically interact with the actin in the thin filaments as ATP is hydrolyzed (see, e.g., [81X Various mathematical models of muscular contraction have been constructed based on this simple concept. Huxley in [ll] introduced the idea of modeling the myofibrils as a continuum of active cross-bridges that create force through extension and contraction of an elastic tail anchored in the thick filament while the head of the cross-bridge is attached to the thin filament. He formulated the kinetics of attachment and detachment in a straightforward way, proposed rules for the coefficients governing these attachments and detachments, and deduced predicted behaviors of the muscle. Most of the mathematical work that has been done since could reasonably be described as correction and elaboration of Huxley’s work. In particular, the basic formulation that he used was modified by Eisenberg and coworkers (e.g., [5]), by Huxley and Simmons [12, 131, by the Pavia group, Comincioli and Torelli et al. (e.g., [3, 4]), by Pollack (e.g., [17]), and MATHEMATICAL


OElsevier Science Publishing 655 Avenue of the Americas,




Co., Inc., 1990 New York, NY 10010





by Lacker and Peskin [15, 201. In all of these the essential idea remains the same, in that the number of attached cross-bridges is modeled as distributed over the variable x, representing displacement of the cross-bridge from its rest or force-free position. This leads to some problems: The treatment of isometric contractions requires a distortion of the theory, and it does not lend itself naturally to the use of time-governed phenomena. Here I introduce a model in which the cross-bridge number is distributed over time, which I believe leads to a more natural formulation and should more naturally lend itself to those further elaborations of the model that are necessary to applications to in uiuo modeling. In [12, 131, Huxley and Simmons introduced an important conceptual modification of the Huxley model. Huxley’s original picture of the action of the cross-bridge required that it be possible that it attach to the actin filament while in a displaced position (allowing for initial generation of force) and that further development of force occurred due to mechanical transport of the myocin head by the relative velocity of the thick and thin filaments. Huxley and Simmons suggested a model in which the myocin head first attaches to the thin filament with its elastic tail in a relaxed position; rotation of the head about the head-tail junction then stretches the tail, creating an elastic force even in the absence of relative motion of the fibers. They assumed that the rotation was staged and introduced a kinetics based on rate-law transitions between two states. Here I extend their model to a continuous rotation. I choose to picture the myosin head as a cogwheel attached to an elastic element that is anchored to the thick filament (see Figure 1). Upon attachment (engagement of the cog with the rack attached to the thin filament), the cogwheel begins to rotate at a constant speed (corresponding to the maximal possible speed of contrac-


I FIG. 1.





tion), stretching the elastic element and/or transporting the rack and hence the thin filament. The cross-bridge remains attached, in this mechanical model, until it has completed its rotation or until the elastic element is forced beyond its limits in extension or in contraction. The fixed time of attachment allows modeling of hydrolyzation of a fixed quantum of fuel (ATP) per cross-bridge. Likewise, by allowing a fixed rule of detachment, I eliminate the stochastic detachment rule used in previous theories and simplify the model. I do not presume to claim that the simple mechanical description that I have sketched here and describe more fully below describes what actually occurs in muscle. Most fundamentally, it has not been clearly established experimentally that any sort of systematic rotation of myocin heads actually occurs during the contraction process [19]; beyond this, no model as elementary as this can reasonably be suggested to be an exact model of the interacting molecules in the muscle.’ Nonetheless, the model of Huxley and Simmons is a plausible and roughly physiological description of the muscle microstructure and can account for much of the observed behavior of muscle fibers. The mathematical embodiment of that theory that I present here I feel captures the essence of the process and offers a very realistic possibility of a model simple enough to figure in more complete models of entire muscle behaviors. In Section 2 I enlarge upon this sketch of the model and describe it mathematically, finally presenting the evolution equations that govern the system. After verifying in Section 3 that these equations are equivalent, for nonzero velocity, to a formulation like that of Lacker and Peskin, in Sections 4 and 5 I introduce constitutive equations, first specifying the detachment function and then using the Hill curve to determine the elastic response function. In Section 6 I discuss the length-step experiment of Ford et al. [6], using their experimental results to justify a form for the attachment function and then showing how the predictions of the theory compare to those results. Of particular interest is that the model predicts the two-stage exponential recovery found in the experiments without direct introduction of reaction rate equations. Section 7 has a discussion of the mathematical well-posedness of the evolution equations. In the Appendix, I show how to introduce a viscoelastic response in place of the elastic one and show how this affects the predictions for the length-step experiment.

‘A beginning toward a realistic model based on the rate analysis of the linked biochemical reactions involved is given in [211. The resulting model is (appropriately) complicated, so much so that it probably is not useful for embedding in a more general global model of muscle mechanics.

252 2.








We model a half-sarcomere; since a muscle filament consists of a series arrangement of such units and a muscle fiber of a family of parallel filaments, the general response of a fiber will be functionally the same, assuming that the myofibrils that make up the fiber are uniformly activated. Likewise, if we presume an entire muscle uniformly activated, then the model will describe the response of the entire muscle, within the admittedly narrow range where the muscle architecture does not interfere. Huxley and Simmons propose a clever mechanical model for the interaction of the two filaments that make up the half-sarcomere. The cross-bridges are pictured as rotating pinions attached by an elastic tail to the thick filament; these pinions then periodically engage the rack attached to the thin filament, as illustrated in Figure 1. They propose that the pinion flip between two possible positions (a three-tooth pinion) according to a rate law. Here we shall instead model the process as being a continuous one, with the pinion rotating counterclockwise at a steady rate from the time it engages until the time it detaches. The recovery period, which would represent another phase of the full chemical process of the ATP-ADP cycle, in this model would enter as a time of rerotation of the detached pinion. As suggested by recent studies (see [l]), it would be reasonable to suppose that in states of submaximal, or in particular zero, activation, a certain fraction of the cross-bridges remain attached to the thin filament in a passive mode (here, in the initial position without rotation); we will presume throughout that the system is in a state of maximal activation, so that we do not have to introduce this complication and consider only fully active cross-bridges. Looking at Figure 1, one sees that the rotation of the pinion is compensated by elastic extension of the attaching link or by translation of the rack or by a combination of the two. The steady angular velocity of the pinion we measure by the velocity of the teeth, and thus we can describe it as u,, the maximum velocity of contraction of the filament. If the velocity of the thin filament is u (measured as positive if the sense of motion is to the left in the figure), then the velocity of extension of the connecting link is u. We suppose that the process of rotation occupies a total time (Yand “?n that the pinion disengages within a shorter time only if the link is forced to exceed its maximal extension Y or is forced to shorten to less than its initial length.* Within the limitations imposed by these constraints we may con-

‘This is chosen as the simplest possibility. It could be pictured that the cross-bridge tail could accept some compression; in the Appendix we remark on the effect this might have on the predictions for the length-step experiment. Compression would not enter in any of the other experiments we discuss.



sider all possible velocities u either positive or negative. Finally, let us note that we shall replace Huxley and Simmons’ assumed linear elastic response by a nonlinear one. The linear elastic law does not allow easy replication of the classical hyperbolic force-velocity law of Hill without the introduction of elaborate demographics of the cross-bridge population. One may argue that full description of experimental behavior might require viscoelastic behavior in the model at some point. As we shall see, this model exhibits exponential decay to steady-state values without the addition of “parallel viscoelastic elements” or the series viscoelastic effect often introduced in whole muscle models to represent the tendons that secure the muscle. However, we find that the purely elastic model leads to some slightly anomalous behavior, which could be partly relieved by replacing the elastic element in the thin filament by a viscoelastic one. It is also at least appealing to prejudices that biological materials should always have a damped response. For these reasons, in the Appendix we will consider the effect of introducing an elementary viscoelastic response into the model. Our description of the demographics of the population of attached cross-bridges is similar to that used by Huxley and many others, except that we shall represent the number of attached cross-bridges as an integral in time rather than in displacement. Thus, if N(t) denotes the number of cross-bridges attached at time t we suppose that

with b identified as the time of attachment. Alternatively, we could follow the nomenclature of population dynamics and call b the birth time so that the elapsed time of attachment ’ a:=t-b

could be called the age of the cross-bridge. The natural domain of the function p, the number density, is {(t, b)lb < t}, but we extend it to all of the plane by letting it be zero outside this region. We expect p to be smooth (C’) except on a finite set of smooth curves, where it may suffer jump discontinuities. An obvious first example of such a curve is the attachment line z = b at which p may jump from zero to a new value:

Udb, b)ll ’ 0. Here we use the convention that the jump 11.1as t increases is positive if the quantity increases with t. It may also occur that sudden changes occur at various other times. In particular, we shall regard each cross-bridge as




having a fied

maximal time of attachment LY.This represents the physiological exhaustion of fuel available, and, in the mechanical model, completion of the allowable rotation of the pinion gear. Thus we suppose

p(t,b) = 0

if t>b+a.

But in addition we shall allow only a limited range of motion of the attached cross-bridges. Let us introduce the velocity of contraction u(t), which we presume to be a piecewise smooth function of time, the maximal velocity v,, and the corresponding displacement function y(t, b): y(t,b):=/bf[u,-v(T)]dr.


that there is a maximal displacement Y that cannot be by any cross-bridge, and that none are allowed to be carried into

we suppose


negative displacement.


dt,b) = 0

if y(t,b)


p(t, b) = 0

if y(t,b)

dt> 6) = 0

if t>b+a,



if y(t,b)


p(t, b) = 0

if y(t,b)




If we let G(t,c,b):=exp then the solutions


p(t,b) = cp(b)G(t,O,b)


and =f(l-


Each solution persists until it reaches a detachment of notation, let us set

curve. For convenience

D(t)={blt-a = f(lwith corresponding



N(t) = _/$ The force term transforms


Of course,


in a similar manner:


These equations were considered in [20] with a nonlinear elastic law determining the relation of p and y, and existence and uniqueness and a one-to-one relation between force and velocity were established. The use of the discontinuity attachment condition is replaced by Huxley and most of his followers by a stochastic attachment term similar to the term involving g, so the equation becomes




thus allowing attachment to occur at various values of displacement. A difficulty with this approach is that one must postulate a mechanism for displacement of the elastic element prior to attachment. (The hypothesis that such displacement is created by thermal agitation (see [ll]) leaves one still with the problem of explaining why this displacement is not symmetric about zero.) In the mechanical model we introduce here, we have chosen to follow Lacker and Peskin in introducing a jump condition for the attachment rule. 4.





By constitutive equations we mean those equations that specify the response of a particular material as opposed to the general equations that apply to all materials. Here these relations would include the specification of the attachment function f, the detachment function y, the attachment limits (Y and Y, and the form of the force-density versus displacement relation. The stochastic detachment term y is not directly measurable. Huxley took it to be constant and large for negative displacement and linear but small for positive displacement; this was motivated within his model by the requirement that detachment must occur even after motion ceases. In Lacker and Leskin’s modification, which is closer in spirit to our model here, g takes on a strong S shape, with high detachment rates far from zero displacement and very small ones near zero (they do not consider negative values of displacement). We achieve more or less the same effect as Lacker and Peskin with our abrupt detachment at a prescribed time or displacement, so it is not clear that g need even enter into the model. In some sense, the term is inconsistent with the purely mechanical model, and since its deletion does not seem to change the predictive ability of the model, at least for the simple experimental situations we consider, we shall omit it.3 Thus we shall confine our attention to the case y := 0.

This leads to a trivial equation

of evolution:

p(t,b) =P(b+,b)


pt = 0, so


31t might be expected that such a term would be of significance in conditions of fatigued muscle, when it is not reasonable that each cross-bridge could remain active for the entire time a, or when the muscle is subject to extreme oscillation.



and N(t)




The magnitude of the time limit (Y in principle could be obtained from knowledge of the chemical reaction that occurs in the cross-bridge cycle, and checked by looking at the time-course experiment of Section 6. Similarly, it should be possible in principle to obtain the limit on the extendability of the myocin tail from experimental observations. Such data do not seem to be available, so we deduce them from the force-velocity experiments discussed in Section 5. What we find there in particular is that in the normal course of contraction the time limit predominates over the displacement limit. 5.





In this section we investigate a general form for the force-density function. Without choosing particular values, we show how the Hill curve, a hyperbolic relation between force and velocity in steady contraction, dictates a form for our force-density function. Hill [9] discovered this functional relationship in his early experiments, and while it is no longer regarded as precise, it certainly should be accurate enough for the purposes of our model. We shall assume, in correspondence with Figure 1, that the force in the cross-bridge is generated by a nonlinear elastic response to displacement. There is some question, which cannot be fully answered by our study here, as to whether a viscous contribution should be present: While any biological material may be expected to show some rate-dependent behavior, it is not at all clear whether the damping effect observed in muscle response significantly reflects effects in the microstructure of the muscle or is solely a response of the tendons that connect the muscle to the bone. To maintain generality (and to illustrate that it is not difficult to do), we will illustrate in the Appendix that it is possible to replace the elastic response of the myosin cross-bridges by a viscoelastic one without major changes in the modeling process. Here we will assume that the specific force n- is engendered by the displacement through a nonlinear elastic law @-:

Steady Motion. The force-velocity experiments that we consider are based on a presumed steady-state condition: The velocity is assumed to be constant in time and the experiment conducted over a period in which the muscle fibers remain in a position of maximal overlap, so that the number




of active cross-bridges at each level of activation (measured by the time of attachment) can also be taken to be constant. For our model this means that the velocity should be taken to have been constant for all time: LJ(t) = ug and that p(t,,bJ

if t,-bb,=t,-bZ.


Clearly, there is an advantage attachment or age. We define

to reformulating

in terms of the time of

A(a) := p(b + a,~). The support of A clearly lies in [0, A], where A:=min{a,&},

and the total number

The assumption

of attached


is independent

that y is zero now leads to A’= 0, so A=A,

and hence N,, = A,A. Now, if we define M,:=lto represent

the number

of unattached

N, cross-bridges,

A, = f(Mo) and M,=l-L+(M,)da=l-Af(Mo),

we have

of time:




The displacement

y is now given by


=(,bm- uo) = (urn- uo)(t -b)

so that

where P= G, so that F is the potential corresponding to 6. It is instructive to examine the form taken by P in the case that e is linear. If e(y) = Ey, then F(y) = iEy2, and if A = (Y, then


PO= while if A = Y/(u,

- u,),




Thus in the case in which the cross-bridges complete their cycle, the response would be linear and decreasing in velocity, uo, while in the case in which stripping of the cross-bridges due to overextension occurs, the response would be increasing in uo, with an asymptotic constant value. Anticipating a general behavior (Hill’s hyperbola) as in Figure 2, we are led to assume that only in the case that u. is negative does stripping occur, that is, that in general A=a.

In terms of the mechanical model, this means that the rotational extension u,cy is less than Y. Now let us give a functional form to F, and hence to .+, by introducing the Hill equation. We deal with the situation where the muscle is in a



FIG. 2.

steady contraction

so that

~O=&~[(u,-%)4~ 0 m

while the classical Hill hyperbolic

PO =

law is

4 u, - uo) b+u, ’

u1 Q ug d Urn,

with u1 a negative value. Thus we deduce a form for F,

Cd2 F”(d>=B_dv


c :=

a Ao’Y






FIG. 3.

From this the nonlinear

elastic law is CB=




which has the general form shown in Figure 3. As mentioned, we will adduce a separate form for the response Pa in the case uc < ui. Since the same elastic law is presumed to govern, but now with the stripping rule (A = Y) in force, we find that

PO= um so the total picture




i (B-Y)2


for PO is as in Figure 4.

Remark 1. We have assumed that the steady density of attached crossbridges is independent of the velocity uo, which may be questioned. Lacker and Peskin assume that it does vary with uo. However, the present model is based on the concept that the attachment of the cross-bridges is governed purely by activation, and it is reasonable that the muscles in the experiments studied were fully activated. I do not know of data that establish




FIG. 4.

directly the variation of the number attached with the velocity,4 indeed, even measurements of the number of cross-bridges attached at a given time do not seem to be reliable [191. Remark 2. The presence of the stochastic detachment function y would affect the determination of these curves; assuming that the detachment is a function only of the time of attachment, y(t, 6) = P(t - 6), one finds

and it is not difficult to obtain a value for the elasticity function if P is known. To obtain a form for the function P from experimental evidence is not so easy, however. Lacker and Peskin use Hill’s energetic formula to determine its form; such data are somewhat questionable [lo], as the sources of heat in contraction are many and are not easily related to

4The only measurements that I know of are indirect. For example, Julian and Sollin’s [14] argument that N, = A,A should decrease with velocity is based on transient force measurements showing a decrease of stiffness with velocity. But this is dependent upon their particular constitutive assumptions. In the present model, with A, constant, the approximate stiffness at velocity ua is h,,(e[(u, - u,)A]l/(u, - ~a), which, since .z. is convex, always decreases as uc increases from 0 to u,.


cross-bridge turnover. One presumed virtue of the present is simple enough to evade such difficulties. 6.

265 model is that it


In this section we will consider modeling another classical experiment. Podolsky [16] introduced an experiment involving a sudden release of an isometrically held stimulated muscle; later experimenters [6] refined the experiment to create, instead of this force step, a length step: a nearly instantaneous reduction in length. An approximate reproduction of the resulting force-time plot is shown in Figure 5; it is characterized by a drop in tension concurrent with the length step followed by a two-stage recovery of tension. Such a situation can present a good test of our mathematical model, as it involves a type of behavior quite different from the characterizing one. In addition, it is of interest because it presents an opportunity to examine whether this model can produce such a damped response without the addition of the exponential decay due to the stochastic term or the addition of viscoelastic response or linear rate laws. Finally, it allows us to examine the question of the appropriate form for the attachment function f, which we have not yet considered in detail. Let us assume that the muscle is in a steady isometric state. The density of attached cross-bridges is A0 and has been so for a time period greater than (Y, and the velocity of contraction is zero. At time zero the muscle is


Time FIG.5.



allowed to contract instantaneously by an amount 6. (It is important to recognize that this is not a steady contraction, which would be governed by the maximum speed U, but instead is a release of elastic extension, which occurs instantaneously.) Referring to our picture of the mechanical model, we suppose that 6 does not exceed u,cy, the maximum extension of the cross-bridge tails, and recognize that the effect of the contraction is to cause detachment of those cross-bridges whose current extension is less than 6 and to reduce the extension of those whose extension is greater than 6.’ Thus the distribution of number density with time t and time of attachment b is as shown in Figure 6. Also graphed in Figure 6 is the density distribution as a function of “age” a and the displacement as a function of a, all at the time t = 0+ . Note the parameter (~a := 6/u,, which is of relevance in the later computations. Clearly the steady-state conditions are eliminated at this time, as the number of, unattached cross-bridges suffers a jump increase, so that we must recompute the number and the force. It is possible to do so analytically but is exceptionally tedious, so a numerical calculation was done. At this stage, it is necessary to discuss the form of the ‘attachment function f, as it becomes significant for the first time in our considerations. Consider first the balance equation. Let us reintroduce M as the number of unattached cross-bridges, so M(t):=l-N(t). Then, since y is zero, the equation

M(t) = l-

of balance


As noted earlier, the time-independent


of cross-bridges



db. MO will satisfy

or f(M,)

= q.

‘The actual length of extension of a single cross-bridge, of course, cannot compare to the extension S of the entire muscle. However, the filaments of the muscle are connected in series so that the displacements are additive and presumably the displacement 6 actually distributes over the entire collection of filaments. The model involves such an implicit addition of lengths and velocities throughout.


We wish to examine

the stability of this solution.

M(t):=M(t)-Mo=-/f and consider





the linearization

so that

m(t) = -/I



Pm(b) db.

m(t) = exp(At), we find the eigenvalue




The standard assumption for f, that it is linear, that is, that f(M) = PM, leads in particular to the requirement that /3 > 0. It is easy to see, though, that if /3 is positive there is no real solution to the eigenvalue equation and hence that oscillatory behavior predominates. This would lead to an oscillatory recovery in the length-step experiment, which in general is not observed; compare Figure 5, which is a reproduction of results from [6]. On the other hand, one can see that if p is negative and greater in magnitude than l/a, then the equilibrium point is not stable. Hence we are left with the conclusion that /3 must be negative and less in magnitude than l/a, which ensures that there is a real negative root to the eigenvalue equation. Since the uptake must be zero at M = 0, we are led to hypothesize a behavior as shown in Figure 7. Notice that we have chosen the form of f to ensure that there is a single equilibrium point MO for the system, in this case at Ma = 0.4. With this form for the attachment function and the form of the elastic law previously discussed, a computation was made of the response of the muscle to the length-step experiment. The values of the various parameters were chosen essentially arbitrarily, as it is rather difficult to find sufficient data for a single muscle to begin to deduce physiological values for them.6 In Figures 8 and 9 we illustrate the results of the computations for a typical case. (Note that this is an experiment with a large step, as the value of force, P, drops almost to zero.) In Figure 9, the correct sort of two-stage recovery of the force is predicted, with an exponential final recovery. It

FIG. 7.





FIG. 8.

should be emphasized that there is no extraneous damping introduced into the model or into the computation; the exponential behavior is purely a consequence of the dynamics of the cross-bridge demographics fundamental to the model. Let us examine these dynamics. The process, beginning with a

distribution of cross-bridges as in Figure 6, divides into three stages. For times between 0 and CY- LYE,the equation governing A4 is

M(t)=l-/-UoA,,db-jdf(M(b))db 1-a =l+h,,[t-(ru-a,,)]-jdf(M(b))db.

The same equation, less the A, term, governs between (Y- CQ and LY,while the period of normal recovery beginning at time cx is governed by the

6Values of the Hill parameters and of the maximum velocity of contraction are obtainable from the standard Hill experiments. As observed before, I do not know of experiments that establish the value of M,. Likewise the precise form of f would not be easy to obtain, although a value for p could be found from determination of the coefficient of decay of the force in the length-step experiment (see below), and this could suffice for the present use.















Time FIG. 9.






In this final stage of recovery note that

M’(t)=f(M(t))-f(M(t-a))--p[M(t)-M(t-a)], which yields an exponential solves

decay to M,, with an exponent

of p, where p

/A=-_p[ePLn-l]. Of course, this is just the eigenvalue calculation that was involved in the analysis of the behavior of f. In the final time period, P is governed by the equation p(t)




so, since 40) = 0, P’(t)=-,(u,&f(M(t-a))+/’


=/I e(~,(t-b))f’(M(b))M’(b)db. t--a






Since we can approximate

f’(M(~))M’(b) = P/-+fW - MO1 = P.f(M(b))> it follows that

and so P decays to its equilibrium value at the same exponential rate as A4 does. The exponent of decay depends upon LYand p through the eigenvalue equation. It should be noted that since p is approximately l/a, PLYis nearly constant, and hence the half-life of the decay is approximately linear in (Y. In the earlier time stages, the equations for P are straightforward. If 0 c t Q LY- ao, then P(t)=~[.7(c,,(u-a,,))-.Y(o,t)]+~~r(~~(t-h))f(M(b))db, m while the same equation governs in (Y- (~a < t < a, but without the first term. The prediction for the behavior of A4 is curious, perhaps, in that M rises, and hence the number of attached cross-bridges continues to drop, in the period 0 Q t 1, so that the response is initially asymptotic to the G = 1 elastic form. With this and the indicated P, we find the form for P predicted for the experiment: For t Q 0:


for O


stiff but form for




and for a

A purely mechanical theory of muscular contraction.

A simple mathematic model is constructed for the cross-bridge dynamics that govern muscular contraction, based on ideas introduced by Huxley and Simmo...
1MB Sizes 0 Downloads 0 Views