Non-axial Muscle Stress and Stiffness

Non-axial Muscle Stress and Stiffness

J. theor. Biol. (1996) 182, 59–84 Non-axial Muscle Stress and Stiffness G I. Z Department of Mechanical Engineering, Washington University...

491KB Sizes 14 Downloads 31 Views

J. theor. Biol. (1996) 182, 59–84

Non-axial Muscle Stress and Stiffness G I. Z Department of Mechanical Engineering, Washington University, St. Louis, MO 63130, U.S.A.

(Received on 27 November 1995, Accepted in revised form on 10 May 1996)

A generalization is developed of the classic two-state Huxley cross-bridge model to account for non-axial active stress and stiffness. The main ingredients of the model are: (i) a relation between the general three-dimensional deformation of an element of muscle and the deformations of the cross-bridges, that assumes macroscopic deformation is transmitted to the myofibrils, (ii) radial as well as axial cross-bridge stiffness, and (iii) variations of the attachment and detachment rates with lateral filament spacing. The theory leads to a generalized Huxley rate equation on the bond-distribution function, n(w, u, t), of the form 1 1n 1n + (D n) + m{lz Dzz + 2erDzr } = fn˜ − (f + g)n, 1w 1t 1u ur where the Dij are the components of the relative velocity gradient and r and n˜ are functions of the polar angle, u, and time that describe, respectively, the deformation of the myofilament lattice and the distribution of accessible actin sites (both of these functions can be calculated from the macroscopic deformation). Explicit expressions, in terms of n, are derived for the nine components of the active stress tensor, and the 21 non-vanishing components of the active stiffness tensor; the active stress tensor is found to be unsymmetric. The theory predicts that in general non-axial deformations will modify active axial stress and stiffness, and also give rise to non-axial (e.g., shearing) components. Under most circumstances the magnitudes of the non-axial stress and stiffness components will be small compared with the axial and, further, the effects of non-axial deformation rates will be small compared with those of the axial rate. Large transverse deformations may, however, greatly reduce the axial force and stiffness. The theory suggests a significant mechanical role for the non-contractile proteins in muscle, namely that of equilibrating the unsymmetric active stresses. Some simple applications of the theory are provided to illustrate its physical content. 7 1996 Academic Press Limited

in the fiber direction (Caputi et al., 1989). Another example is the myocardium. It is well known that the direction of cardiac muscle fibers rotates by approximately 180° through the thickness as one proceeds from epicardium to endocardium (Streeter et al., 1969), although the details of this structure are still under active investigation (Hunter & Smail, 1988). This ‘‘composite’’ construction of the myocardium implies that in any deformation of an element of the ventricular wall most fibers must experience non-axial strains. Thus it is important for understanding the function of muscular organs to know how non-axial deformations may affect the stresses in muscle tissue.

Introduction Although muscle is usually regarded as a uniaxial force generator, producing force and experiencing length changes in the fiber direction, it is often subjected to non-axial deformations. Changes in concentrations of osmoticallly active solutes can compress or expand the tissue. Even under normal physiological conditions the constraints of muscle architecture can impose non-axial strains on actively contracting muscle. For skeletal muscles with pennate architecture a change in overall muscle length will involve some shear strain in addition to normal strain 0022–5193/96/170059 + 26 $18.00/0

59

7 1996 Academic Press Limited

. . 

60

osmotic expansion and compression indicates that the axial muscle stress decreases markedly as the lateral myofilament spacing becomes either large or small, dropping to zero at lateral stretch-ratios of roughly 0.7 and 1.5, respectively (Howarth, 1958; Koch-Weser, 1963; Blinks, 1965; Gordon & Godt, 1970; Krasner & Maugham, 1984; Kawai & Schulman, 1985). Not only the non-axial strains, but also the non-axial strain rates can be expected to affect muscle stress. For example, consider a sample of muscle tissue subjected to a macroscopic time dependent shear strain that is transmitted to the myofilaments. The situation is illustrated in Fig. 1(b), which shows a single thick filament and surrounding thin filaments with which it is interacting. For the thin filaments ab and gh the motion appears to be that of axially lengthening muscle, as cross-bridges are being instantaneously stretched. Conversely, for ef and kl the motion appears to be that of axial shortening, with cross-bridges being pushed towards the M-line. Finally for thin filaments cd and ij the situation appears to be that of isometric muscle, where the thick and thin filaments remain stationary with respect to each other. These various relative sliding

At the other end of the geometrical scale, a study of non-axial deformations in muscle could contribute important information about the molecular mechanisms of contraction, the role of the non-contractile proteins in the sarcomere, the effects of fiber deformation vs. calcium activation, and similar questions. Some simple considerations of cross-bridge theory suggest that it is not unreasonable to expect that non-axial deformations may affect the stress in muscle. For example it is generally accepted that the myofilament lattice of intact muscle fibers changes shape during shortening or lengthening so as to maintain constant volume (Elliot et al., 1963). Based on this incompressibility condition, Fig. 1(a) shows, to approximately correct scale, the extent of lateral cross-bridge movement when a sarcomere is shortened from zero overlap to the length at which the isometric force drops to zero. The total displacement of the cross-bridge is approximately 18 nm which is large compared with a typical separation between the binding region on the cross-bridge and the actin filament. Such large relative motions could produce significant variations in the intermolecular forces, bonding rates, and stresses. Indeed a considerable body of experimental literature on the effects of s = 1.05 µm

s = 2.2 µm

s = 3.65 µm

(a)

. γ k

l

a

b

C d i j g

h

e

f (b)

F. 1. (a) Lateral motion of a cross-bridge, approximately to scale, as a function of sarcomere length. (b) Illustration of the effects of shearing on relative myofilament motion.

-     motions of the myofilaments should produce variations in the axial stress, in accordance with the well known force-velocity characteristic of muscle. But they may also produce non-axial shear stresses on planes parallel to the fiber axis; the latter would be absent in the case of purely axial deformation. Do non-axial deformations actually affect muscle stress? More specifically, do variations of the lateral spacing of the myofilaments influence cross-bridge dynamics? Huxley (1980) thought that they do not, basing his opinion primarily on two phenomena on the descending limb of the length-tension curve: (i) a strict linear variation of isometric force with sarcomere length, and (ii) a lack of variation of the maximum shortening speed. Edman & Hwang (1977) reached the same conclusion, based on their measurements of a constant maximum shortening speed over a limited range of sarcomere lengths (Edman, 1979) in intact frog muscle fibers. These authors believed that the changes in contractile behavior which accompany osmotic perturbations of the fiber cross-section were more likely attributable to changes in ionic concentration than changes in lateral myofilament spacing. But several more recent studies on skinned fibers, where ionic concentrations are controlled, suggest a direct effect of lateral spacing on cross-bridge interactions. Metzger & Moss (1987) found a weak effect on force and a strong effect on maximum shortening speed of osmotic compression in chemically skinned rat fibers, and interpreted these effects as due to lateral-spacing dependence of the actin-myosin interaction rates. Goldman (1987) measured similar osmotic influences on shortening speed in mechanically skinned frog fibers, but interpreted them as the effects of filament spacing on effective cross-bridge stiffness. Krasner & Maugham (1984) observed a large decrease of both isometric force and ATP hydrolysis with osmotic compression in chemically skinned rabbit psoas fibers, and attributed these observations to hindrance of cross-bridge movement with decreasing filament spacing. Goldman & Simmons (1986) believed that a change of myofilament separation accounted for the effects of osmotic perturbations on stiffness that they observed in mechanically skinned frog fibers. Thus, it appears at present that a direct effect of filament spacing on cross-bridge dynamics is experimentally plausible. Indeed, Bagni et al. (1990) have also drawn this conclusion from their experiments on intact fibers. Further, Edman & Reggiani (1984, 1987) have published observations that call into question two cornerstones of the argument against lateral-spacing effects: (i) the existence of a ‘‘plateau’’ in the

61

length-tension curve, and (ii) the linearity of the descending limb; this will be discussed in more detail in the Results section to follow. It should be emphasized that osmotic perturbations are the only type of experimental non-axial deformations that have been published so far, and these are somewhat limited in that the fiber cross-section deforms isotropically—that is, the strain perpendicular to the fiber axis is independent of direction. To the best of my knowledge there have been no experiments reported on the effects of non-axial deformation rates, comparable to the large number of axial force-velocity studies. Given the likelihood that non-axial strains and strain rates affect the active muscle stresses, and the relative paucity of experimental data addressing these effects, it seems desirable to have an appropriate quantitative theory that would permit at least order-of-magnitude estimates. The object of this paper is to develop such a theory. I will present a generalization of the basic Huxley (1957) two-state model which allows for three-dimensional motions of the cross-bridges, abandoning the traditional restriction to purely axial motions. Although many refinements of the Huxley theory have appeared since its introduction (including larger numbers of states, ‘‘self-consistency’’ constraints, competition among myosin heads, multiple actin sites, and nonlinear cross-bridge elasticity) this appears to be the first attempt to account for non-axial cross-bridge strains and strainrates. The initial restriction to a two-state model is reasonable, as one should solve this simpler problem (which, as the sequel will show, is by no means trivial) before proceeding to more complex multi-state cases. Further, as Brenner (1990) has argued, two-state models can still be quite useful as a first-order approximation of contraction dynamics if one interprets the ‘‘attached’’ and ‘‘detached’’ states to represent, respectively, the sets of strong-binding and weak-binding states. The goals of this paper are to formulate the theory mathematically, apply it to a few simple illustrative problems, and draw some general inferences from its structure. More detailed calculations and predictions must be reserved for future publications. The task of formulating the theory can be divided into the following three steps: (1) Establish a relation between the macroscopic deformation of the muscle and the relative motion between a cross-bridge and its actin binding site. In Huxley’s 1957 theory this is accomplished by assuming that the relative sliding velocity between thick and thin filament,

. . 

62

v, is equal to (sV/2L) where s is the sarcomere length, V the muscle’s velocity, and L its length. (2) Establish an equation for the evolution of the actin-myosin bond-distribution function, n. For this purpose Huxley developed his classic rate equation 1n 1n − v(t) = f(x)(1 − n) − g(x)n 1t 1x

(1)

where f(x) and g(x) are, respectively, the bonding and unbonding rate functions. (3) Establish relations between the bond-distribution function and measurable macroscopic variables such as stiffness, K and stress, s. The Huxley 1957 theory may be shown to lead to K=

and

s=

0 1g ms 2k 4l

n(x, t)dx

−a

0 1g msk 2l

a

a

xn(x, t)dx,

(2)

−a

where the factors in parentheses are combinations of microstructural parameters, defined in Huxley’s paper. I will follow these three steps in the present development although, naturally, the details will be more complicated for three-dimensional cross-bridge motions. Due to space restrictions certain calculations will not be exhibited explicitly, but these details are available in a technical report (Zahalak, 1995). The Relation Between Macroscopic Deformation and Cross-bridge Strain The following arguments use certain standard concepts from the theory of continuous media as explained, for example, in Eringen, 1962. Figure 2(a) shows a muscle in a reference state, where it is by definition ‘‘undeformed’’. Figure 2(b) shows the same muscle at some later time as it is contracting—or, more generally, deforming. If the muscle is viewed as a continuous medium, and u represents the absolute velocity of its particles, with respect to a fixed laboratory frame, as a function of position and time, it is well known from continuum theory that the relative velocity between two neighboring points separated by a differential position vector dx can be calculated as du = d · dx + BV · dx,

(3)

where d = (u9 + 9u)/2 is the rate-of-deformation tensor, and V = (u9 − 9u)/2 is the material’s

rate-of-rotation tensor; 9 is the gradient operator. This is the relative velocity in a fixed frame. Suppose, however, that the vector dx emanates from the origin of a local moving frame of reference whose rate of rotation in the fixed laboratory frame is given by the antisymmetric rotation tensor V . Then the relative velocity of two neighboring points in the moving frame will be dv = D · dx, where D = v9 = d + V − V

(4)

is the (conjugate) gradient of the muscle velocity field with respect to the moving frame. The moving reference frame (x1 , x2 , x3 ) will be defined so that its origin is attached to and moving with a material point of the muscle, the x3 direction coincides at all times with the local direction of the myofilaments, and the x1 and x2 directions are mutually perpendicular and lie in a plane perpendicular to x3 . The precise orientations of x1 and x2 will be specified presently. Thus, according to a continuum model the relative velocity of two neighboring points can be expressed in terms of the relative velocity gradient D as dv = Dij dxj ei , where ei is the orthonormal base vector associated with the xi coordinate system. (The summation convention for repeated indices is used throughout this development.) In qualitative terms one can express this by saying that the spatial rates-of-change of the velocity evaluated at any point in the muscle define completely the relative motions of all material particles near that point. It will now be assumed that this relative continuum velocity manifests itself as a relative velocity of cross-bridges and their actin binding sites, but this manifestation is constrained by the assumed rigidity of the myofilaments. The specific implications of this last statement are illustrated in Fig. 3(a), which shows a single thick filament surrounded by a cage of thin filaments. The four velocity-gradient components in a plane perpendicular to x3 , D11 , D22 , D12 , and D21 , are assumed to be transmitted to the myofilaments without modification, it being required only that dx1 and dx2 represent the relative position of the filaments. Similarly, the axial shear components, D31 , and D32 , and are assumed to translate directly into relative myofilament motion [Fig. 3(b)]. This is not the case, however, for the transverse shear components, D13 , and D23 . As illustrated in Fig. 3(c), to first order in the time increment these modes of deformation produce neither axial nor transverse relative displacements between points fixed, respectively, on the thick and thin filaments. Indeed, as x3 is constrained to remain along the fiber directions, these velocity-gradient components

-    

3

dX

3

X

1

X

dX (a)

63

Undeformed (reference) configuration

dR

dR

3

dx

x1

dx



z dX 3

x3

X2

Deformed (contracted) configuration

λxdR

λydR

F. 2. (a) The undeformed reference state of a muscle, showing the material coordinates, XI . (b) The deformed current state of the muscle, showing the spatial coordinates, xi .

represent rigid body rotations that produce no strain. Therefore, D33 , and D23 , will simply be dropped from the velocity gradient when computing relative cross-bridge displacements. Finally, the component D33 , produces a relative axial velocity D33 dx3 . If dx3 is set equal to a half-sarcomere length, s/2, this will give the relative velocity between the center of a thick

filament and the z- line. As in Huxley’s original theory, it will be assumed that this velocity, D33 (s/2), is equal to the axial rate of separation of all cross-bridges and their actin-sites on the ‘‘right’’ half of a sarcomere. This completes the translation of the macroscopic relative velocity field into cross-bridge motion, in terms of the relative velocity gradient v9.

. . 

64

dx1

D21dx1 + D22dx2

dx2 D11dx1 + D12dx2 (a)

D31dx1 dx1 (b)

D13dx3

(c) dx3

F. 3. Relations between macroscopic deformation rates and cross-bridge deformation rates (a) Transverse deformations. (b) Shear parallel to the fibers. (c) Shear perpendicular to the fibers.

The Rate Equation for the Bond-Distribution Function The next step is to develop a mathematical representation of the cross-bridge bonding kinetics. For this purpose I adopt a simplifying assumption

from Huxley’s original model: at any time instant each cross-bridge has a significant probability of interacting with only one actin binding site, which can conveniently be labeled the ‘‘nearest’’. Consider an ensemble of cross-bridges, and assume that for accounting purposes each cross-bridge has been

-     transported without change of orientation to the origin of a (x, y, z) Cartesian coordinate system. For each cross-bridge there will be, by assumption, one actin binding site to which the cross-bridge may or may not be attached. Thus, the accessible binding sites form a moving ‘‘swarm’’ about the cross-bridges at the origin, and this construct will be referred to as cross-bridge space, and is illustrated in Figure 4(a). (The ensemble of superimposed cross-bridges at the origin, each with its correct orientation, will be called the representative cross-bridge.) The (x, y, z) axes of cross-bridge space are aligned with the (x1 , x2 , x3 ) axes of the fiber coordinate system. It will be convenient to differentiate between position in cross-bridge space, denoted by x = (xex + yey ) + zez = rer + zez

(5)

Theorem of continuum mechanics (Eringen, 1962) to write

g$ dV

%

DN + N9 · v dV = Dt

g

dV

N dV =

g

[F(N − N) − GN] dV

dV

where DN/Dt = 1N/1t + v · 9N is the material derivative of N. The Transport Theorem effectively accounts for the cross-bridge fluxes mentioned in the preceding parenthetical comment. Finally, as eqn (8) must hold for an arbitrary volume, dV, we may conclude that the integrands on both sides of this equation are equal, so that 1N 1N + N9 · v + v · 9N = + 9 · (Nv) 1t 1t = F(N − N) − GN.

(6)

although ex = e1 , ey = e2 , and ez = e3 . The kinetic rate equation describing cross-bridge attachment and detachment is formulated most easily for a group of actin sites of fixed identity—that is, for a group of sites ‘‘moving with the swarm’’. Thus, we can conceive of the actin sites as particles in a continuum (like a fluid) moving with a velocity v. Let N (x, y, z, t) dV be the total number of accessible actin sites in a volume element of cross-bridge space dV at time t, and N(x, y, z, t) dV be the number of those sites that are bonded to myosin. Let F(x, y, z,) and G(x, y, z,) represent, respectively, the attachment and detachment rates, considered to be functions of position in cross-bridge space—that is, relative position of cross-bridges and actin sites. Then if dV represents a volume in cross-bridge space bounded by a surface dS that moves with the swarm velocity v, no actin sites cross dS and one may write d dt

g

(8)

and position within the muscle, denoted by x = x1 e1 + x2 e2 + x3 e3

65

This is the general form of the rate equation governing the evolution of the bond distribution function, N. Next, the specific form of the relative velocity field developed in the last section will be introduced into eqn (9). Expanding into components the relation between the velocity and the velocity gradient, one may write v = Dxx xex + Dxy yex + Dyx xey + Dyy yey s +Dzx xez + Dzy yez + Dzz ez 2 s = D · r + e z · D · ez e z 2

(10)

The velocity components proportional to Dxz and Dyz are absent for reasons explained previously. Instead of Cartesian (x, y, z) components the velocity gradient can be expressed in terms of polar (r, u, z) components, in which case

0

1

s v = rDr er + rDur eu + rDzr + Dzz ez 2 [F(N − N) − GN] dV

(7)

dV

This simply states that the rate of increase of attached cross-bridges within dV is the difference between the attachment and detachment rates. (If dV were a fixed volume in cross-bridge space one would have to consider the fluxes of attached and detached cross-bridges into and out of dV—clearly a more complex calculation.) As dV is a region moving with the continuum we may invoke Reynolds’s Transport

(9)

(11)

where er is the radial unit vector, eu is the circumferential unit vector, and Drr = Dxx cos2 u+Dyy sin2 u+sin u cos u(Dxy+Dyx ) Dur = (Dyy − Dxx )sin u cos u + Dyx cos2 u −Dxy sin2 u Dzr = Dzx cos u + Dzy sin u Dzz = Dzz .

(12)

. . 

66

The relations, eqn (12), between the Cartesian and polar components of D are standard properties of a second rank tensor and are discussed in any comprehensive book on continuum mechanics. Note that the components of D are in general functions of position within the muscle and time. In order to justify subsequent averaging operations it must be assumed that in each elementary volume of muscle the Cartesian components of the velocity gradient, Dxx , Dxy , . . . , Dzz are constants at any time instant; they are not functions of position in cross-bridge space. By eqns (12) this means that the polar components of the velocity gradient, Drr , Dur , Dzr , and Dzz , are functions of time and polar angle at each point in the muscle. Inserting eqn (11) into the general rate equation, eqn (9), one obtains

0

1

1N 1N 1 1N s 1N + rDrr + rDur + rDzr + Dzz 1t 1r r 1u 2 1z

0

+ 2Drr +

1

1Dur N = F(N − N) − GN. 1u

(13)

At this point it must be recognized that the actin sites accessible to the representative cross-bridge are not distributed throughout cross-bridge space, but rather are confined to a two-dimensional surface—the cylinder generated by the thin filaments surrounding a thick filament. The cross-bridges on a particular thick filament have actin sites accessible only at six discrete angles. But if the orientations of actin lattices in the separate fibrils of a fiber are uncorrelated, as suggested by Fig. 4(b), then a representative cross-bridge picked at random from any fibril can have its binding site oriented at any angle in an instantaneously fixed frame, presumably with equal probability. This motivates the following assumption: in an undeformed reference configuration all the actin binding sites accessible to the representative crossbridge are uniformly distributed on a circular cylinder of radius a, where a is the length of a hexagonal side of the actin lattice. Further, the assumption that a cross-bridge interacts only with the ‘‘nearest’’ site implies that the axial extent of this cylinder is − (l/2) Q z Q (l/2), where l is the axial distance between successive binding sites, as shown in Fig. 4(c). As the muscle contracts and deforms this initially circular cylinder will deform into a cylinder of elliptic cross section, at all times parallel to the fiber axis. Let r = r0 (u, t) be the polar equation of its cross-section. As shown in Zahalak (1995) this equation can be calculated from the macroscopic deformation of the muscle tissue. The main results are summarized in

Appendix A, and are as follows. A small unit circle perpendicular to the fiber axis in the undeformed configuration is deformed into an ellipse, which is in general not perpendicular to the fiber axis in the deformed configuration. But the projection of this deformed ellipse on a plane perpendicular to the fiber axis is also an ellipse. Let the x1 and x2 (and therefore the x and y) coordinate directions coincide at each point in the muscle with the principal axes of the projected ellipse. Further, let lx and ly denote the lengths of the principal axes of the projected ellipse; as the undeformed circle was assumed to be of unit radius these represent principal stretch ratios of the projection of the macroscopic deformation on a plane perpendicular to the local fiber direction (see Fig. 2). This particular choice of coordinates permits the theory to be expressed in its simplest mathematical form. It follows that the cross-section of the deforming actin cylinder can be described by x = alx cos v and y = aly sin v—the standard parametric equation of an ellipse—and therefore r0 (u, t) = azlx2 cos2 v + ly2 sin2 v ,

(14)

where tan u = (ly /lx )tan v. Note that lx and ly are in general functions of time at each point in the muscle. It is also shown in Appendix A that if two neighboring generators of a circular cylinder in the undeformed configuration are separated by a distance dS then in the deformed configuration they will be separated by a distance ds given by ds = zlx2 sin2 v + ly2 cos2 v dS

(15)

Let m* represent the areal density of accessible actin sites, one for each cross-bridge on the undeformed actin cylinder in cross-bridge space. This density can be related to the volumetric density of cross-bridges, m, as follows: if the ensemble is taken as the set of cross-bridges/actin sites in one half-sarcomere then bpa 2(s0 /2)m = 2palm*, or m* = (bas0 /4l)m, where s0 is the sarcomere length in the undeformed reference state and b = 0.827 . . . is a numerical constant equal to the ratio of the area of a regular hexagon to its circumscribed circle. This calculation assumes that all the cross-bridges are interacting with their actin binding sites. Then by eqn (15) the areal density of binding sites on the deformed actin cylinder is m* dS dz m* . = ds dz zlx2 sin2 v + ly2 cos2 v

(16)

Not all the cross-bridges may be participating in the contraction, due to reduced myofilament overlap at long sarcomere extensions or steric hindrance at short

-    

67 v

x

χ (a)

z

y Fibril 2 Fibril 1

(b)

Fibril 3

x,x1

a z, x3

(c) a

l

a

F. 4. (a) Cross-bridge space, showing the representative cross-bridge surrounded by a swarm of actin sites. (b) Orientations of actin cages in separate fibrils are assumed uncorrelated. (c) The undeformed actin cylinder, showing a distribution of actin sites about the representative cross-bridge.

lengths. In that case the constant m* in eqn (16) must be multiplied by a factor a equal to the fraction of participating cross-bridges; a is a function of fiber length. Thus, the areal density of accessible binding sites depends on polar angle u and time, but not on axial position, z. With this in mind, as accessible

cross-bridges are found only on the surface of the deforming actin cylinder, their distribution must have the particular form N (r, u, t) = m*an˜ (u, t )

d[r − r0 (u, t)] r0 (u, t)

(17)

. . 

68

where d is the Dirac delta function, and n˜ is a reduced distribution function of accessible actin sites, depending only on polar angle and time. As the number of accessible sites in a wedge du dz is equal to du dz f0a N r dr, and also to dsdz times the r.h.s. of eqn (16) a brief calculation shows that this requires n˜ = a[r0 (u, t)/a]2/lx ly . Similarly, as attached cross-bridges are confined to sites on the deforming actin cylinder, we can introduce a reduced bond-distribution function n(z, u, t), by the equation N(r, u, z, t) = m*an(z, u, t)

d(r − r0 ) . r0

(18)

To derive the rate equation governing n, simply multiply eqn (13) through by r and integrate with respect to r over the interval (0, a), remembering that the Dij are functions of u but not of r. Further, it is assumed at this point that the attachment and detachment rate functions, F and G, depend on r and z, but not on the angular coordinate u; this is a reflection of the assumed rotational symmetry in the orientation of the cross-bridges. The calculation is straightforward (Zahalak, 1995) and yields 1n 1n 1n s 1n 1D + Dur + Dzr r0 + Dzz + ur n 1t 1u 1z 2 1z 1u =f (n˜ − n) − gn

1n 1n 1n 1n s = − v = f(z) − [ f(z) + g(z)]n. + D 1z 1t 2h zz 1w 1t Equation (20) applies to the cross-bridges on one side of the M-lines—the ‘‘right’’ half-sarcomeres. A similar equation must be written for the ‘‘left’’ half-sarcomeres. For this purpose one of two coordinate systems can be established on the left side, as shown in Fig. 5(a); these are obtained by rotating the original fiber coordinate system on the right through 180° about either the x1 or x2 axes. These rotations produce two new coordinate systems in cross-bridge space, (x', y', z') and (x0, y0, z0), and we denote the reduced bond distribution function in the left half-sarcomere as either n'(x', y', z', t) or n0(x0, y0, z0, t). It is not necessary to actually write and solve a rate equation in the left halfsarcomere, as it can be shown (Zahalak, 1995) that, in addition to the periodicity condition n(w, u 2 2p, t) = n(w, u, t), the bond-distribution function satisfies the following symmetry conditions based on the form of eqn (20): n'(w', u', t) = n(w', p − u', t)

(19)

where f(z, u, t) = F[r0 (u, t), z] and g(z, u, t) = G[r0 (u, t), z]; note that the radial velocity gradient, Drr has disappeared in the integration. This equation may be put into a convenient dimensionless form if we define h as Huxley’s bond-length scaling parameter (a measure of the range of bond lengths over which there is a significant probability of attachment), w = (z/h) as the dimensionless axial distance, and the parameters m = (s0 /2h) and e = (a/ s0 ). Then eqn (19) becomes 1n 1 1n + (D n) + m{lz Dzz + 2erDzr } 1t 1u ur 1w =f (n˜ − n) − gn

If one assumes no non-axial deformation then lx = ly = r = n˜ = 1, Dur = Dzr = 0, and one recovers the classic Huxley rate equation

(20)

where r(u, t) = r0 (u, t)/a and lz is the axial stretch ratio, implying s = lz s0 . Recall that n˜ (u, t) = ar 2(u, t)/lx ly , which includes the effect of overlap in the factor a. Equation (20) is the main result of this analysis. It generalizes Huxley’s original rate equation so that the bond-distribution function n(w, u, t) at each point in the muscle tissue depends not just on the axial deformation and deformation rate, lz and Dzz , but also on the non-axial deformation through r, and the non-axial deformation rates through Dur and Dzr .

n0(w0, u0, t) = n(w0, −u0, t).

(21)

Stresses and Stiffnesses The third and last step of the theoretical program outlined in the Introduction is to develop relations between the bond-distribution function, n, and the instantaneous state of active stress and stiffness in the muscle. Brief reflection indicates that a precise calculation of stress and stiffness would be very difficult, and would require detailed statistical information about the geometrical distribution of fibrils and non-contractile material. As the present theory is based on a simplified two-state cross-bridge model, it is doubtful that such a high degree of rigor is warranted. Therefore, the stress and stiffness calculations will be based on the behavior of a single ‘‘representative’’ sarcomere, which is probably sufficiently accurate at the present level of modelling. Further, in deriving the relations to follow, it will be assumed that the hexagonal actin unit cell can always be approximated by its circumscribing cylinder, thus effectively setting b = 1; this may actually compensate to some extent for the fact that muscle contains some proportion of non-contractile material between the fibrils.

-     x'

69

x

y'

z (a)

z', z" y

y"

x" z

(b)

r

r*

λ xa

l F31

l F33

r F33

r F32

l F32

r F31

l F23

(c)

λ ya

r F21

l F22

l F21

r F23

r F22

r F11 l F11

l F13

l F12

r F13

r F12

λzso F. 5. (a) Coordinate systems for the right and left half-sarcomeres. (b) Cross-bridge and actin binding site, showing schematically axial and radial cross-bridge elasticity. (c) Various sections of the actin cylinder surrounding a representative myosin filament that are used in the computations of the stress components, showing the resultant forces acting on these sections.

To introduce forces into the model it is assumed, as shown in Fig. 5(b), that a cross-bridge is elastically restrained against both axial and radial motions, with spring constants kz and kr , respectively. Despite the

suggestive diagram, no specific assumptions have been made about the structure or kinematics of the cross-bridge, except that there is no significant cross-coupling between the axial and radial modes of

. . 

70

deformation (otherwise yet another parameter would enter the model). The radial cross-bridge position at which the radial force vanishes is denoted by r* and, for convenience, z is measured from the zero-force axial position, With this model the force of a single attached cross-bridge is fCB = kr (r − r*)er + kz zez .

(22)

The instantaneous state of stress at any point within the muscle can be described by the Cauchy stress tensor, s = sij ei ej , whose nine components represent the forces per unit deformed area acting on planes perpendicular to the coordinate axes. Consider first s33 , the direct axial stress [Fig. 5(c)]. The resultant axial force on the right myosin halffilament is

whereas the corresponding force (in the opposite direction) in the left half-sarcomere is 1 F13 = kz

0

ma 2s0 h 2 4l

1g 6g a

z'

−a

(27) The shear stress is the difference of these two forces assigned to the area of intersection of an axial plane with the deformed actin cylinder—that is r 1 − F13 )/2(ly a)(lz s0 ). Using the symmetry s13 = (F 13 relations, eqns (21), this can be reduced to s13 =

01 g 6 g a

1 p a s ly lz 2 s0 0

1 2p

z

−a

F = kz zN dV

= kz

0

1g

a

g

w dw

−a

2p

n(w, u, t) du.

(23)

0

This force is assigned to the deformed elliptic cross-section of the actin cylinder, p(lx a)(ly a). Therefore, s33 =

s0 lx ly

g 6 g a

1 2p

w

−a

7

2p

n(w, u, t)du dw,

0

F = kz

0

1g 6g a

−a

z

p/2

−p/2

7

n du dz.

p/2

(28)

r F 12 = kr

0

1g 6g

ma 2s0 h 4l

a

p/2

−a

−p/2

[r0 (u, t) − r*]

7

× sin un(z, u, t) du dz

(29)

(25)

is a normalizing constant with dimension of stress that depends on the microstructure and axial cross-bridge stiffness. It follows from Eqns (21) that 1 F33 the total axial force on the left myosin r half-filament is exactly the same as F 33 of eqn (23), so that the myosin filament is in axial equilibrium. Next the shear stress s13 will be calculated: this is the force per unit deformed area in the axial three-direction acting on a plane perpendicular to the transverse one-direction. Figure 5(c) illustrates the basis of this calculation. The resultant force exerted by all the cross-bridges crossing the x2 − x3 plane in the right half-sarcomere is ma 2s0h 2 4l

g

3p/2

whereas on the left

s0 = (ms0 kz h 2/2l)

r 13

n du

−p/2

For the shear stress acting in the x2 direction on the x2 − x3 plane, the algebraic sum of the appropriate forces acting on the two half-sarcomeres must again be determined [Fig. 5(c)]. On the right

(24)

where

p/2

1 2p



bma 2h 2s0 4l

n'(z', u', t) du' dz'.

−p/2

g

r 33

7

p/2

7

n(z, u, t) du dz (26)

1 F12 = kr

0

1g 6g

ma 2s0 h 4l

a

p/2

−a

−p/2

[r0 (u', t) − r*]

7

× sin u'n(z', u', t) du' dz'.

(30)

Again employing the symmetry relations, eqns (21), r 1 s12 = (F 12 − F12 )/2(ly a)(lz s0 ) becomes

0 10 10 1 g 6 g

1 p a ly lz 2 s0

s12=

× sin un du − where r* = r*/a.

a h

1 2p

kr s kz 0

g

3p/2

p/2

a

−a

1 2p

p/2

(r − r*)

−p/2

7

(r − r*)sin un du dz

(31)

-     In a similar fashion the remaining six components of the Cauchy stress tensor can be calculated. Due to space constraints these calculations can not be exhibited here, but the results are presented in Appendix B, and the details are available in (Zahalak, 1995). To summarize the stresses, and the stiffness to be developed, it is convenient to introduce the notation u1 , u2 =p, q =k, m, s =

g 6 g a

zs

−a

1 2p

(r − kr*)m

u1

7

(32)

With this notation we may write

01

− p/2, 3p/2=0, 0= 0, 0, 1}

1Dxi 1t

.

(34)

xj

But the relative velocity field has already been described in terms of the relative velocity gradient D [eqn (4)] so that

0 1

1 1Dxi 1xj 1t

ei ej .

(35)

xk

Dij =

010 1 1 1t

xk

1Dxi . 1xj t

(36)

Suppose at t = 0 a rapidly varying velocity gradient— in the limit a Dirac d-function of time—is imposed on the muscle element: Dij = gij d(t). Inserting this into eqn (36) and integrating with respect to time yields

0 10 10 1 kr kz

× {−p/2, p/2=0, 1=1, 1, 0 − p/2, 3p/2=0, 1=1, 1, 0}.

0 1

The order of space and time differentiation can be interchanged so that

1 p a {−p/2, p/2=0, 0=0, 0, 1 ly lz 2 s0

a h

vi =

D = Dij ei ej =

1 (s33 /s0 ) = 0, 2p =0, 0=0, 0, 1 lx ly

1 p a (s12 /s0 ) = ly lz 2 s0

posed on an element of muscle tissue one can proceed as follows. Starting at some arbitrary time instant, labelled t = 0 without loss of generality, let Dxi represent the displacements of material particles in the muscle relative to the origin of the moving fiber coordinate system. These displacements are functions of position xi at t = 0, and time, and the relative particle velocity is

u2

× cosp u sinq un du dz.

(s13 /s0 ) =

71

(33)

The remaining stress components are listed in Appendix B. Note that sij $ sji —that is, the active stress tensor is not symmetric. Finally, the instantaneous stiffness associated with cross-bridge interactions will be calculated. This stiffness is defined as the ratio between a change in stress and a rapid—in theory, instantaneous—change in strain. As both the stress and strain are second-rank tensors, each with nine components, the stiffness is a fourth-rank tensor with a total of eighty-one components. For example, in general a quick axial stretch gives rise not only to a change in axial stress (associated with an axial stiffness component) but also to a change in shear stress, producing a shear stiffness component. Fortunately not all these components are independent, and many turn out to be zero, as shown below. To describe mathematically an instantaneous deformation im-

gij =

0 1

1Dxi . 1xj t

(37)

So gij , which is not symmetric, is a relative material position gradient, based on the positions of particles in the fiber coordinate system at an arbitrary time instant. We may write Dxi = gij dxj (noting that Dxi vanishes at the origin xj = 0), or Dx = g · dx.

(38)

Thus, the relative displacements in a continuum are related to the incremental position vector dx in exactly the same way as the relative velocities, with D replaced by g. But it has been argued previously that not all the components of the velocity gradient produce relative cross-bridge displacements, and the same must be true of the displacement gradient. Therefore, based on eqn (10) we may express the

. . 

72

incremental changes of the position vector in cross-bridge space, Dx as s Dx = g · r + gzz ez 2

DT33 =

0

1

s = rgrr er + rgur eu + rgzr + gzz ez 2

(39)

0

1

s DF = kr rgrr er + kz rgzr + gzz ez 2

DF = kz

0

1g 6g 0 a

−a

2p

0

1 7

and on the left half-filament is

0

ms0 a 2h 4l

1g 6g 0 a

−a

2p

r0 g'zr

0

1 7

s + g'zz n' du' dz'. 2

(42)

But gij is a second -rank tensor obeying the standard transformation rules (Zahalak, 1995) g'zr (u') = −gzx cos u' + gzy sin u' gzr (u) = gzx cos u + gzy sin u

a

−a

g'zz (u') = gzz

gzz (u) = gzz .

(43)

s0 2a

1 2p

a

−a

1 2p

1 2p

2p

n du dz

0

2p

r cos un du dz

0

g 6 g a

−a

(41)

1 DF33 = kz

+ g31

= +g32

s r0 gzr + gzz n du dz 2

1

ms0 kz h 2 1 lx ly 2l

$ 0 1g 6 g 7 g 6 g 7

(40)

From this point on the calculation of the stiffness coefficients is mechanical, albeit lengthy. One rather subtle point should be made, however. If a muscle element is rapidly deformed then both the force on a selected area element and the area of the element itself will change. The change in the Cauchy stress sij (force per unit current area) is affected by both of these factors. It is easier to express the stiffness in terms of a Lagrange stress, Tij , defined as the force after the deformation is imposed per unit area before the deformation is imposed. Thus, DTij is affected only by the change in force, and not the change in area. Consider a force change in the axial direction. The resultant force change on the right myosin halffilament is ms0 a 2h 4l

0 10

r 1 − DF33 a DF33 = h p(lx a)(ly a)

× g33 lz

which, as the circumferential stiffness has been assumed negligible, produces an instantaneous change of cross-bridge force equal to

r 33

These relations plus the symmetry conditions, eqns (21), lead to

2p

7 %

r sin un du dz

0

=DT0 (B3333 g33 + B3331 g31 + B3332 g32 ).

(44)

where DT0 = (a/h)s0 is a normalizing constant with dimensions of stress, and, using a previously introduced notation B3333 = (lz /lx ly )(s0 /2a)0, 2p =0, 0=0, 0, 0 B3331 = (lx ly )−10, 2p =1, 0=0, 1, 0 B3332 = (lx ly )−10, 2p =0, 1=0, 1, 0.

(45)

The form of this last result is quite general. Thus, for any ij and pq we can write (DTij /DT0 ) = Bijpq gpq .

(46)

As noted previously the fourth-rank active stiffness tensor, Bijpq , has 81 components, but it can be shown that there are only 21 non-vanishing independent components. These are listed in Appendix C, and the details of the calculations are provided in Zahalak, 1995. Note that this active stiffness tensor does not satisfy the usual symmetries of linear elasticity theory: Bijpq = Bjipq , Bijpq = Bijqp , and Bijpq = Bpqij.

Applications The main objective of this paper is to present the development of the theory for non-axial effects on muscle contraction that is embodied in eqn (20) and Appendices B and C. Solutions of the rate equation for various specific experimental conditions are not trivial, and must be reserved for future communications. In order to illustrate the physical content of

-     the theory, however, three simple applications will be discussed.

A considerable literature exists concerning the effects of osmotic perturbations on contractile behavior in both intact and skinned fibers. For present purposes I will focus on the experiments of Brenner & Yu (1991), and interpret the results in terms of the present theory. When subjected to osmotically active solutes that are excluded from the myofilament lattice the cross-section of skinned fibers deforms isotropically without appreciable change of sarcomere length; thus the sarcomere volume changes. Under these conditions eqn (20) has the steady-state isometric solution r2 f[r(u), z] a 2 f[r(u), z] + g[r(u), z] lr

(47)

where lr = lx = ly is the radial stretch ratio and a is the fraction of participating cross-bridges, which may be assumed equal to one at full overlap and maximal activation. To proceed further some definite assumption must be made about the rate functions, and for simplicity it will be assumed— following Huxley’s original formulation—that f = 0 outside the interval 0 Q z Q 1, whereas f(r, z) = f1 (r)z and g(r, z) = g1 (r)z within this interval. Further, assuming that r0 = a in the undeformed reference state, lr = r0 /a = r, so that eqn (47) simplifies to n = [1 + g1 (r)/f1 (r)]−1 = F(r)

(48)

which is a pure function of the radial strain. From Appendix B one finds that the only non-zero stress components are szz = s33 =

1 sQ lr2 0 1

srr = s11 = s22 =

0 10 10 1

1 a lr lz s0

a h

the stresses by Fz = szz pa 2lr2 and Fr = 2pas0 lr lz srr so according to eqns (49) Fz = pa 2s0 Q1

      

n=

73

kr s Q (r − r*) kz 0 0 (49)

where Q0 = f n dz = F(r) and Q1 = f zn dz = F(r)/2 are the first two moments of n. Brenner & Yu worked with chemically skinned rabbit psoas fibers near 0°C, and reported their results in terms of Fz and Fr , respectively the axial and radial ‘‘forces per single thick filament.’’ These are related to

Fr = 2pa 2lr2

0 10 1 a h

kr s Q (r − r*). kz 0 0

(50)

Brenner & Yu (1991, fig. 2) found that Fz was almost constant with compression, varying by about 25% about a mean value. This suggests, by the first of eqns (50), that in this preparation the rate function ratio g1 /f1 was almost constant and insensitive to lateral filament motion, over the experimental range. A numerical value can be assigned to this ratio using Brenner’s (1988) measurements of the apparent attachment and detachment rate constants in rabbit psoas under isometric conditions: g1 /f1 = gapp / fapp = 1s−1/4s−1 = 1/4. Muscle fibers expand laterally when they are skinned and then contract laterally when they are activated (Krasner & Maugham, 1984; Matsubara et al., 1984; Goldman & Simmons, 1986; Yu & Podolsky, 1990). For purposes of this calculation it will be assumed that the expansion and contraction are each 15% so that for a fully activated skinned fiber, unperturbed osmotically, r0 = a and r = lr = 1. Brenner & Yu found that under these conditions the myosin lattice spacing was d10 = 38 nm, and therefore a = (2/3)d10 = 25.3 nm. The sarcomere length is taken as s0 = 2.35 mm. Further, these authors found that the axial isometric stress was approximately szz = 1 kg cm−2c1 M dyne cm−2 [which is close to 1.24 M dyne cm−2 given by Kawai & Schulman (1985)]. With these numbers the first of eqns (49) gives s0 = 0.25 pN(nm)−2 = (ms0 kz h 2 )/2l. An estimate of the axial stiffness, kz , can be calculated from the experiments of Finer et al., (1994), which are basically corroborated by Ishijima et al. (1994), that under isometric conditions the average values of cross-bridge force and movement are respectively 3.5 pN and 11 nm, giving an estimated kz = 0.32 pN(nm)−1 . If the distance between successive actin binding sites is set equal to a half-turn of the actin helix, l = 35.8 nm, then the above value s0 of results in an estimate of the bond-length scaling parameter h = 18 nm

(51)

assuming that the cross-bridge concentration is 1.2 × 10−4 M (Tregear & Squire, 1973). This value is in the range 12–26 nm suggested by quick release experiments (Ford et al., 1981). The second of eqns (50) states that Fr varies linearly with lateral compression, which is essentially the

. . 

74 (P/P1)

1.0

1.66

λz F. 6. The plateau and descending limb portions of the isometric tetanic length-tension curve of frog tibialis anterior. Dashed lines: classical piecewise linear variation. Filled circles: data of Edman & Reggiani (1987). Solid lines: prediction of the theory presented in this paper, with the parameter k set equal to 0.8, based on osmotic compression of intact fibers (Blinks, 1964; Edman & Hwang, 1977; Mansson, 1993, 1994).

behavior found by Brenner & Yu (1991, Fig. 6). The radial force was found to vanish at d10 = d* 10 = 34 nm, giving r* = 0.9. Further, as r − r* = (d10 − d* 10 )/d* 10 , the second of eqns (50) gives

0 10 1

Fr 2pa 2 a = d10 − d* d* h 10 10

kr sQ, kz 0 0

(52)

and this ratio was experimentally found to be approximately 100 pN(nm)−1 . Using the previously calculated values of a, h, s0 , and Q0 , this gives kr /kz = 3.0

(53)

confirming Brenner & Yu’s conclusion that the radial and axial stiffnesses are of comparable magnitude. (Note that in a two-state model kz should probably be taken as closer to the T2 stiffness rather than the higher T1 stiffness.) Finally, returning to eqns (50) with all the established parameters we can calculate at maximum activation in the absence of osmotic perturbations, with r = 1, that Fr = 339 pN and Fz = 201 pN. Assuming 300 cross-bridges per thick filament this gives 1.13 pN and 1.34 pN, respectively, as the radial and axial forces per single cross-bridge.

Taking note of the Finer et al. (1994) value of 3.5 pN per cross-bridge, this suggests that only 1.34/ 3.5 = 38% of the cross-bridges are attached at any instant during an isometric contraction, a number which is not inconsistent with biochemical estimates (Squire, 1986) and nano-manipulation of single actin filaments (Ishijima et al., 1994).  -     Edman & Reggiani (1984, 1987) have published experimental results that contradict certain commonly accepted ideas about the tetanic isometric length-tension curve of intact frog skeletal muscle fibers. They found that (i) the sarcomere length-tension curve has no ‘‘plateau’’, or constant portion, and (ii) the descending limb is not linear, but rather concave-up. To observe these features Edman & Reggiani used a preparation with shorter marked segments of a fiber than had been used in previously published studies, and this enabled them to obtain ‘‘creep-free’’ tension-time records not requiring extrapolations or corrections. The theory of the present paper can predict these two features of the

-     length-tension relation from the response of intact fibers to osmotic perturbations, as the following arguments will demonstrate. Edman & Hwang (1977) found that hypo- and hyper-osmotic solutions of 0.81R and 1.22R (where 1.0R represents standard isotonic Ringer’s solution) produced, respectively, isometric force changes of +10 and −12% and area changes of +11 and −13%, in intact frog semitendinosus fibers near 0°C. They stated that doubling the osmotic concentrations doubled the force and area changes. These observations appear to be consistent with those of Mansson (1993, 1994) on intact frog tibialis anterior fibers, although area change are not explicitly reported in these papers. According to the first eqns (49) szz = l−2 r s0 F(lr )/2, as r = lr in this case. Then the muscle force is P = Al−2 r s0 F(lr )/2, where A is the deformed area, and the force at any degree of osmotic compression or expansion relative to P1 , the force in the reference state 1.0R, is

01

P Al−2 s F(lr )/2 F(lr ) = , = r 0 P1 A0 s0 F(1)/2 F(1)

(54)

(assuming that lr = 1 in the reference state) because the deformed and undeformed cross-sectional areas are related by A = lr2 A0 . But the l.h.s. of eqn (54) is measured in osmotic perturbation experiments and, as noted above, found to be a linear function of the cross-sectional area. Thus, we may write DA F(lr ) =1+k , F(1) A0

(55)

where k is an experimentally measured proportionality constant with a magnitude of the order of unity. The previous calculations assumed full overlap, with a = 1. If we now consider length-tension measurement on intact fibers a will be one only for sarcomere lengths less than a critical value (taken here as 2.25 mm) which defines the axial reference length, with lz = 1. For longer sarcomeres a will decrease linearly due to reduced filament overlap, according to a(lz ) = (l0 − lz )/(l0 − 1),

(56)

where l0 is the axial stretch ratio at zero overlap, taken here to be l0 = 1.66, implying a sarcomere length of 3.65 mm. Further, as deformations of intact fibers are volume conserving, lz lr2 = 1, which implies that (DA/A0 ) = (A/A0 ) − 1 = lr2 − 1 =

75

l−1 z − 1, which can be inserted into eqn (55). With a(lz ) in general different from one for intact fibers

01

P a(lz ) F(lr ) = P1 a(1) F(1)

8

1 + k(l−1 z − 1)

=

$

l − lz [1 + k(l − 1)] 0 l0 − 1 −1 z

%

for lz Q 1 for lz q 1

9

.

(57) This is the normalized isometric length-tension curve (excluding the ascending limb) that the present theory predicts for intact fibers. It is completely determined, if the myofilament lengths are known, by a single constant, k, which can be measured in osmotic compression experiments on skinned fibers. Equation (57) is plotted in Fig. 6 together with Edman & Reggiani’s data for frog tibialis anterior. [A value of k = 0.8 has been used in this figure, based on estimates from the Mansson papers cited previously; this is slightly lower than a value of 0.9 which appears appropriate for the semitendinosus of Edman & Hwang (1977).] While the correspondence with experiment is not perfect, the theory clearly exhibits the two features mentioned previously: rising force on the putative ‘‘plateau’’, and a concave-up descending limb. Edman & Reggiani (1987) have proposed a model that attributes their length-tension observations to assumed statistical variations of sarcomere length over a fiber cross-section. Edman & Hwang (1977) had previously attributed the osmotic perturbations of isometric force in intact fibers to changes in ionic concentration. The present theory attributes both sets of observations to a single mechanism— variation of cross-bridge attachment and detachment rates with filament spacing—which seems more parsimonious, although it certainly does not rule out some effects of sarcomere dispersion and ionic concentration.       The two preceding applications have been concerned with isometric steady-state conditions and therefore do not exhibit the effects of deformation rates. Although at present there are no published data on non-axial deformation rates, I will provide the solution to one simple but instructive problem. Suppose a block of contractile tissue is sheared at a constant rate, g˙ , with the shear planes parallel to the fiber axes, as shown in Fig. 7(a). The cross-section of the actin cylinder does not deform in this situation,

. . 

76 (a)

. γ

X1,x1

coefficients shown in Fig. 8 must be multiplied by appropriate scaling factors, listed in Appendix D, to calculate the actual stresses. The scaling factors for s11 , s22 , and s13 contain the ratio (a/s0 ) = O(10−2 ), implying that these stresses are of the order of one hundredth of the axial stress [in spite of the fact that, as shown in Application (1) above, the radial and axial forces per cross-bridge may be of comparable magnitude.] The stress s31 does not have this factor, but Fig. 8 shows that nevertheless this stress is small compared with the axial stress in this particular case. The axial stiffness coefficients (i.e. those that

X3, x3

(b)

σ11 X1,x1

Q0 n0

1

σ13

0

0.8 0.6 0.4 0

–1

σ31

(a)

σ33

–2 θ 0.5 1 °δ

(c)

X3,x3 1.5 2

–3

σ22 X2,x2

0 –1 –2 θ

0.5 1

δ°

1.5 2

–3

F. 7. (a) Simple shearing of a block of active contractile tissue parallel to the fibers. (b) Variation of the steady-state axial zero-order moment, Qo , with polar angle, u, and dimensionless shear rate, d . (c) Variation of the steady-state axial first-order moment, Q1 . Both Qo and Q1 are even functions of u. (see Appendix D.)

and therefore there are no radial variations of the rate functions f and g. The solution—which can be worked out, to a large extent, in closed form—is described in Appendix D. Figure 8(b) and (c) shows the variation of the stresses with dimensionless shear rate, d = (2a/h)(g˙ /( f1 + g1 )), for two sets of parameters characterizing the attachment and detachment rate functions: the first set (b) is appropriate for frog sartorius at 0°C, and the second set (c) is for cat soleus at 37°C, a muscle that shows a very pronounced tendency to ‘‘yield’’ on axial stretch. The shearing deformation decreases the three direct stresses (s11 , s22 , s33 ) and gives rise to two new active shear stresses (s13 and s31 ). Note that the stress

0.6

Stress coefficients

0.8 0.6 0.4 0.2 0 0

0.5

S33

S13

0.4 0.3

S22

0.2

S11

0.1

S31

0 –0.1 (c)

0

2

4

6

8

10

0.7 0.6

Stress coefficients

Ql n0

(b)

0.5 0.4

S33

0.3

S13 S22 S11

0.2 0.1

S31

0 –0.1

0

2

4

6 . Shear rate δ

8

10

F. 8. (a) Steady-state stresses produced by simple shearing. Variations of the stresses with dimensionless shear rate for the following values of the rate parameters (in inverse seconds: (b) f1 = 43.3, g1 = 10, g2 = 209, g3 = 19.6, and (c) f1 = 50, g1 = 14, g2 = 600, g3 = 41.7.

-     involve axial forces) are displayed in Fig. 9(b) and (c), and their physical interpretation is illustrated in Fig. 9(a). Again, the deformation rate decreases the

(a)

77

stiffnesses that existed under isometric conditions (K3333 , K2332 , K1331 ) and gives rise to new stiffnesses (K1333 and K3331 ). An examination of the scaling factors

K1333γ33

K1331γ31

K3333γ33

(b)

K3331γ31

1

0.8

Stiffness coefficients

K3333 0.6

0.4 K2332

K1331 0.2 K1333

K3331

0.0 0 (c)

2

4

6

8

10

1

Stiffness coefficients

0.8

0.6

K3333

0.4 K2332

K1331

0.2

K1333 K3331

0.0 0

2

4

. Shear rate δ

6

8

10

F. 9. (a) Illustration of stress and strain increments used for calculation of axial stiffness components during steady-state simple shearing. (b) and (c): variations of selected steady-state stiffnesses with dimensionless shear rate during simple shearing, for the rate parameter values listed in the caption of Fig. 8.

. . 

78

in Appendix D shows that all the stiffnesses are small compared to the axial stiffness, K3333 due to the presence of ratio (a/s0 ). It is also worth remarking that for the values of rate parameters chosen for this example the shear rate would have to be quite large to produce significant effects on the axial stress and stiffness. The ratio (2a/h) is approximately 3 whereas f1 + g1 is about 50 s−1 , so d = 1 implies that g˙ = 17 s−1—a high deformation rate. Clearly, the effect of deformation rate would be more pronounced for slower muscles. Thus, for the rabbit psoas at 0°C considered in Application (1), where f1 + g1 = 5 s−1 , d = 1 implies g˙ = 1.7 s−1—a much more modest deformation rate. Discussion To the best of my knowledge this is the first attempt to present a general theory for estimating the effects of non-axial deformation and deformations rates on muscle stress and stiffness. There appears to be no reason why the ideas used in the present development could not be applied to more accurate multistate models, but the moderate complexity of the present results suggest that a fuller exploration of the two-state 3-D Huxley model would be worth while. It may be that theoretical solutions for modes of deformation so far not much employed experimentally—such as shearing, torsion, or compression between rigid surfaces—could yield additional information on the molecular mechanisms of contraction. Certainly the present theory can serve as a basis for modelling muscular function in muscles with complex architecture, like the heart. Even without specific solutions, the structure of the model suggests some general characteristics of non-axial deformations. All variables and parameters in eqn (20) are dimensionless, except time. The stretch ratio lz and the function r are both 0(1), whereas the parameter e = 0(10−2 ); this suggest that the effect of the shear rate Dzr is small compared with the axial elongation rate, Dzz . Further, m = 0(102 ), which again suggests that Dzz dominates the transverse shear rate, Dur . Thus, it appears that these non-axial deformation rates will significantly influence the bond-distribution function only if the axial deformation rate is low, or if they are themselves very large—for example, in high-frequency vibrations. Turning to the expressions for stresses listed in Appendix B, it can be seen that all stresses except s33 , s31 , and s32 are proportional to (a/s0 ) and therefore expected to be small compared with the axial stress, s33 . While this general constraint does not apply to s31 , and s32 , it can be seen from the simple shear example of Appendix D that these

stresses may also be small compared with the axial stress. This argument holds with even more force for the stiffness coefficients listed in Appendix C, which shows that the ratio of every stiffness coefficient to the axial one, B3333 , is proportional either to (a/so ) or its square. Thus, the non-axial stiffnesses can be expected to be small compared with B3333 . In certain situations, such as unloaded shortening, the axial stress may vanish and non-axial stresses become relatively important. This remains for future studies to investigate. Although the preceding considerations suggest that, except perhaps for slow muscles, non-axial deformation rates will not have a strong effect on the axial stress and stiffness, this is not to say that the non-axial deformations themselves may not have such an effect. An examination of the literature on osmotic perturbations in muscles suggests that the axial force drops to zero when the cross-section is expanded to lr = 1.5 or is compressed to lr = 0.7. The important point is that the axial stress decreases markedly both as the filament spacing increases and decreases. This, in turn, suggests that any substantial deformation of the cross-section, isotropic or not, will decrease the axial stress. [This need not be true for sufficiently small transverse deformation, as in Application (2).] An analysis of this question with the present model will require a specification of the radial variations of f and g for moderately large transverse strains. This theory, in agreement with Brenner & Yu (1991), predicts comparable axial and radial crossbridge forces, but it predicts very different axial and non-axial stresses, and it is the stresses that determine the macroscopic mechanical behavior of muscle tissue. But the theory predicts only the active stress, s(a) , associated with cross-bridge interaction. The full stress tensor, s, acting at any point of the muscle is the sum of this active stress and a passive stress tensor s (p) , associated with deformation of the non-contractile matrix in which the myofilaments are embedded. Thus, we may write s = s (a) + s (p) . In the absence of body forces equilibrium requires 9 · s = 0 (Eringen, 1962), so 9 · s (p) = −9 · s (a); the divergence of the active stress acts as a body force loading the passive matrix. More importantly, the formulas of Appendix B show that s(a) is not symmetric, so that I×· s (a) $ 0 (I is the identity tensor). The quantity I×· s (a) is equivalent to a body-moment per unit volume and this moment acts on, and must be equilibrated by, the passive matrix. This suggests that the passive matrix should be regarded as a so-called ‘‘micropolar continuum’’ capable of generating unsymmetrical stresses and possibly stress-couples (moments per unit area), denoted by a second-rank tensor m(p) , that equilibrates

-     the active stress by satisfying the equation 9 · m (p) + I×· s (p) = −I×· s (a) . (Eringen, 1962). The important qualitative point being made here is that, except in special circumstances like uniaxial stretch, the active stresses produced by the cross-bridges are unsymmetric, and therefore not in rotational equilibrium. Passive stresses are required to equilibrate internally these active stresses, and perhaps this structural role is an important mechanical function of the non-contractile proteins in muscle. Exactly how this is accomplished without significantly impeding axial force generation should be an interesting question for future research. (But it should also be borne in mind that, as the preceding discussion has shown, these unsymmetric shear stresses are usually small.) These continuum mechanical consideration have a counterpart at the myofilament level. The calculations leading up to the stress formulae of Appendix B (Zahalak, 1995) show that the myosin filament is in force equilibrium under the action of the cross-bridge forces; that is, the cross-bridge forces on both halves of the filament sum to zero. But the filament is not in moment equilibrium. There are, in general, non-zero moments about the x and y axes. While these moments may be balanced to some extent by local moments exerted by the cross-bridges, which are not considered in this theory, it seems likely that passive internal resistance is necessary to achieve rotational equilibrium for the myosin filament, thus again suggesting an important mechanical role for the non-contractile proteins within the sarcomere. Conclusions A generalization of the two-state Huxley theory of contraction has been developed to account for non-axial active muscle stress and stiffness. The main ingredients of the theory are a relation between local macroscopic deformations of the muscle and the deformations of the ‘‘cage’’ of actin filaments surrounding a myosin filament, radial as well as axial cross-bridge stiffnesses, and a variation of attachment and detachment rates with lateral filament spacing. In general, non-axial deformations affect the direct stresses and give rise to active shear stresses. The theory suggests that under most circumstances non-axial stresses and stiffnesses will be small compared with the axial and, further, that the effect of non-axial deformation rates will be small compared with the axial rate. Substantial non-axial deformations may, however, reduce the axial stress and stiffness markedly. Further studies of this model, which could suggest useful new experiments, will

79

require specific quantitative hypotheses about the variations of the actin-myosin interaction rates with lattice spacing.

This research was supported by the National Science Foundation under Grant No. INT 9315743.

REFERENCES B, M. A., C, G. & C, F. (1990). Myofilament spacing and force generation in intact frog muscle fibers, J. Physiol. 430, 61–75. B, J. R. (1965). Influence of osmotic strength on cross-section and volume of isolated single muscle fiber, J. Physiol. 177, 42–57. B, B. (1988). Effect of Ca2+ on cross-bridge turnover kinetics in skinned single rabbit psoas fibers: Implications for regulation of muscle contraction, Proc. Natl. Acad. Sci. U.S.A. 85, 3265–3269. B, B. (1990). Muscle mechanics and biochemical kinetics, In: Molecular Mechanisms in Muscle Contraction, (Squire, J. M., ed.) pp. 77–142. Boca Raton, FL: C. R. C. Press. B, B. & Y, L. C. (1991). Characterization of radial force and radial stiffness in CA2+-activated skinned fibers of the rabbit psoas muscle, J. Physiol., 441, 703–718 C, A. A., H, J. A. & P, I. E. (1989). Muscle fiber lengths, pinnation angles, and deformation of aponeuratic sheets in the cat medial gastrocnemius muscle during normal movement, Soc. Neurosc. Abs., 15, 213–13. E, K. A. P. & H, J. C. (1977). The force-velocity relationship in skeletal muscle fibers at varied tonicity of the extracellular medium, J. Physiol., 269, 255–272. E, K. A. P. & R, C. (1984). Absence of plateau of the sarcomere length-tension relation in frog muscle fibers, Acta Physiol. Scand., 122, 213–216. E, K. A. P. & R, C. (1987). The sarcomere length-tension relation determined in short segments of intact muscle fibers of the frog, J. Physiol., 385, 709–732. E, K. A. P. (1979). The velocity of unloaded shortening and its relation to sarcomere length and isometric force in vertebrate muscle fibers, J. Physiol., 291, 143–159. E, G. F., L, J. & W, C. R. (1963). An x-ray and light diffraction study of the filament lattice of striated muscle in the living state and rigor, J. Mol. Biol., 6, 295–305. E, A. C. (1962). Nonlinear Theory of Continuous Media, New York: McGraw-Hill. F, J. T., S, R. M. & S, J. A. (1994). Single myosin molecule mechanics: piconewton forces and nanometer steps, Nature, 368, 113–119. F, L. E., H, A. F. & S, R. M. (1981). The relation between stiffness and filament overlap in stimulated frog muscle fibers, J. Physiol., 311, 219–249. G, Y. E. & S, R. M. (1986). The stiffness of frog skinned muscle fibers at altered lateral filament spacing. J. Physiol., 378, 175–194. G, Y. E. (1987). Measurement of sarcomere shortening in skinned fibers from frog muscle by white light diffraction, Biophys. J., 52, 57–68. G, A. M. & G, R. E. (1970). Some effects of hypertonic solutions on contraction and excitation-contraction coupling in frog skeletal muscles, J. Gen. Physio., 55, 254–274. H, J. V. (1958). The behavior of frog muscles in hypertonic solution, J. Physiol., 144, 167–175. H, P. J. & S, B. H. (1988). The analysis of cardiac function: a continuum approach, Prog. Biophys. Molec. Biol. 52, 101–164. H, A. F. (1957). Muscle structure and theories of contraction, Prog. Biophys. Biophys. Chem., 7, 255–318.

. . 

80

H, A. F. (1980). Reflection on Muscle, Liverpool: Liverpool University Press. I, A., H, Y., K, H., T, T., H, H. & Y, T. (1994). Single-molecule analysis of the actomyosin motor using nano-manipulation, Biochem, Biophys. Res. Com., 199, 1057–1063. K, M. and S, M. I., (1985). Crossbridge kinetics in chemically skinned rabbit psoas fibers when the actin-myosin lattice spacing is altered by dextran T-500, J. Musc. Res. Cell. Mot., 6, 313–332. K-W, J., (1963). Influence of osmolarity of perfusate on contractility of mammalian myocardium, Amer. J. Physiol., 204, 957–963. K, B. and M, D. (1984). The relationship between ATP hydrolysis and active force in compressed and swollen skinned muscle fibers of rabbit, Pflugers Arch., 400, 160–165. M, A. (1993). Tension transients in skeletal muscle fibers of the frog at varied tonicity of the extracellular medium, J. Musc. Res. Cell Mot., 14, 15–25. M, A. (1994). The tension response to stretch of intact skeletal muscle fibers of the frog at varied tonicity of the extracellular medium, J. Musc. Res. Cell Mot., 15, 145–157. M, I., G, Y. E. & S, R. M. (1984). Changes in the lateral filament spacing of skinned muscle fibers when cross-bridges attach, J. Mol. Biol., 173, 15–33. M, J. M. & M, R. L. (1987). Shortening velocity in skinned single muscle fibers, Biophys. J., 52, 127–131. S, J. M., (1986). Muscle: Design, Diversity, and Disease, Menlo Park, CA: Benjamin/Cummings. S, D. D., J., S, H. H., P D. P., R, J. J. & S, E. H. (1969). Fiber orientation in the canine left ventricle during diastole and systole, Circ. Res. 24, 339–347. T, R. T. & S J. M. (1973). Myosin content and filament structure in smooth and striated muscle, J. Mol. Biol., 77, 279–290. Y, L. C. & P, R. J. (1990). Equatorial x-ray diffraction studies of single skinned muscle fibers, In: Molecular Mechanisms in Muscle Contraction (Squire, J. M., ed.), pp. 265–286. Boca Raton, FL: C. R. C. Press. Z, G. I. (1995). A three-dimensional generalization of the Huxley cross-bridge theory, Technical Report No. WU 9523, Dept. Mech. Eng., Washington University, St. Louis, U.S.A.

APPENDIX A The relation between the macroscopic deformation of the muscle and that of the actin cylinder is derived in this Appendix. Referring to Fig. 2 let XI represent material coordinates in the undeformed configuration; and xi spatial coordinates in the deformed configurations; in both configurations the three-direction is along the fibers. Define also a material projection tensor, P = I − E3 E3 , on planes perpendicular to the undeformed fiber direction, and a corresponding spatial projection tensor, p = I − e3 e3 , on planes perpendicular to the deformed fiber directions; I is the identity tensor, and E3 and e3 are unit vectors in the undeformed and deformed fiber directions, respectively (see Fig. A1). The equation of a small circular cylinder, of radius dR, parallel to the fibers in the undeformed configuration is (dX − dX · E3 E3 ). (dX − dX · E3 E3 ) = dR 2 which

e1

(a)

λ x dR ds

dst

e3

Deformed

λ y dR

E1

(b)

dR dS

E3

∆ dst· X·P

∆ dst· X

Undeformed (c)

λ xa

y Undeformed

λ ya

a

ro

ex

θ ey

Deformed x

ro(θ) = aλ xcos(ω)ex + aλ y sin(ω)ey λ x tan(θ) = λ y tan(ω) F. A1. (a) Elliptic deformed cylinder with axis parallel to the current direction of the muscle fibers. (b) Circular-cylinder undeformed image of the cylinder in (a), with axis parallel to the original undeformed direction of the muscle fibers. (c) Comparison of the elliptic deformed and circular undeformed cross-sections of the cylinder.

may be written as (dX · P) · (dX · P) = dR 2, and further as dX · P · dX = dR 2

(A.1)

because, P is symmetric, and further P · P = P. But dX = dx · 9X = X9 · dx, where 9X = XI, j ej EI is the material position gradient (the comma represents partial differentiation with respect to xj ). Therefore, dx · 9X · P · X9 · dx = dR2

(A.2)

is the deformed spatial image of the original circular cylinder. Because both coordinate systems have been chosen with the three-direction along the fibers, which are material lines consisting of the same particles, it follows that X1.3 = X2.3 = 0 for the muscle element at the origin of coordinates. This in turn implies that undeformed straight lines parallel to the X3 axis are mapped into straight lines parallel to the x3 axis. In

-     particular a brief calculation (Zahalak, 1995) shows that 9X · P · X9 = cˆ = cˆab ea eb

(A.3)

with cˆab = XD,a XD,b , where Greek indices range over (1, 2), and summation over repeated indices is assumed. The two-dimensional spatial tensor cˆ is the Cauchy deformation tensor of the projection of the deformation on planes perpendicular to the fiber direction. It is related to the projection of the full Cauchy deformation tensor, c = 9X · X9, by cˆ = p · c · p − X3,a X3,b ea eb .

(A.4)

Thus, the deformed image of the original circular cylinder is, according to eqn (A.2) cx ab dxa dxb = dR 2.

cˆxx x 2 + 2cˆxy xy + cˆyy y 2 = a 2.

(A.6)

The simplest representation is obtained if the x and y axes are chosen along the principal axes of the ellipse, which are the eigendirections of cˆ. Further it is well known (Eringen, 1962) that the eigenvalues of the Cauchy deformation tensor are the inverse squared stretch ratios; in this case the stretch ratios are the ratios of deformed to undeformed principal axes, denoted by lx and ly . With this particular orientation of the x1 and x2 axes, which is employed throughout the paper, the equation of the actin cylinder’s cross section becomes −2 2 2 2 l−2 x x + ly y = a ,

As t=

dr dr dv dv = = {−alx sin vex + aly cos vey } ds dv ds ds

and ds2 = dx 2 + dy 2 = a 2{lx2 sin2 v + ly2 cos2 v}dv 2, it follows that t=

(A.7)

or, in parametric form, x = alx cos vt and y = aly sin vt [see eqn (14)]. To calculate how the deformation affects the spacing of the axial generators on the actin cylinder consider Fig. A1. If dst is a differential vector measuring the spacing between two generators on the deformed cylinder, then its image on the undeformed cylinder is dst · 9X, and the undeformed spacing is given by the projected vector dst · 9X · P. Therefore, the magnitude of the undeformed spacing, dS, is dS 2 = ds 2t · 9X · P · X9 · t = ds 2t · cˆ · t. (A.8)

−lx sin vex + ly cos vey . zlx2 sin2 v + ly2 cos2 v

Given the choice of −2 cˆ = l−2 x ex ex + ly ey ey , so

coordinate

ds = (t · cˆ · t)−1/2 = zlx2 sin2 v + ly2 cos2 v dS

(A.9) system,

(A.10)

which appears in the main body of the paper as eqn (15).

(A.5)

As dx3 does not enter the above quadratic form it represents, as expected, an elliptic cylinder. If it is assumed that the macroscopic deformation is transmitted to the actin cylinder in cross-bridge space, then the cylinder’s cross-section is described by the equation

81

APPENDIX B This Appendix lists expressions for the nine components of the active Cauchy stress tensor in terms of the bond-distribution function, n(z, u, t). The details of the calculations are available in Zahalak (1995), and the following expressions use the notation of eqn (32) and the normalization constant s0 = (ms0 kz h 2/2l). (s11 /s0 ) =

0 10 10 1

1 p a ly lz 2 s0

a h

kr kz

× {−p/2, p/2=1, 0=1, 1, 0 − p/2, 3p/2=1, 0=1, 1, 0} (s22 /s0 ) =

0 10 10 1

1 p a lx lz 2 s0

a h

kr kz

× {0, p =0, 1=1, 1, 0 − p, 2p =0, 1=1 1, 0} (s33 /s0 ) =

1 0, 2p =0, 0=0, 0, 1 lx ly

(s12 /s0 ) =

1 p a ly lz 2 s0

0 10 10 1 a h

kr kz

× {−p/2, p/2=0, 1=1, 1, 0 − p/2, 3p/2=0, 1=1, 1, 0} (s21 /s0 ) =

0 10 10 1

1 p a lx lz 2 s0

a h

kr {0, p =1, 0=1, 1, 0 kz

− p, 2p =1, 0=1, 1, 0}

. . 

82 (s23 /s0 ) =

01

1 p a {0, p =0, 0=0, 0, 1 lxlz 2 s0

B2212 = B2221 = B2122

0 10 1

1 a lx ly h

kr 0, 2p =0, 1=1, 1, 0 kz

0 10 1

1 a (s31 /s0 ) = lx ly h

(s13 /s0 ) =

kr 0, 2p =1, 0=1, 1, 0 kz

01

kr kz

× {0, p =1, 2=0, 1, 0 − p, 2p =1, 2=0, 1, 0}

− p, 2p =0, 0=0, 0, 1} (s32 /s0 ) =

0 10 1

1 p a lx lz 2 s0

01

B3333 =

lz s0 0, 2p =0, 0=0, 0, 0 lx ly 2a

B3331 =

1 0, 2p =1, 0=0, 1, 0 lx ly

B3332 =

1 0, 2p =0, 1=0, 1, 0 lx ly

B1222 =

1 p a ly lz 2 s0

1 p a {−p/2, p/2=0, 0=0, 0, 1 ly lz 2 s0

− p/2, 3p/2=0, 0=0, 0, 1}

0 10 1

kr {−p/2, p/2=0, 3=0, 1, 0 kz

− p/2, 3p/2=0, 3=0, 1, 0} APPENDIX C The non-zero components of the active stiffness tensor have been calculated in Zahalak (1995) and are listed below. The stiffness tensor Bijpq is defined by eqn (46), and the notation of eqn (32) is employed.

0 10 1

1 p a B1111 = ly lz 2 s0

kr {−p/2, p/2=3, 0=0, 1, 0 kz

B2111 =

kr {0, p =3, 0=0, 1, 0} kz

− p, 2p =3, 0=0, 1, 0} B2333 =

1 p {0, p =0, 0=0, 0, 0 lx 4

−p, 2p =0, 0=0, 0, 0}

− p/2, 3p/2=3, 0=0, 1, 0}

0 10 1

1 p a B1122 = B1212 = B1221 = l y lz 2 s 0

kr kz

× {−p/2, p/2=1, 2=0, 1, 0 − p/2, 3p/2=1, 2=0, 1, 0} B1112 = B1121 = B1211 =

0 10 1

1 p a lx lz 2 s0

kr kz

× {−p/2, p/2=2, 1=0, 1, 0

01

1 p a {0, p =1, 0=0, 1, 0 lx lz 2 s0

− p, 2p =1, 0=0, 1, 0} B2332 =

0 10 1

1 p a l y lz 2 s 0

B2331 =

01

1 p a {0, p =0, 1=0, 1, 0 lx lz 2 s0

− p, 2p =0, 1=0, 1, 0} B3211 = B3112 = B3121 =

− p/2, 3p/2=2, 1=0, 1, 0} B2211 = B2112 = B2121 =

0 10 1

p 1 a 2 lx lz s0

kr kz

× {0, p =2, 1=0, 1, 0 − p, 2p =2, 1=0, 1, 0}

0 10 1

1 p a B2222 = lx lz 2 s0

B3222 =

− p, 2p =0, 3=0, 1, 0}

01

1 kr 0, 2p =0, 3=0, 1, 0 lx ly kz

B3212 = B3221 = B3122 =

kr {0, p =0, 3=0, 1, 0 kz B3111 =

01

1 kr 0, 2p =2, 1=0, 1, 0 lx ly kz

01

01

1 kr 0, 2p =1, 2=0, 1, 0 lx ly kz

1 kr 0, 2p =3, 0=0, 1, 0 lx ly kz

-     B1333 =

1p {−p/2, p/2=0, 0=0, 0, 0 ly 4

This gives d = (g˙ /2)(e1 e3 + e3 e1 ) and V = (g˙ /2) (−e1 e3 + e3 e1 ). As, in this simple case, the x3 axis remains parallel to the X3 axis, and any transverse direction is an eigendirection of cˆ, the spatial x1 coordinate axes can be chosen parallel to the material XI axes. Then the spatial fiber coordinate system does not rotate in the laboratory frame, so V = 0, and the relative velocity gradient has only one non-vanishing component, D = d + V − V = g˙ ez ex . This gives, by eqns (12), Dzr = g˙ cos u, Dur = 0, and Dzz = 0. Thus, the rate equation, eqn. (20) becomes

− p/2, 3p/2=1, 0=0, 1, 0} B1331 =

01

1 p a {−p/2, p/2=1, 0=0, 1, 0 ly lz 2 s0

− p/2, 3p/2=1, 0=0, 1, 0} B1332 =

83

01

1 p a {−p/2, p/2=0, 1=0, 1, 0 ly lz 2 s0

− p/2, 3p/2=0, 1=0, 1, 0}

1n 1n + 2meg˙ cos u = f(z) − [f(z) + g(z)]n(z, u, t). 1t 1z (D.2)

APPENDIX D Suppose that a block of muscle starts in an undeformed configuration, with r0 = a, and is subjected to a shearing deformation parallel to the fibers at a constant rate g˙ , as shown in Fig. 7(a). The macroscopic deformation is described by the equations x1=X1 , x2 = X2 , x3 = X3 + g˙ X1t or inversely X1 = x1

X2 = x2

X3 = x3 − g˙ x1 t.

Note that because r remains constant in this problem the rate functions f and g are independent of u and t. Equation (D.2) admits a steady-state solution, n(z, u), if (1n/1t) is set equal to zero. This results in an ordinary differential equation on n, with u as a parameter. Indeed, this equation is identical with that of Huxley’s 1957 paper if the axial velocity (−sV/2) in the latter is replaced by 2meg˙ cos u. Choosing piecewise linear rate functions as in Huxley’s original theory (modified slightly to give a more realistic stretch behavior), f(z) = f1 z{H(z) − H(z − 1)} and g(z) = g2 H(−z)+ g1 z{H(z) − H(z − 1)} + g3 zH(z − 1), where H(z) is the Heaviside step function and f1 and g1 are constants. The steady-state solution is easily found to be the following:

(D.1)

From this it follows immediately that l = X3,3 = 1, XD,a = dDa , cˆab = dab , so that lx = ly = 1. This means that r = 1, so that the actin cylinder cross-sections remain circular and undistorted, and further n˜ = r 2/ lx ly = 1, if it is assumed that all the cross-bridges participate in the contraction during the deformation (a = 1). The absolute particle velocity associated with the deformation eqns (D.1) is simply u = gx1e3 . −1 z

for (−p Q u Q −p/2) or (p/2 Q u Q p)

F 1 2z Gn0 1 − exp d cos u exp −cd cos u n(z, u) = g Gn0 1 − exp − 1 (z2 − 1) d cos u f

0 0

6 6

71 6

71

7

for z Q 0 for 0 Q z Q 1

and zero otherwise, whereas for (−p/2 Q u Q p/2)

F z2 Gn0 1 − exp −d cos u n(z, u) = g 1 Gn0 1 − exp − 1 exp − (z 2 − 1) d cos u bd cos u f

0 0

6 6

71 71 6

7

for 0 Q z Q 1 for 1 Q z

and zero otherwise. The parameters appearing in the above equations are defined as n0 =

01

f1 2a g , d = , f1 + g 1 h f1 + g1

c=

g1 + f1 , g2

b=

g1 + f1 . g1 + g3

. . 

84

An examination of the formulae for the stresses listed in Appendix II shows that in this case there are five out of nine non-vanishing stress components, which may be expressed as s33 = s0 n0 S33 , s13 =

s31 =

01

p a snS , 2 s0 0 0 13

n0 S11 =

0 10 1 a h

kr (1 − r*)s0 n0 S31 , kz

n0 S22 =

0 10 10 1

s11 =

p a 2 s0

s22 =

p a 2 s0

a h

kr (1 − r*)s0 n0 S11 , kz

0 10 10 1 a h

kr (1 − r*)s0 n0 S22 kz

where the Sij are the angular averages of the zero-and first-order axial moments of n, specifically n0 S33 =

n0 S13 =

n0 S31 =

1 p

6g

p

7

1 p

6g

p/2

Q1 du ,

0

Q1 du −

g

7

p

Q1 du ,

p/2

0

1 p

6g

p

1 p

6g

p/2

1 p

6g

p

Q0 cos u du

0

Q0 cos u du −

g

p

7

Q0 cos u du ,

p/2

0

7

Q0 sin u du ,

0

a with Ql (u; d ) = f−a z ln(z, u; d ) dz. These integrals, defining the moments Q0 and Q1 , can be evaluated in closed form for the steady-state n(z, u) given above (Zahalak, 1995) but to conserve space the results will not be listed. Graphs of Q0 and Q1 are given in Fig. 7(b) and (c); these can be interpreted, respectively, as showing the angular distribution of axial stiffness and force. To obtain the stress coefficients, Sij , these moment were integrated numerically with respect to u, and the results are shown in Fig. 8(b) and (c) as functions of the dimensionless shear rate, d . Finally, the formulae given in Appendix C were used to compute the stiffness coefficients. Restricting attention to coefficients involving axial forces and displacements we may write

F a 2 G p K1331 0 F G DT13 J G G s0 G G G 2 1 a G DT23 G = G 0 p K2332 s0 s0 G G G n0 s 2h 0G f DT33 G j G G 2 a K3331 0 s0 f

01

01

p a K 2 s0 1333

01

01

7

0

01

K3333

J G GF g31 J GG G GG g32 G GG G Gf g33 j G j

where pn0 K3333 =

g

p

Q0 du, pn0 K2332 =

0

pn0 K1331 =

g

p/2

0

g

p

Q0 sin u du,

pn0 K3331 =

0

Q0 cos u du −

g

p

p/2

Q0 cos u du,

g

p

Q0 cos u du,

0

pn0 K1333 =

g

p/2

0

Q0 du −

g

p

Q0 du.

p/2

The stiffness coefficients were again evaluated by numerical integration of the moments illustrated in Fig. 7(b) and (c), and their variation with shear rate is displayed in Fig. 9(b) and (c).