IV. Implementing the Model - Washington State University

eatablesurveyorInternet και Εφαρμογές Web

14 Δεκ 2013 (πριν από 3 χρόνια και 7 μήνες)

206 εμφανίσεις


1

VARIABLY ACTIVATED S
ARCOMERE MODEL FOR M
USCLE
CONTRACTION STUDIES
AND APPLICATIONS

Kenneth B. Campbell
1,2
, Maria V. Razumova
3
, Robert D. Kirkpatrick
1
,

Bryan K. Slinker
1

1
Department of VCAPP

2
Department of Biological Systems Engineering

Washington State Un
iversity

Pullman, WA 99163

3
Department of Physiology

University of Wisconsin
-
Madison

1300 University Ave

Madison WI 53706




Running Title


Sarcomere Model for Muscle Studies.



Supported by NIH Grant RO1 HL21462
-
20



Address for correspondence:

Kenneth
B. Campbell

Department of VCAPP

Washington State University

Pullman, WA 99163


Phone: (509) 335
-
8011

FAX (509) 335
-
4650


E
-
mail: cvselkbc@vetmed.wsu.edu



2

Abstract

A mathematical model of muscle contraction is constructed from myofilament structural featu
res
and kinetic processes that take place within the muscle sarcomere. Sliding filament processes are
considered with respect to their effect on filament overlap and potential number of interaction sites
between thick and thin filaments. The kinetics of
thick
-
thin filament interaction are considered in terms
of the kinetics of thin filament regulatory units and of crossbridge, XB, cycling. A differential equation
expressing the dynamics of the average distortion among XB in each of two attached XB states

allows
formulation of sarcomere force production in terms of stiffness and distortion. With this expression for
force, interdependencies within the myofilament kinetic system such as various forms of cooperativity
and of length
-
dependent kinetics are int
roduced into the model. The full model consists of five
differential equations: three to describe dynamics of XB states and two to describe dynamics of XB
distortion. In addition, the model contains the following nonlinear relationships: 1) the relation
of
filament overlap

to
sarcomere length
; 2) saturation effects due to finite numbers of XB within
the
overlap zone

and finite number of binding sites for activator Ca
2+

on Troponin C; 3) the dependence of
kinetic rate coefficients on distortion, neighbor i
nteractions, and
sarcomere length
; and 4) the
appearance of products of state variables in the stiffness
-
distortion expression for force. In its full
nonlinear form, the model contains 22 parameters: 16 associated with kinetic coefficients; 5 associated
w
ith structural factors; and one linking structural and kinetic factors. A critique of the model is given
with respect to assumptions made in its construction. Suggestions for model implementation are also
given.





3

I.

Background

Mathematical modeling of m
uscle has been central to understanding the biophysics of muscle
contraction and has also been important in integrating muscle dynamics into body movement functions.
The pursuit of these two objectives has largely followed two different model
-
building pat
hs. Those
interested in contractile biophysics build models based on molecular mechanisms of contraction (often
patterned after the landmark work of A.F. Huxley [1957]) whereas those interested in integrating
knowledge of muscle into body movement functio
n most often base their models on phenomenologic
behaviors [Brown et. al., 1996; Curtin et. al., 1998; Dorgan and O’Malley, 1997; Kaufman et. al., 1991;
Krylow and Rymer, 1997; Schoechting and Flanders, 1997; Schultz et. al., 1991; Shue and Cargo, 1998;
We
xler et. al., 1997; Williams et. al., 1998; Winters, 1995; Zahalak, 1992; Zajac, 1989]. In the latter
case, the principal phenomenologies used are isometric length
-
tension and isotonic force
-
velocity
relations where the force
-
velocity relationship is repr
esented by the classical hyperbolic relation of A.V.
Hill [1938]. Unfortunately, isometric length
-
tension and isotonic force
-
velocity relations are not general
descriptors of muscle behavior and pertain only to the specific experimental conditions in whic
h data for
the respective relation have been collected. In fact, the force, produced by a muscle shortening against
an arbitrary mechanical load, depends on the histories of both activation and shortening/force production
as well as on the instantaneous v
elocity and length. Length
-
tension and force
-
velocity relations do not
possess expressions for history dependent effects and thus, cannot serve as an adequate basis for a model
of a dynamic system whose behavior is heavily influenced by its immediate hist
ory. Those who use the
phenomenologic approach do so largely by default because it is generally accepted that models based on
molecular mechanisms have not reached either a status of adequate predictive capability or of user
friendliness that allows ready

application to practical problems [Winters and Stark, 1987; Zajac, 1989;
van den Bogert et.al., 1996; Wu et.al., 1997].


Our interest in presenting the current model is motivated by a desire to document a concise,
contractile
-
mechanism inspired sarcomer
e model that may be used as a central component of a muscle
model to: 1) clarify the role of the various contractile regulatory processes in normal, adapted, and
diseased muscle function and 2) integrate muscle dynamics into a broad range of application ta
sks that
involve muscles with different dynamic attributes. Towards this end, we have proposed an approach
that utilized the muscle sarcomere, with essential properties of stiffness and distortion, as the elemental
building block of a muscle model (Razumo
va et.al., 1999). Both sarcomeric stiffness and distortion
dynamics arose from a cycling
-
crossbridge, XB, kinetic scheme. In that study, the model sarcomere,
under conditions of constant calcium activation, was evaluated with respect to:
1
) the frequenc
y
-
dependent complex stiffness,
2
) time course of transient force response to step and ramp changes in
length, and
3
) the force
-
velocity relationship. We found that the model reproduced all relevant
experimentally
-
measured facets of these responses and pro
vided a sound basis for interpreting
mechanisms of these phenomenologic behaviors. Further aspects of the behavior of this stiffness
-
distortion sarcomere model were investigated in a second study [Razumova et. al., 2000] where the
effects of neighbor inte
ractions among thin filament regulatory units and XBs were evaluated with
respect to their effect on:
1
) magnitude, steepness, symmetry, and position of the steady
-
state force
-
pCa
relationship and
2
) the characteristic time of the dynamic force redevelopme
nt process. We found that
different kinds of neighbor interactions had profound, unexpected, and unique effects on these measures
of muscle function. From these two studies, the groundwork was laid for a model that could be used to
describe, predict, and

explain the complexities of the many behavioral features of muscle contraction.

Mathematical models of imperfectly understood biologic processes can never be considered
finished. Thus, the current sarcomere model is a work in progress that has achieved
some level of
validity and, therefore, some degree of usefulness. In order to be tractable, the current model has

4

incorporated many assumptions and simplifications that some may argue are poorly considered. While
remaining respectful of such criticism,
we defend the decisions that were made as necessary at this time
as an early step towards realizing an objective of connecting overall contractile system behavior with
underlying processes. Recognizing that these simplifying decisions will be modified as
more
information is obtained, we present the model as it currently exists and then critique it with respect to
some of the most far
-
reaching assumptions to indicate where future improvements may be made.


II.

Contractile System Model

A.

Muscle Structure and Slidi
ng Filaments

Skeletal muscle fibers may be several cm long and are typically 100


200

m diameter. These
fibers contain bundles of parallel, 2

m diameter myofibrils. Myofibrils are constructed of protein
filaments that are organized into a chain of se
ries
-
connected sarcomere units; each sarcomere is on the order
of 2 to 3

m long and is bounded by a pair of narrow Z
lines (Fig. 1). Within the sarcomere is a central, optically
dense A
-
band flanked by translucent I
-
bands. The I
-
bands
are transected by t
he Z
-
lines. The A
-
band contains
approximately 540 parallel thick filaments (of length
L
M

~1.6

m) made primarily from myosin. Thin filaments
(constructed on an F
-
actin backbone and of length
L
A

~1.1


1.3

m) originate from the Z
-
line and extend through

the
I
-
band and into the A
-
band. A half sarcomere is that part
between the midline of the A
-
band (the M
-
line) and the Z
-
line. Within a half sarcomere there are twice as many thin
filaments as there are thick filaments [Squire, 1981].

Individual myosin

molecules within the thick
filament consist of a head region connected by a hinge to a flexible tail. The head and a portion of the
tail project out of the thick filament into the lattice space that accommodates the interdigitating thin
filaments. This
projection is called the crossbridge, XB. Force is generated as a result of interaction
between thick and thin filaments through XB attachment to the actin binding site on the thin filament.
Because XB along a thick filament are arranged in parallel with

one another when attached to thin
filaments and because thick filaments within a sarcomere are also parallel to each other, attached XB
within a half sarcomere are all in parallel. Thus, the force generated in a half sarcomere equals the sum
of forces gen
erated by all force
-
generating XB in that half sarcomere. The two halves of a sarcomere
are mirror images, such that the force generated in one half just balances the force generated in the other
half and there is no tendency for the thick filaments to mo
ve preferentially towards either Z
-
line.


XBs bear force only when they can attach to the thin filament. Attached XBs are restricted to
that part of the filament overlap zone exclusive of the bare zone of the thick filament,
B

(0.2

m), and
exclusive of the double
-
overlap zone in which the ends of the thin filaments from opposite sides of the
sarcomere overlap [Trombitas and Tigyi
-
Sebes, 1984; 1985]. The zone in the half sarcomere,
Ov
, in
which force
-
bearing XBs may form between thic
k and thin filaments varies with sarcomere length, SL,
according to three overlap conditions as shown in Fig. 2:

Condition I:
Ov

increases with
SL

according to



(1)

Condition II:
Ov

does not change with
SL

and equals


5




(2)

Condition III:
Ov

decreases with
SL

according to




(3)

If
R
T

is the total number of XBs on all thick
filaments in a half sarcomere (roughly, 540 thick
filaments/sarcomere x 300 XB/half
-
thick filament =
1.62 x 10
5

XB/half sarcomere), then, the number of
XBs that are in
Ov

and capable of attaching to the
thin filament is given by


(4)

Considering the condition of optimal
overlap, i.e., Condition II where
R
Ov

= R
T
, and
assuming

F/

L

of a single

XB to be 3 x 10
-
4

N/m
(based on reports of Finer et. al. [1998] that indicate
3 x 10
-
12

N force with a step on the order of 10 nm),
then, when in rigor with all XB attached, the

F/

L

of the collection of all XB in a half sarcomere would
be 48.6 N/m. Fur
ther, assuming the cross
-
sectional
area of a representative sarcomere to be 3 X 10
-
12

m
2

and the length of the half sarcomere for

L

normalization to be 10
-
6

m, the half sarcomere in
rigor and at optimal overlap would exhibit a
stiffness of 16.2 MN/m
2
.

Th
is theoretical rigor stiffness compares
favorably with values obtained experimentally in
rabbit psoas muscle (15
-
20 MN/m
2
) using sinusoidal
perturbation techniques [Kawai and Brandt, 1976; Murphy et. al., 1996]. Thus, from a purely structural
standpoint,
we have satisfactorily accounted for all stiffness elements and, by commonly accepted
assumptions, all potential force
-
generating XB elements.

B.

Thin Filament Regulatory Kinetics and Crossbridge Cycling Kinetics

Sarcomere force generation depends on the numb
er of XB in a force
-
bearing state, which, in turn
and for a given
Ov
, depends on the kinetics of Ca
2+
-
activation, thin filament regulation and XB cycling.
These interacting kinetic processes are represented in Fig. 3. Thin
-
filament regulation occurs by
s
witching the tropomyosin
-
troponin regulatory unit, RU, between an
off

state that does not permit XB
attachment and an
on

state that does. This
on
-
off

switching is strongly affected by whether Ca
2+

is
bound to the troponin C subunit (TnC) of the RU troponi
n complex. When Ca
2+

is not bound to TnC,
the preferred RU state is '
off
'; when Ca
2+

is bound to TnC, the preferred RU state is '
on
'. However,
whether or not Ca
2+

is bound, there remains a finite probability that the RU is in a non
-
preferred states.
Th
us, we identify four rate coefficients for RU
on
-
off

state transitions:
, where the
superscript indicates whether Ca
2+

is bound to the RU. In accord with our understanding of the
influence of Ca
2+
, the following inequalities hold:

;
; and
.


6

If Ca
2+

binding to, and dissociation from, TnC is fast relative to all other kinetic events

and
independent of other steps in the contraction sequence, then Ca
2+

binding ma
y be considered to be in
equilibrium and for time
-
varying Ca
2+
,
Ca(t)
, we may write,


(5)


(6)

where the ratio,
Ca(t)/[Ca
50

+ Ca(t)],

on the right
-
hand side is the TnC Ca
2+

binding isotherm for a
single low
-
aff
inity Ca
2+

binding site and
Ca
50

(
) is the calcium concentration of half saturation
of the TnC Ca
2+

binding site [Campbell, 1997].

For two low
-
affinity TnC Ca
2+

binding sites (binding to the first site precedes binding to the
second
site), the binding isotherm for occupancy of both sites is given by


where
K
1

is the Ca
2+

concentration of half occupancy of the first site and
K
2

is the Ca
2+

concentration of
half occupancy of the second site. This expression may b
e substituted for the single occupancy binding
isotherm in Eqs. (5, 6) as appropriate.




Figure 3


Three kinetic processes in myofilament interaction showing: 1) Ca
2+

binding to the tropomyosin
-
Troponin
complex that acts o
f the thin
-
filament regulatory unit (RU); 2) switching between ‘off’ and ‘on’ states of the RU; and 3)
myosin cross bridges (XB) cycling between attached and detached states. The forward,
k
+
, and backward,
k
-
, Ca
2+

binding
constants are taken to be large
with respect to other kinetic constants and independent of other steps in the contraction
process such that Ca
2+

binding is always in equilibrium. The three connected circles represent the thin filament. The bar
along the thin filament represents the RU.

If in the down position, the RU is
off

and XBs cannot attach. If in the up position,
the RU is
on

and XBs can attach. The rate coefficients regulating the transition between ‘on’ and ‘off’ states,
k
on
(Ca)

and
k
off
(Ca)
, are functions of Ca
2+

bound to RU

as given in the text. The darkened ellipse with the coiled tail represents the XB
which may be in one of three states: detached,
D
, attached pre
-
power
-
stroke,
A
1
, or attached and post
-
power
-
stroke,
A
2
.


A complication in thin
-
filament activation is that
the
RU
off

-

RU
on

transition is affected by the
state of the XB. As shown in Fig. 3, there are three XB states: detached,
D
; attached pre
-
power stroke,

7

A
1
; and attached post
-
power stroke,
A
2
. If the XB is in either the
A
1

or
A
2

attached states, the
RU
on

m
ay
not turn '
off
' and transition to
RU
off

must await the formation of the XB detached state,
D
. For
simplicity, we assume a 1:1 stoichiometry between RU and XB. This is done even though it is
recognized that the span of a single RU is over 7 actin monome
rs and, within this span, as many as 3 XB
on the thick filament may attach [T.L. Hill, 1985]. However, considering that RU:XB stoichiometry
greater than 1:1 introduces mathematical complications that detract from our current overall objectives.
Alternat
ive stoichiometries represent an item for future model refinement.

Within the XB cycle (Fig. 3), we refer to the transition between state
D

and state
A
1

as the
attachment step; the transition between states
A
1

and
A
2

as the power stroke; and the transition

between
state
A
2

and
D

as the detachment step. Because chemo
-
mechanical energy transduction takes place in
the power stroke, the pre
-
power stroke
A
1

state does not bear force except in non
-
isometric conditions
when the sarcomere undergoes length change.

In contrast, the post
-
power stroke
A
2

state

is force
-
bearing in both isometric and non
-
isometric conditions. The attachment step and the power stroke are
treated as reversible but the detachment step is treated as irreversible.

Steps in the XB cycle are
regulated by the rate coefficients
f, f', h, h'
, and
g
, which determine the rate
of change between states. With these, inspection of Fig. 3 allows the following differential equations to
be written for state transitions:


(7)


(8)


(9)

By applying conservation principles to the number of sites potentially available for actin
-
myosin
interaction (i.e., those sites in the filament overlap zone as given by Eq. (4)),



(10)

and
all variables in Eqs. (7


9) are specified.

It will be convenient in the developments that follow to define fractional distribution of states among
all sites in terms of probability of finding a given site in a given state. Assuming random distributions,

the probability within the population of a given site being in a given state may be written as:


is the probability that a given site possesses an RU that is turned ‘on’.


is the probability that a given sit
e possesses an RU that is turned ‘off’.


is the probability that a given site possesses a XB in the force
-
bearing,
A
2

state.


C.

Cross Bridge Distortion


XB are assumed to be linearly elastic and, therefore, generate force when distort
ed. XB are
distorted by two mechanisms:
i
) the power stroke and
ii
) filament sliding. We treat these as independent
mechanisms. In the power stroke, chemo
-
mechanical energy transduction causes conformational change
in the head around the hinge region an
d introduces a distortion,
x
0
, in the spring
-
like portion of the XB
molecule (Fig. 4). This
x
0

distortion is induced during isometric contraction as well as during conditions
when sarcomere length is undergoing change and filaments slide past one another.



8

Figure 4


XB distortion as a result of power stroke. Prior to the power stroke, the XB in the
A
1

state is undistorted and does
not bear force. The coiled part of the XB of length
l
0

is taken to be slack. Chemo
-
mechani
cal energy transduction of the
power stroke results in conformational change of the head and
x
0

distortion of the tail as the XB enters the
A
2

state. The
introduction of this
x
0

distortion causes the elastic XB to become force generating.


During filament

sliding as when sarcomeres are stretched, the sarcomere may be thought to
undergo repeated, sudden incremental stretches of 2

. In the half sarcomere, a sudden stretch causes an
instantaneous sliding of the thick and thin filaments past one another by an

incremental distance,

, and
this causes a shear
-
like distortion in the XB. Thus when filaments slide, all XBs in both
A
1

and
A
2

attached states undergo an instantaneous distortion equal to

Fig. 5). Incremental distortion,

,
induced by filament slid
ing is not sustained but is dissipated according to the rate at which the distorted
XBs are replaced by newly
-
formed, undistorted XBs and according to the persistence of the velocity of
stretch.


State

Pre
-
stretch

Post
-
stretch


Pre
-
power
Stroke

A
1




Post
-
power
Stroke

A
2



Figure 5


XB distortion as a result of filament sliding. An incremental stretch of the half sarcomere,

, causes the thick and
thin filaments to slide past one another. XBs in the attached
A
1

and
A
2

states both undergo shear and experience consequent
distortion.


Stretch =


Stretch =



9

Let
x
1
(t)

be the time
-
varying average distortion among
A
1

XBs and
x
2
(t)

be the time
-
varying
av
erage distortion among
A
2

XBs. The differential equations for
x
2
(t)

may be derived by performing a
macroscopic, distortion balance over all parallel XBs in the
A
2

state. The collective distortion among
the parallel
A
2

XBs is given by


(11)

where the upper case
X
2
(t)

is used to designate the collective distortion and the small case,
x
2
(t)

,
designates the average distortion among the
A
2
(t)

XBs.
X
2
(t)
at some
t +

t
, can be written as

X
2
( t +

t)

=
X
2
(t) + [added distortion due to shear from change in SL over

t]

+ [added distortion due to formation of new A
2

XBs via power stroke over

t]

-

[lost distortion due to transition of distorted A
2

XBs to other states o
ver

t]

(12)

where:

[added distortion due to shear from change in SL over

t] =

[(externally imposed

SL/2) ( # A
2

XBs existent at t)] =

(

SL/2) A
2
(t)


(13)

[added distortion due to formation of new A
2

XBs over

t]=

[(number of newly formed A
2

XBs)(disto
rtion of newly formed XBs)]=

[h A
1
(t) ]

t (x
0

+ x
1

+

SL/4)

(14)

[lost distortion due to transition of distorted A
2

XBs to other states over

t] =

[(number of A
2

XBs lost)(average distortion of these lost A
2

XBs)]=

[(g + h') A
2
(t) ]

t x
2
(t)

(15)


It is a
ssumed that XBs enter the
A
2

state from the
A
1

state via the power stroke that induces
x
0

distortion.
In addition, because we treat the two causes of distortion as independent events, these XBs just entering
the
A
2

state also bring with them whatever dist
ortion they may have possessed in the
A
1

state prior to the
transition to
A
2
.

Substituting Eqs (13


15) into Eq. (12) gives:



(16)

Rearranging


(17)

Taking the limit as

t


0
, yields


(18)

Now,

may be eliminated by noting that differentiation of Eq. (11) yields


(19)

Equating (11) and (19), making appropriate substitutions for
fro
m Eq. (9), and solving for
gives the desired differential equation,


(20)


10

In like manner, it can be shown that


(21)

where XBs enter the
A
1

state from
D

with no distortion but return to
A
1

from
A
2

via the reverse power
stroke with residual distortion not removed by loss of
x
0
.

Eqs. (20
-

21), describing the dynamic changes for average distortion of the post
-
power stroke
and pre
-
power stroke states, are the key developments that allow a sim
plified model of sarcomeric
dynamics. We know of only one other approach [Thorson and White, 1983] that approximates that given
here. However, the results from derivations in Thorson and White apply only to small
-
amplitude
changes in
SL

whereas those give
n in Eqs. (20) and (21) apply generally to all length perturbations of
any amplitude.

D.

Stiffness
-
Distortion Concept


In general, the force generated by all XB within a half sarcomere is given by the sum of forces
from all XB in
A
1

and
A
2

states as


(22)

where


is the stiffness of a single XB (= 3 x 10
-
4

N/m). Because


is a constant and equal for all XBs
in attached states, the

A
1,2
(t)

products represent the stiffness of the respective attached
-
state
populations,



(23)



(24)

Thus, Eq. (22) may be rewritten in a stiffness
-
distortion format as


(25)


E.

Variable
-
dependent kinetic coefficients

The preceding development did not take into account consideration of the poss
ible dependence
of state transition kinetic coefficients (
f
,
f’
,
h
,
h’
,
g
,
k
+
,
k
-
,
k
on

, and
k
off

) on state variables or system
mechanical state. In fact, these coefficients may have such a dependence leading to nonlinearities that
have fundamentally imp
ortant roles in contraction dynamics. The functional form that expresses the
dependence of a kinetic coefficient on a variable may be taken either from fundamental chemical
-
kinetic
considerations or from phenomenology observed from experiments. In our si
mplified, population
-
average model, we combine both approaches and rely on chemical kinetic expressions as a starting point
for functional relations that then may be altered to accommodate experimental observations. It is
understood that the need for such

accommodations arise not because of the inapplicability of chemical
kinetics but because assumptions, that are necessary to allow population average variables to be used,
prevent the straightforward application of fundamental chemical kinetic principles t
o model
construction.

1.

Reaction energy conventions for variable
-
dependent kinetic coefficients

It is assumed that all variable
-
dependencies of the kinetic coefficients arise because the variable
in question has an effect on the generalized reaction energy p
rofile in the transition from say state


to
state


as shown in Fig. 6. Let state


have an associated Helmholtz free energy of
B


and state


have
a lower Helmholtz free energy,
B

. Thus, there is a tendency for


states to make a transition to


sta
tes. However, a reactant, undergoing the transformation from


to


, follows the reaction coordinate
and, at every point along the reaction coordinate it is required that the reactant have a certain amount of

11

energy given by the reaction profile. This p
rofile
shows that there is an energy barrier that needs to be
overcome for the reaction to proceed. The reaction
coordinate associated with the barrier energy is the
transition state with a free energy of
B


+
B

. The
activation energy for an


to


t
ransition is the
difference between the free energy of the


state and
the barrier energy imposed by the free energy of the
transition state, i.e.,
B

. The reverse


to


transition requires the larger activation energy,
B


=
B


-

B


+
B


. The ra
te of reactions in both
direction obey Boltzman statistics.

Take the rate of


to


transition to be
regulated by a forward coefficient
k


and the rate of


to


transition to be regulated by a backward
coefficient
k


. According to Boltzman statistic
s,


(26)


(27)

where the pre
-
exponential terms,
k
a

and
k
b
, are attempt frequencies;

and

are the respective
activation energies;


is the Boltzman constant; and
T

is the absolute temperature. The exponential
term,
, expresses the probability that an attempt to make a transition will be successful. The
higher the activation energy (i.e.,
B
ij
), the smaller the probability of success.

It is im
possible to know how factors that change the activation energy impact the reaction
profile. Thus, operational assumptions must be made. One employed here is that multiple factors,

k
,
are assumed to have independent effects on the activation energy, i.e.,


(28)

where

is a reference value independent of the effect of these multiple factors.

Having adopted the convention of Eq. (28) for

calculation of the rate coefficient for a forward
reaction, the rate coefficient of the reverse reaction in the presence of factors changing the reaction
profile could be calculated, in principle, from the equilibrium relation defined by


(29)

For reasons that will be rationalized in several instances below, we will not always follow this principle.


2.

Neighbor
-
dependent kinetic coefficients

Cooperative behavior is an integral part the myofilament contractile process [Brandt et.al., 1
984;
1987; Bremel and Weber, 1972; Murray and Weber, 1980; T.L. Hill, 1985; Tobaccman, 1996; Squire
and Morris, 1998; Lehrer, 1994; McGillop and Geeves, 1993; Moss, 1992; Solaro and Rarick, 1998].

12

Most likely, this cooperativity is the result of interacti
ons between neighbor locations along the
filaments. Although there are many possibilities, we adopt the approach of Razumova et. al. [2000], and
include three kinds of cooperative neighbor interactions in the current model: 1) RU
-
RU interactions; 2)
XB
-
XB

interactions; and 3) XB
-
RU interaction (Fig. 7). Because they influence the activation energy
for selected state transitions, the application of our conventions results in neighbor interactions having
complex effects on specific kinetic coefficients. As

an example, we repeat the development for RU
-
RU
neighbor interaction as given in Razumova et. al. [2000] and refer the reader to that reference for
detailed development of XB
-
XB and XB
-
RU interactions formulations.



Figure 7


Three di
fferent kinds of nearest neighbor interactions that lead to cooperative behavior in the myofilament system.


a)

RU
-
RU Interaction

Regulatory units are aligned head to tail along a thin filament. Any unit, whether in the ‘off’ or
‘on’ position, may have 4 pos
sible nearest neighbor configurations as inferred from Fig. 7: 1) both
neighbors ‘off’, (= 11); 2) left neighbor ‘off’ and right neighbor ‘on’, (= 12); 3) left neighbor ‘on’ and
right neighbor ‘off’, (= 21); 4) both neighbors ‘on’, (= 22).

Turning ‘on’

En
d
-
to
-
end interactions between adjacent RU (perhaps through some mechanical coupling due to
overlapping ends of tropomyosin and/or overlapping Tn
-
T) result in the state of the neighboring unit
influencing the propensity of an ‘
off
’ unit to make a transition

to the ‘
on
’ state. Let this influence be
exerted through the activation energy. Thus, for an RU in the ‘
off’

state, the activation energies
associated with a state transition to the ‘
on
’ state may be ordered:


(30)

where the sup
erscripts refer to the states of the left and right neighbors and the ordering is the result of
the influence of the left and right neighbor states on the activation energy required to make the
transition. Thus, the success frequency (i.e., the rate const
ant for transition) from an ‘
off’
to an ‘
on’

state
is greater as more neighbors assume the ‘
on’

state because the activation energy that must be overcome
for this transition declines as more neighbors turn ‘
on
’.

We may address this quantitatively by consid
ering, for the whole population,


(31)


13

where, in accord with the Boltzman expressions above,
k
a

is an attempt frequency and the term in pointy
brackets is the average probability over the whole population that an attempt will be suc
cessful. This
average probability was estimated by taking the sum of weighted probabilities as follows. Weights were
assigned using the assumption that events were randomly distributed along the length of the
myofilament. (This is a far
-
reaching assumpt
ion with important consequences that are examined in
detail in the Critique of Model Section.) With random distribution, the likelihood that a neighboring site
will be in the ‘
off
’ state is

and the likelihood that it will be in the


on
’ state is
. Joint likelihood
are given by the appropriate products as, for example, the likelihood that both neighboring sites will be

off
’ is
, that the right one will be ‘
on
’ and the left one ‘
off
’ is
, etc. Therefore, Eq (31)
may be written in quantitative terms as follows:


(32)

Extracting

out of the expression in brackets gives


(33)

Interaction between adjac
ent RUs impact the activation energy differences

and
. Let this interaction be such that it reduces the activation energy required for an ‘
off
’ to ‘
on

transition by an amount
U
. Then,


(34)

Substituting these into Eq (33), yields


(35)

or


(36)

where

is a reference
k
on

coefficient for the condition where both neighbors are ‘
off
’.
Note, because the effects of Ca
2+

on

k
on

were taken independent of the effects of neighbor interactions,

incorporates the Ca
2+

effect as given by Eq. (5).


14

Because

is simply a number, it may be given the value

such that i
f there is no
effect from neighbor interaction,
U

= 0 and
u

= 1. Further, because

on

+

off

= 1, Eq. (36) can be
rewritten as


(37)

The term in square brackets on the rhs indicates that if there are neighbor interactions (i.e.,
u

>
1), the value of
k
on

increases with increasing number of regulatory units in the “
on
” position.

Turning ‘off’

A similar analysis may be done to determine the average rate constant for the ‘
on
’ to ‘
off

transition. The relations are slightly reordered
for the reverse transition such that


(38)

and, finally,


(39)

where

is the reference
k
off

coefficient for the condition where both neighbors are ‘
off
’.
Note, because the effects of Ca
2
+

on
k
off

were taken independent of the effects of neighbor interactions,

incorporates the Ca
2+

effect as given by Eq. (6).

Interaction between neighbors during the ‘
on
’ to ‘
off
’ transition have the effect that, if most RU
are ‘
off
’ (i.e.,

on

is small), the effect of
u

is to increase
k
off
. One the other hand, if most RU are ‘
on

(i.e.,

on

is large) the effect of
u

is to decrease
k
off
.

Net Steady
-
state Effect

The net steady
-
state effect of
u

on the transitions between ‘
off
’ and

on
’ states can be appreciated
from examining the ratio


(40)

If

on

is small (i.e.,

off

is large) as during low Ca
2+
, increasing
u

decreases the
k
on
/k
off

ratio and
RU tend to be held in the ‘
off
’ position as Ca
2+

increases. Howeve
r if

on

is large as during high Ca
2+
,
increasing
u

increases the
k
on
/k
off

ratio and RU tend to be held in the ‘
on
’ position as Ca
2+

decreases.

b)

RU
-
RU, XB
-
RU, and XB
-
XB Interactions

Let
v

be a parameter that grades the strength of XB
-
XB interactions, and
w

be a parameter that
grades the strength of XB
-
RU interactions. Similar to
u
, the
v

and
w

parameters are such that if they
have the value 1, there is no interaction. Values progressively greater than 1 express progressively
stronger interaction. Utiliz
ing the convention for independent effects on activation energy, analytic
developments in a vein similar to that presented for RU
-
RU interaction lead to the following expressions
for the impact of neighbor interactions on the respective kinetic coefficient
s:


15


(41)



(42)


(43)



(44)


If there are no neighbor interactions, i.e.,
u

=
v

=
w

= 1, then
;
;
; and
.

3.

Distortion
-
dependent kinetic coefficients

A central tenet in the XB theory of muscle contraction is that rate coefficients governing state
transitions depend on the mechanical energy (strain
energy) imparted to the XB by

shearing motions between
thick and thin filaments [Dijkstra et.al., 1973; Homsher et.
al., 1997; Huxley, 1957; McMahon, 1984; Eisenberg et.al.,
1980; Pate and Cooke, 1989]. Distortions occur in
attached states and influence the value of kinetic
coefficie
nts governing reaction paths leading away from
the respective attached state. Energies associated with
distortion are the elastic energies resultant from stretching
the XB spring. These add to the energy of state as shown
in Fig. 8.

a)

Effects of Distortio
n on Energy
of XB States

The XB is assumed to be linearly elastic and thus, the force
due to distortion is given by
. Consequently, the
mechanical energy associated with XB distortion,
. This mechanical energy

adds to the
chemical free energy of a given state to determine its total
free energy. For instance, if we consider
A
1

and
A
2

as the
two XB states that may undergo distortion, we may take
the free energy of these states to depend on distortion
according t
o Fig. 8:


16




(45)



(46)

The first term on the rhs in Eqs. (45


46) is the free energy associated with the non
-
distorted condition
and the 2
nd

term is the mechanical energy component.

Given the above, the mech
anical work done in an
A
1

to
A
2

transition is given as





(47)

and, for the reverse transition,




(48).


b)

Effect of Distortion on Power Stroke Coefficients


h and h’

To approximate the effects of distortio
n on the power stroke coefficients, we perform a thought
experiment in which we consider two extremes of filament sliding during the power stroke (Fig. 9):
extreme 1, no filament sliding as in an isometric contraction; and extreme 2, free filament sliding
with
no load on the thin
-
filament. In the first extreme, head rotation during the forward power stroke
stretches the molecular spring and stores mechanical energy within the tail segment of the XB to bring
about chemical to mechanical energy transduction,

whereas head rotation in the reverse power stroke
relieves the stretch in the molecular spring and affects, one assumes, mechanical to chemical energy
transduction. In the 2
nd

extreme, the thin filament is free to slide past the thick filament and head
r
otation moves the thin filament with no stretch of the spring and no energy transduction in either the
forward or reverse power stroke (this assumes that the thin filament has no mass and its motion does not
give it any kinetic energy).



17

The critical point in our argument hinges on the notion
that the speed of the power stroke under
extreme 1 (isometric), when there is potential energy transduction, is slower than the speed of the power
stroke under extreme 2 (no load) when t
here is no
potential energy transduction. If that notion is
accepted, then one must conclude that the
activation energy to be overcome for
A
1
-
2

state
transition in extreme 1 must be greater than that
for extreme 2 (Fig. 10). Because the free energy
of th
e
A
1

state is identical in both extremes, these
two extremes must represent two different barrier
energy conditions. Thus, it is reasonable to
assume that the barrier energy is affected by the
magnitude of the energy transduction (i.e., the
work) that is
performed. It follows that




(49)

where

is the lowest value of the barrier
energy during no
-
load sliding and
W
1
-
2

is the
work done in the
A
1
-
2

transition (and is also a
measure of the energy transduced fr
om chemical
to mechanical form), and


is a constant of
proportionality that grades the effect of work to
impact the barrier energy.



c)

General case of power stroke during ongoing filament sliding

In the most general case, the power stroke occurs in the pre
sence of ongoing thin
-
filament sliding
against a load as when the sarcomere is shortening. In accord with Eqs. (20
-
21), during shortening,
A
1

XB are compressed and
x
1

< 0 whereas
A
2

XB are distorted less than during an isometric contraction and
x
2

<
x
0
.

T
hen, in accord with our thought experiment, we let the activation energy of the forward power
stroke depend on the work that must be done on the molecule to make that transition in the following
manner:


(50)


(51)


(52)


18


(53)


(54)

Substituting the activation energy expression into the Boltzman expression for the rate coefficient


(55)

where the pre
-
exponential term,
h
nL
, represents the maximal possible value for
h

when there is no energy
transduction. Because under most conditions that may be encountered
x
2

>
x
1
, the exponential term
represents a factor that causes
h

to be <
h
nL

as a result of energy transduction.

For
the reverse direction, the transduction event is a mechanical to chemical transduction. We
assume that the effect on the barrier energy is independent of whether the transduction event is chemo
-
mechanical or mechano
-
chemical. In a manner similar to the a
bove, we find that


(56)

where analogous interpretations are given to the pre
-
exponential and exponential terms except that,
because
x
2

>
x
1

generally, the exponential term represents a factor that causes
h’
>
h’
nL

.

Some feeling for

these distortion effects are obtained from the following analysis.
Because
,

T

= 3.77 x 10
-
21

J

@ 0º or 4.14 x 10
-
21

J

@ 37º. Further, because


= 4
x 10
-
4

Nm
-
1

,
. Now, suppose
x
0

= 10
-
8

m
. Then,
x
0
2

=

10
-
16

m
2

and
.
Therefore, during isometric conditions,
x
1

= 0 and
x
2

=
x
0
. Thus,

and
. Therefore, the isometric values of each of these coefficients are respectively
smaller and la
rger than their free
-
sliding, no
-
load values and the magnitude of difference between
isometric and no
-
load values depends on the strength of the distortion effect represented by the
parameter

.

Care must be exercised in applying these interpretations to the unloaded myofilament system of
the sarcomere. Because of the uni
-
directional detachment step, an unloaded sarcomere will shorten
when XB are cycling. This directionality to unloaded fil
ament sliding, means that individual XBs will
be loaded; XB in the
A
1

state will be compressed (
x
1

< 0) and XB in the
A
2

state will be extended (
x
2

>
0). With no external load, the compression force from
A
1

XB will just equal the tension force from
A
2

XB.

In this situation, we can safely predict that the effect of shortening, in the presence of distortion
effects (


> 0), will be to increase
h

above its isometric value. However, it is not so clear what effect
this will have on
h’
.



19

d)

Effect of Distortion on Detachment Coefficient, g

The effect of distortion on the detachment coefficient,
g
, is more problematic. Foll
owing the
activation energy argument from above, we can write


(57)

where

is the activation energy for an
A
2

to
D

transition when the
A
2

state is unstrained,

is
the elastic energy impa
rted to the
A
2

state as a result of elastic,
x


distortion, and

(x
2
)

is an unknown
functional effect of
x
2

distortion on the reaction profile. Thus,


(58)

where
g
0

is the coefficient regulating XB detachment when
A
2

is in the unstrained condition,


is half
the value of the sp
ring constant in

T

units, and

(x
2
)

is an unknown function.

Complications in applying Eq. (58) to population
-
average differential equations, as we are
attempting to develop here, arise when:
1
) it is necessary to consider non
-
random distributions of
x
2

wi
thin the populations; and
2
) when
g
0

becomes a function of
x
2

for reasons unrelated to energy
considerations. An example of the latter is the manner in which Huxley assigned
g1

and
g2

values as a
function of
x

in his classic 1957 model [Huxley, 1957]. Ap
proaches for dealing with these
complications are given in the Critique of Model section.

4.

Length
-
dependent kinetic coefficients


Length
-
dependent Ca
2+

activation is a well
-
known property of both cardiac and skeletal muscle.
There now seems to be a consens
us that this length dependence is secondary to changes in filament
lattice spacing, whereby increasing SL causes compression of the filament lattice which reduces the
energy required for and facilitates XB attachment [Godt and Maughan, 1981]. In our model
, such an SL
effect increases the probability of a
D

to
A
1

transition. To consider the impact of SL on
D

to
A
1

transitions, we defined a normalized SL variable as



(59)

where
SL
0

is a reference SL. Let decreased lattice spacing due

to increased SL reduce the activation
energy,
B
, needed for a
D

to
A
1

transition as


(60)

where


is a coefficient expressing the strength of the activation energy reduction and the 0 superscript
on
B

refers the first term on the rh
s to the activation energy at
SL = SL
0
. Then, according to our
conventions, the rate coefficient for a
D

to
A
1

transition is enhanced with


according to


(61)

where
f
a

is an attempt frequency. Coincidentally, following the convent
ion given by Eq. (28), the
probability for an
A
1

to
D

transition is reduced according to


20


(62)

The outcome from this treatment of the effect of
SL

on the kinetics of the attachment step is that the
strength of these effects in the m
odel can be graded with the single parameter,

.

F.

Activator Ca
2+

Dynamics

Rather than attempt to model the several aspects of activator Ca
2+

dynamics which would
involve at least as much complexity as already considered for the myofilament system [Ashley et. al.,
1991], the purposes of this modeling
effort are suitably met by treating Ca
2+

as a time
-
varying input
function,
Ca(t)
, acting to drive the myofilament system. This input driving function may be modified as
appropriate for a wide variety of model applications including use of the model to stu
dy the twitch
where
Ca(t)

may be given a pulse
-
like character as occurs following a single excitation event or use of
the model to study mechanodynamics during constant activation in skinned fibers where
Ca(t)

would
simply be made a constant. In addition,

varying time courses may be ascribed to
Ca(t)

in accord with
the situation to be studied or emulated.


G.

Model Summary and Reduced Form for Isometric Conditions

To summarize, the model consists of five differential equations: three (Eqs. (7
-
9)) describe
dynamics of XB states
D
,
A
1
, and
A
2

and two (Eqs. (20, 21)) describe dynamics of XB distortion (
x
1

and
x
2
)
. In addition, there are nonlinear relationships, including: the relation of
Ov

to
SL

(Eqs. (1
-
3)); the
saturation effects due to finite numbers of X
B within
Ov

(Eq. (4)) and finite number of binding sites for
Ca
2+

on TnC (Eq. (5
-
6); the dependence of rate coefficients on distortion (Eqs. (55, 56, 62)), neighbor
interactions (Eqs (41
-
44)), and
SL

(Eqs. (47, 48)); and the appearance of products of state

variables in
the stiffness
-
distortion expression for force (Eq. (22)). In its full nonlinear form, the model contains 22
parameters: 16 associated with kinetic coefficients (
f
0
, f
0
', h, h', g
0
,


Ca
50
,

,
u, v, w,

,

); 5 associated with structural factors (
L
A
, B, L
M
, R
T
, x
0

); and one linking structural and kinetic
factors (

).

Isometric contraction allows certain model simplifications. Specifically, for isometric twitches,
and, according to Eqs. (20, 21), there are no variations in average distortions. Thus, throughout
the isometric contraction episode,


(51)


(52)

As a result, according to Eq. (22),


(
53)

Further, by using the isometric condition as a reference, distortion dependence of kinetic coefficients is
never exhibited and, lacking cooperative effects on the power stroke and detachment,
h
,
h’
, and
g

are
constants equal to their isometric value.

A
s a consequence during isometric contraction, dynamic dimensions of the model are reduced
from 5 to 3 and the nonlinear behaviors associated with distortion dependence never become expressed.


III.

Model Critique

Like all models, the current model is a simplifi
cation. For instance, we considered RUs and XBs
interacting only along a single thin filament and ignored interactions that may occur among multiple
Figure A.2

AAAAA.1


21

thick and thin filaments. A partial list of factors that were ignored include: 1) end
-
effects (we assumed

an infinitely long thin filament); 2) restrictions due to spacing between actin attachment sites along thin
filaments and XBs spacing along thick filaments; 3) compliant properties of thick and thin filaments; 4)
steric relations that allow one RU, in spa
nning 7 actin monomeres, to regulate the availability of as
many as 3 XB attachment sites on the thin filament; 5) more than one Ca
2+

regulatory binding site on an
RU; 6) states of an RU other than ‘
on
’ or ‘
off
’; 7) the multiple states within the XB cycle;

8) influence of
the titin filament to transmit forces directly to the thick filament, to interact with the thin filament in the
I
-
band, and to exert regulatory effects of thick
-
thin filament interactions ; and 9) many other
complicating factors that are k
nown to exist within the myofilament system and the regulated interaction
between actin and myosin. All assumptions were made to allow focus on specific phenomena and in
order to keep the problem tractable.

Some assumptions have more far reaching conseq
uences than others and these are discussed
below.

A.

Number of States

1.

Cross Bridge Cycle

The question of how many states need to be represented in a model of the XB cycle has been
vigorously debated since the publication of A.F. Huxley’s 1957, 2
-
state XB mode
l. Fundamentally,
there are two questions: 1) how many states need to be represented to account for all the distinct
biochemical steps in the hydrolysis of ATP and 2) how many states can be observed from mechanical
perturbation studies. From the biochemi
cal side, the issues concern ATP binding to myosin, ATP
hydrolysis, and product release from the adenine nucleotide binding site all as related to concurrent actin
binding and chemo
-
mechanical energy transduction. Most kinetic schemes have been elaboratio
ns of
the basic 4
-
state scheme originally proposed by Lymn and Taylor [1971] in which both hydrolysis
products are released during a single step in the attached state followed by ATP binding and then
detachment. However, it is now known that the two hydro
lysis products, P
i

and ADP, are released
sequentially, P
i

is released following actin binding and ADP is released following S1 conformational
change. ADP must be released forming the non
-
nucleotide rigor state before ATP may bind. After ATP
binding, the
S1 myosin head detaches from actin and hydrolysis takes place while myosin is detached
from actin. When all biochemical steps are accounted for, the number of states associated with actin
-
activated, myosin catalyzed ATP hydrolysis is quite large.

It has b
een a challenge to associate the biochemical states with chemo
-
mechanical events.
Following the work of Huxley and Simmons [1971] suggesting that there was rapid equilibrium between
two attached states to account for rapid mechanical transients, a 3
-
state

XB cycle model (two attached
states, one detached states) was proposed by Julian et. al. [1973] and then, a 4
-
state model (two attached
states, two detached states) by Eisenberg et. al. [1980]. In a comprehensive evaluation of nucleotide
ligand binding a
nd contractile system mechanics, Pate and Cooke [1989] made use of a 5
-
state model
(three attached states, two detached states) to interpret the effect of varying ATP, ADP, and P
i

concentrations on force and velocity and the force
-
velocity relationships.
Recently, Thomas and
Thornhill [1996] proposed a 6
-
state mechanodynamic model.

However, for fixed concentrations of nucleotide ligand, Kawai and coworkers [Kawai and
Brandt, 1980; Murphy et. al., 1996; Wang and Kawai, 1996; Kawai et. al., 1993; Kawai, et.

al., 1990;
Zhao and Kawai, 1994; Kawai and Halvorson, 1989] have shown repeatedly that no more than three
dynamic processes may be identified from mechanical perturbation experiments in muscles receiving a
variety of treatments. These findings have been
substantially confirmed and extended by Maughan and
coworkers [Maughan et. al., 1998; Blanchard et. al., 1999; Moore et. al., 1999; Dickinson et. al., 1997]
who investigated a wide range of muscles that had undergone genetic modification to express differe
nt

22

isoforms of contractile proteins. Therefore when there is no variation in nucleotide, no more than three
mechanically relevant states are warranted in representing the mechanodynamics of the XB cycle.
Indeed, variants of the 3
-
state model of Julian e
t. al. [1973] are widely used to describe and interpret a
broad range of observations involving the mechanodynamics of myofilament systems [White, 1973;
Wang et. al., 1993; Cuda et. al., 1997; Dantzig et. al., 1992; Thorson and White, 1983; Thomas and
Thor
nhill, 1995; Murase et. al., 1986; Herzig, 1977; Blanchard et. al., 1999]. Significantly, Pate and
colleagues now make more use of a reduced 3
-
state XB model than they do of their original 5
-
state
model to interpret experimental mechanodynamic data [Wang
et. al., 1993; Cuda et. al., 1997].

Structurally, our 3
-
state XB model is not significantly different from most others that have been
proposed. However, different behaviors from 3
-
state models arise from different values assigned to the
various kinetic
coefficients and from different nonlinear dependencies of these coefficients on variables.
In this respect, our XB model may be somewhat unique in that we used kinetic coefficient values that
rendered the pre
-
power stroke attached state,
A
1
, neither stron
gly bound (as in the Huxley
-
Simmons
model [1971]) nor weakly bound (in the sense of Brenner’s weak bridges [1988]). Rather, with these
coefficient values,
A
1

dynamics contribute little to myofilament system dynamics at low physiologic
frequencies, contrib
ute small effects at mid
-
frequencies, and contribute modest effects at frequencies in
the high end of the physiologic range [Razumova et. al., 1999]. Given the immense amount of
investigation addressing the issue of minimum representation of number of XB

states, it seems unlikely
that more than 3
-
states will ever be required to characterize mechanodynamics under conditions of fixed
concentrations of nucleotide ligand. However, future experimental work will be required to substantiate
the appropriate dist
ribution of kinetic coefficient values within the 3
-
state cycle for various muscle
types.

2.

Regulatory Units

Further, it must be pointed out that treating the RU as either '
on
' or '
off
' does not incorporate the
full range of subtleties that would come from

direct interactions between RU and XB at a single binding
site as implied in the 3 RU
-
state model ( '
blocked
', '
closed
', and '
open
') of Geeves and coworkers
[McGillop and Geeves, 1993; Leher, 1994; Geeves and Leher, 1994] A rough equivalence between the
3 RU
-
state scheme and the one used here is as follows: '
blocked
' would be equivalent to our
R
off
, '
closed
'
would be equivalent to an equilibrium combination of
D

and
A
1
, and '
open
' would be equivalent to an
equilibrium combination of
D

and
A
2
. The transit
ion between '
blocked
' and '
closed
' would be associated
with Ca
2+

binding and the transition between '
closed
' and '
open
' would be associated with the XB power
stroke. In this scheme, cooperativity occurs within the span of a single RU as a result of a sing
le
closed

to
open

transition from the power stroke of one XB facilitating the attachment of additional XBs within
the thin filament span of the RU.

Thus, the kind of cooperativity whereby one XB holds an RU in an open state allowing other
XB’s to readily

attach and proceed through a power stroke within the span of a single RU (Geeves and
Lehrer, 1994, Lehrer, 1994, Tobacman, 1996) is somewhat different than the XB
-
XB and XB
-
RU
cooperativity we represent here. Further differences would arise between the 3

RU
-
state configuration
and our configuration in that, for instance, RU
-
RU interaction would impact both the
blocked

to
closed

transition and the
closed

to
open

transition rather than just the '
on
' to '
off
' transition. Whether one would
arrive at different

conclusions regarding the relative effects of RU
-
RU, XB
-
XB, and XB
-
RU neighbor
interactions using the 3 RU
-
state model as opposed to the 2 RU
-
state model used in this study is a
question worthy of future investigation.



23

B.

Spatial Averaging


Non uniform
distribution

Another assumption of consequence was that states were assumed to be randomly distributed
along the length of the myofilament allowing the likelihood of finding a neighboring site in any
particular state to be calculated according to its fract
ional occurrence. This assumption (used in
Razumova et. al. 2000 for formulating Eqs. 41
-

44) is at odds with neighbor interactions because such
interactions would tend to cluster states together in non
-
random patterns.
Our assumption is known in
statis
tical physics as the Bragg
-
Williams or mean
-
field approximation [Hill, 1978; Hill, 1985] where it is
often used in solving neighbor interaction problems. An alternative is the Bethe
-
Peierls approximation
[Hill, 1985] or 'quasi
-
chemical method' in which sp
atial independence is not assumed and more exact
solutions are obtained for one
-
dimensional, Ising
-
like problems of the kind we are treating here.
However, in addition to being more complicated mathematically, the Bethe
-
Peierls approximation
requires that
all interactions be of the same strength. In our case, this requirement is not satisfied and,
thus, the extra mathematical complexity is not warranted.

By employing the Bragg
-
Williams approximation, we retain a deterministic structure to the
neighbor inte
raction problem and circumvent the requirement for probabilistic Monte Carlo methods in
its solution. Deterministic model structures are desirable not only because they are less involved
computationally but because they allow more straightforward cause an
d effect explanations. The most
important consequence of the Bragg
-
Williams approximation is that it exaggerates positive cooperativity
[Hill, 1985]. In our case, this exaggeration is of little consequence as the interaction parameters (
u, v,

or
w
) appe
ar as free parameters whose values may be chosen in accord with the behavior that the user
wishes to simulate. That is, the model should be viewed as a tool in which the interaction strengths are
assigned to represent system behavioral consequences and no
t as exact physical entities representing
specific molecular interactions. However, through changes in parameter values, specific interactions
may be graded accordingly.

Our assumption of randomly distributed states should be taken in the same light as

Tobacman’s
comment [1996] about all previous cooperative models



The cooperative mechanisms put into
mathematical form by the models imply specific statistical distributions of myosin along the thin
filament. Until methods are available for measuring t
hose distributions, no model can be well
substantiated and all models must be viewed with caution.
” The single relevant caution arising as a
consequence of the Bragg
-
Williams approximation is that care be exercised to avoid large
u
,
v
, or
w

values that pr
oduce phase transition or critical point in the binding isotherm (i.e., force
-
pCa). Phase
transition leads to a type of hysteresis that is uncharacteristic of the hysteresis in published force
-
pCa
data [Harrison et al., 1988; Brandt et. al., 1985]. With
this caveat, our assumption of random
distribution is perfectly in accord with standard practices in solving neighbor interaction problems.


C.

Averaging Variable Dependence on Distortion

Previous treatments of distortion
-
dependent effects on rate coeffici
ents have taken into account
thermal effects causing a distribution of distortions within a population and the impact of strain energy
on state transitions. One consequence of these considerations is that the differential equation describing
distortion
-
de
pendent kinetic events becomes a partial differential equation in the distortional variable.
The resulting mathematical problem takes on a fair measure of difficulty that detracts from the utility of
the model.

There have been various approaches for simp
lifying the problem of distributed distortion among
XBs. As outlined by Taylor et. al. [1993], the simplest and, probably the most often used approach, is
simply to ignore the spatial derivative term and solve the resulting ordinary differential equations

as if
distortion in each XB state were uniform. More direct have been approaches to formally simplify the

24

problem. Perhaps the first serious attempt was by Deshcherevskii [1971]. Deshcherevskii considered
that above one value of distortion, XBs were pu
lling in the direction of shortening while for distortion
below this value, XBs were hindering shortening. The detachment rate constant was then set to be
single valued in each of these regions. This allowed completing a spatial integration over the dis
tortion
coordinate such that the partial derivative terms were reduced to variables and the equations could be
written in ordinary differential equation format. A related approach was used by Dijkstra et. al. [1973]
where the distortion coordinate was dis
cretized and all XBs within a discrete distortion range
x




x

were grouped into a bin. All XBs within a bin possessed one value of rate coefficient but that value
varied between bins according to assigned discrete quantities. Ordinary differential equations were
written for XBs entering and leaving each bin d
uring filament sliding. In accord with the approach of
Deshcherevskii [1971], the number of bins employed by Dijkstra et. al. [1973] could be reduced to as
few as two and reasonable results were obtained. Another approach is that of Zahalak [1981; 1986]
where a distribution function for attached XBs was assumed and integration was performed to determine
the moments of the bond
-
distribution giving rise to net population stiffness, force, and elastic energy.
Ordinary differential equations for these quanti
ties could then be written and integrated to determine
their respective dynamic variation. All these approaches to simplification possess considerable merit but
suffer in some details when more than one attached XB state needs to be considered.

Because th
e basic premise of our model was that muscle force equaled sarcomeric stiffness times
average XB distortion (Eq. 25), our challenge was to find some reasonable relationship between
population average rate coefficient and average distortion. We recognized t
hat, starting with such
elementary expressions as one may want for the relationship between the rate coefficient and a given
x
, a
strict mathematical derivation for relating the average rate coefficient to the average
x

cannot be found.
However, there is
no doubt that induced variation in average distortion from its isometric value would
be associated with variation in the associated average rate coefficient from its isometric value. To
express this, it was convenient to adopt parametrically succinct fun
ctional forms that satisfied the
experimental fact that the numbers of attached XBs decline during constant shortening [Julian and
Sollins, 1975]. This meant that during shortening, rate coefficients leading away from an attached state
must increase relat
ive to those leading into that state. Additional considerations resulted from adopting
parameter values (Table 1) such that XBs in the
A
1

state were short
-
lived relative to crossbridges in the
A
2

state (this is somewhat, but not totally, analogous with th
e commonly used concepts of weak
-
binding/strong
-
binding states). The consequence was that, for a given
,
x
1

was never forced far
from its isometric value of 0 whereas
x
2

was forced relatively far from its isometric value of
x
0
. Thu
s,
those rate coefficients leading away from the
A
1

state were much less affected by distortion than those
leading away from the
A
2

state.

To summarize, the stiffness
-
distortion approach to muscle modeling leads to development of
ODE expressions for avera
ge distortion of attached crossbridge states. This average distortion may then
be used to introduce distortion
-
dependence into population
-
average kinetic rate coefficients. Our
approach requires a phenomenologic approximation because: 1
-

no account is ta
ken of the role of
distortion dependence in the derivation from the macroscopic balances leading to Eqs. (3
-
5, 11, 12) and
these ODEs are global approximations of the PDE’s that properly describe the detailed behavioral
features of the system; and 2
-

lack
ing a strict mathematical relationship between population average
quantities, the underlying thermodynamic constraints on kinetic relationships can not reasonably be
applied. Thus, the distortion

dependent expressions are arbitrarily formulated to achieve

parametric
succinctness. As a result, model behavior, rather than internal consistency, becomes the arbiter of model
acceptability.


25


D.

Ca
2+

kinetics and cooperativity

Another important assumption is that Ca
2+

binding to TnC and its effect on
k
on

and
k
off

were
assumed to be fast and independent of all other effects and of other actions associated with myofilament
activation and cross
-
bridge cycling. There is both experimental evidence and theoretical reason for
taking note of this assumption. Experiment
ally, Hoffmann and Fuchs [1987a, 1987b] have shown in
cardiac muscle that Ca
2+

binds to myofilaments with less affinity at low force (vanadate inhibition),
when there are fewer XB’s in the force
-
bearing state, than at high force when there are more. While

not
definitive, this evidence at least strongly implies that there is a dependence of Ca
2+

binding on force
-
bearing XB state [Fuchs and Wang, 1995]. Theoretically, a detailed balance, based on thermodynamic
equilibrium considerations, demonstrates that o
ur assignments of

and
, while in
accord with basic understanding of the
influence of Ca
2+

on the
on
-
off

transitions,
are, never
-
the
-
less, inconsistent with one
ratio of
k
+

/
k
-

for both ‘
on
’ and ‘
off
’ states.

Indeed, T.L. Hill’s thermodynamically
consistent model of myofilament activation
and cross
-
bridge cycling [Hill, 1985]
specifically accounts for different Ca
2+
-
binding affinities depending on whether the
RU is ‘
on
’ or ‘
off
’. In Hill’s treatment, Ca
2+

eff
ects on equilibrium constants are not
independent of other effects.

However, by assuming that Ca
2+

binding to TnC was fast and independent of
other events, we were able to utilize the
simple kinetic scheme in Fig. 3 rather than
the more elaborate scheme of

Fig. 11 which
shows XB cycling for RU conditions that are
both Ca
2+
-
bound (left side) and non Ca
2+

bound (right side). The more elaborate
scheme of Fig. 11 would require different
values of
k
+

and
k
-

at each left
-
side/right
-
side
transition and would also

suggest different values for kinetic coefficients at equivalent steps in the XB
cycle on left and right sides of the figure. Such additional complication is not warranted at this time.

While respectful of experimental findings and thermodynamic constrain
ts, we employ the Ca
2+
-
binding assumption in our model because: 1) we desired to examine specific neighbor interactions
without the complications and obfuscations of multiple other effects as would arise if Ca
2+

interactions
were included and 2) the mathem
atics of including Ca
2+

interactions may not be tractable in a model
whose ultimate use will be for non
-
equilibrium applications during transient behaviors such as force
development. Intuitively, we suspect that increases in Ca
2+
-

binding affinity with an

increase in XB
force
-
bearing states will, among the three neighbor interactions investigated here, have its greatest effect
on the prediction of XB
-
RU interactions by shifting the force
-
pCa curve further to the left than would
occur with XB
-
RU interaction

alone.
Suppose there is a certain level of XB
-
RU interaction and a
resultant shift to the left of the force
-
pCa curve.


Now, if cooperative Ca
2+

binding is imposed on top of

26

existing XB
-
RU interaction, the result will be to shift the curve further left.


These two effects are in the
same direction and if one is not considered but is operative, the other effect is over estimated.

This
means that if the effects of XB
-
dependent Ca
2+
-
binding are operative and not considered, a larger value
of
w

will be estim
ated than may actually exist as a result of XB
-
RU interaction alone. Future model
improvements will be required to determine the effects of a dependence of Ca
2+
-
binding affinity on XB
state.

E.

Heterogeneity and Extra
-
Sarcomeric Structures

Force prediction
s from a model of the sarcomere may be simply scaled up according to a
muscle’s cross sectional area and percentage area that is occupied by sarcomeres to predict overall
muscle force providing that all sarcomeres constitute a homogenous population and the
re is no extra
-
sarcomeric series and parallel elasticity. Of course, neither requirement is met in a realistic situation.
Even within a muscle fiber there will be a certain amount of parallel and series inhomogeneity as well as
parallel and series elasti
city borne by structures not in the sarcomere. When multiple muscle fibers are
considered even within a single motor unit, inhomogeneity and parallel and series elasticity contribute
even more to overall muscle dynamics. Finally, the angle of connection
of the muscle fibers to the
muscle’s tendinous origin and insertion as well as the mechanical properties of the tendon itself will
have significant impact on the dynamics of the muscle organ. Sarcomeric heterogeneity and parallel and
series elasticity ar
e additional factors requiring future study and incorporation into the model to
complete the muscle modeling effort.


IV.

Implementing the Model

A.

Input
-
Output Representation

Because of its high dynamic dimensionality (5
th

order) and numerous nonlinearities, the

response
profile of the above model is enormous and equal to the full spectrum of experimental observations that
have been reported for muscle. An approach to using the model that is homologous to many
experimental muscle investigations and that also len
ds itself for coupling the model to neural input and
mechanical environments, is to structure the model (Fig. 12) such that it generates a force output in
response to two inputs: 1) a change in Ca
2+

activation,
Ca(t)
, and 2) a change in sarcomere length,
S
L(t)
.


Figure 12


Input
-
output representation of sarcomeric myofilament system that recieves two inputs, activator calcium,
Ca(t)
,
and sarcomere length,
SL(t)
, and responds with a single output, force,
F(t)
.


B.

Summary Equa
tion List

Model equations are presented below in an order to accommodate numerical solution for this
system of equations when viewed in the format of the input
-
output relations of Fig. 12. Equations are
numbered as they appear in the text.

1.

Input
-
dependent

Coefficients

a)

Activator Calcium dependent coefficients


(5)


27


(6)

b)

Sarcomere Length dependent coefficients

(1)

Filament Overlap

If:
SL < 2LA


B,
then




(1)

If:
2LA


B < SL < 2LA + B,

then





(2)

If:
SL > 2LA + B,

then





(3)


(4)


(9)

(2)

Coefficient of length
-
dependent attachment



(59)


(61)


(62)


2.

State variable
-
dependent coefficients

a)

Neighbor Interactions


(41)



(42)


(43)



(44)



28

b)

Distortion dependent coefficients


(5
5)


(56)


(58)


3.

State Variable Differential Equations


(7)


(8)


(9)


(21)


(20)

4.

Output Equation


(22)


C.

Parameter Values

There is no one set of parameters for a general model of striated muscle because of the great
diversity: 1) among muscles adapted for specific functions in different species, 2) among muscles
adapted for different fu
nctions within a single organism, and 3) within single muscles undergoing
temporal adaptations because of disease, changing use patterns, and mechanical environments. The
parameter set given below is only representative and one can expect 2
-
10 fold differ
ences from these
values for any given muscle.

Kinetic Baseline

f
0

(s
-
1
)

f'
0

(s
-
1
)

h
iso

(s
-
1
)

h'
iso

(s
-
1
)

g
0

(s
-
1
)


(s
-
1
)


(s
-
1
)


(s
-
1
)


(s
-
1
)

Ca
50

(ML
-
1
)

150

50
0

50

50

25

150

50

0

150

1.78 x
10
-
6












Null value for nonlinear coefficients



u

v

w





0

1

1

1

0

0








29


Structural

L
A

B

L
M

R
T

Slack

x
0

SL
0

1.2


0.2


1.6


1.62 x 10
5

1.8


10
-
2


2.4










Mechanical



3 x 10
-
10

N

-
1



D.

Soluti
on Procedures

1.

Initial Conditions

Initial conditions for state variables either to begin a numerical integration procedure or to
represent a given steady state around which behavior is to be examined depend entirely on the situation
that the model user wish
es to emulate. The only situation in which initial conditions can be specified
a
priori

is the condition of zero Ca
2+

activation. Here, there are no cycling XB and the initial conditions
are:
D

=
A
1

=
A
2

= 0 ,
x
1

= 0, and
x
2

=
x
0
. For all other situati
ons, the user has two options for finding
initial conditions: 1) Enter the initial conditions of the input variables, set all the derivatives in the state
-
variable differential equations to zero, and numerically solve the resulting simultaneous set of nonl
inear
algebraic equations for steady
-
state values of state variables; alternatively, 2) Set preliminary initial
conditions to the zero Ca
2+

activation situation, enter the initial conditions of the
Ca(t)

and
SL(t)

input
variables, numerically integrate the

differential equations until the state variable derivatives approach
zero, choose these steady
-
state state variable values as the true initial conditions, then allow the
Ca(t)

and
SL(t)

input variables to change according to the situation to be emulated a
nd numerically integrate to
obtain state variable and output variable predictions.

Our experience is that the second of these two approaches is the more satisfactory from a
programming stand point.

2.

Input Variable values

Having chosen
Ca
50

to be 1.78 x 10
-
6

ML
-
1
,
Ca(t)

may take on values typically found in the cell
or those used during constant activation experiments, i.e., 10
-
8

to 10
-
4

ML
-
1
. For the single low
-
affinity
binding isotherm with no cooperativity as given in Eq. (5 and 6), it may be more conve
nient to rearrange
that equation and express
Ca(t)

in terms of multiples of
Ca
50
. The time function for
Ca(t)

will vary
depending on the application and the muscle as discussed in the Critique of Model section.

Because there is currently no internal loadi
ng forces in the sarcomere other than compression of
pre
-
power stroke,
A
1
, state,
SL(t)

should not be given values much below slack length (1.8

) and
certainly, not below thick filament length of 1.6

. Values of
SL(t)

up to 2
L
A

+
L
M

(3.8

) where
filame
nts no longer overlap are allowed. However, because there are no passive mechanical elements in
the current model, it must be appreciated that predicted force is strictly that due to active mechanisms.
At long
SL
, total muscle force is dominated by force

arising from passive structures.

3.

Numerical Integration

Any standard numerical integration package will provide a solution for the model differential
equations. We commonly use a fixed
-
step size, 4
th

order Runge
-
Kutta routine and find satisfactory
results

with step size of 0.0001 s. However, attention must be given and care exercised when ever the
coefficients of the differential equations are nonlinear as are these coefficients. Highly nonlinear

30

conditions may drive these coefficients to very high value
s and cause the numerical integration
procedure to become unstable.

4.

Output Variable values

F(t)

output from Eq. (22) is in terms of force/sarcomere and, under isometric conditions, can
never be larger than
x
o

R
T

or 4.86 x 10
-
7

N/sarcomere. A representative sarcomere may have a cross
sectional area of 3 X 10
-
12

m
2
. Thus,
F(t)

from Eq. (22) may be converted to cross sectional area of
muscle as appropriate by taking into account only that fraction of muscle cro
ss sectional area that is
myofilaments.




31

References


Ashley, C.C., I.P. Mulligan, T.J. Lea. Ca
2+

and activation mechanisms is skeletal muscle.
Quarterly
Review of Biophysics

24: 1
-
73, 1991.

Blanchard, E., C. Seidman, J.G. Seidman, M. LeW
inter, D. Maughan. Altered crossbridge kinetics in
the

MHC
403/+

mouse model of familial hypertrophic cardiomyopathy.
Circ. Res.

84:475
-
483,
1999.

Brandt, P.W., B. Gluck, M. Mini, C. Cerri. Hysteresis of the mammalian pCa/tension relation is small
and muscle specific.
J Musc Res Cell Motil

6:197
-
205, 1985.

Brandt, P.W., M.S. Diamond, F.H. Sachat. The thin filament of vertebrate skeletal muscle cooperatively
activates as a unit.
J. Mol. Biol.

180:379
-
384, 1984.

Brandt, P.W., M.S. Diamond, J.S. Rutchik, F.H. Sachat. Co
-
operative interactions between troponin
-
tropomyosin units extend the length of the thin filament in skeletal muscle.
J. Mol. Biol.

195:
885
-
896, 1987.

Bremel, R.D., and A. Weber. Cooperation within actin filament in vertebrate skeletal muscle.

Nature
New Biol.
238: 97
-
101, 1972.

Brenner, B.
Effect of Ca2+ on cross
-
bridge turnover kinetics in skinned single rabbit psoas fibers:
implications for regulation of muscle contraction.
Proc. Natl. Acad. Sci. USA

85: 3265
-
3269,
1988.

Brown, I.E. S. H. Scott, G. E. Loeb. Mechanics of feline soleus: II

design and validation of a
mathematical model.
J Musc Res & Cell Motil

17: 221
-
233, 1996.

Campbell, K. B., H. Taheri, R. D. Kirkpatrick, T. Burton, and W. C. Hunter. Similarities between
dynamic elastance of left ventricular chamber and papillary muscle o
f rabbit heart.
Am. J.
Physiol.

264: H1926
-
H1941, 1993.

Campbell, K. Rate constant of muscle force redevelopment reflects cooperative activation as well as
cross
-
bridge kinetics.
Biophys.J.
72: 254
-
262, 1997.

Cuda, G, E. Pate, R. Cooke, and J. R. Sellers.

In vitro actin filament sliding velocities produced by
mixtures of different types of myosin.
Biophys. J.
72: 1767
-
1779, 1997.

Curtin, N.A., A.R. Gardner
-
Medwin, R.C. Woledge. Predictions of the time course of force and power
output by dogfish white mus
cle fibers during brief tetani.
J. Exp. Biol.

201: 103
-
114, 1998.

Daniels, T.L., A.C. Trimble, P.B. Chase. Compliant realignment of binding sites in muscle: Transient
behavior and mechanical tuning.
Biophys. J.

74: 1611
-
1621, 1998.

Dantzig, J. A., Y. E.
Goldman, J. Lacktis, N.C. Millar, and E. Homsher. Reversal of the cross
-
bridge
force generating transition by photogeneration of phosphate in rabbit psoas muscle fibers.
J.
Physiol. (Lond.)

451: 247
-
278, 1992.

Deshcherevski, V.I. A kinetic theory of stri
ated muscle contraction.
Biorheology
, 7:147
-
170, 1971.

Dickinson, M.H., C.J. Hyatt, F
-
O. Lehman, J.R. Moore, M.C. Reedy, A. Simcox, R. Tohtong, J.O.
Vigoreaux, H. Yamashita, D.W. Maughan. Phosphorylation
-
dependent power output of
transgenic flies: An inte
grated study.
Biophys. J.

73: 3122
-
3124, 1997.

Dijkstra, S., J.J. Denier van der Gon, T. Blange, J.M. Karemaker, A. E. J. L. Kramer. A simplified
sliding
-
filament muscle model for simulation purposes.
Kybernetic

12: 94
-
101, 1973.

Dorgan, S./J., M.J. O'Ma
lley. A nonlinear mathematical model of electrically stimulated skeletal
muscle.
IEEE Trans Rehab Eng

5:179
-
194, 1997.

Eisenberg, E., T.L. Hill, Y. Chen. Cross
-
bridge model of muscle contraction.
Biophys. J.

29: 195
-
226,
1980.


32

Finer, J.T., R.M. Simmon
s, J.A. Spudich. Single myosin molecule mechanics: piconewton forces and
nonometre steps.
Nature

368: 113
-
118, 1994.

Fuchs, F and Y.P. Wang. Sarcomere length versus interfilament spacing as determinants of cardiac
myofilament Ca
2+

sensitivity and Ca
2+

binding.
J. Mol. Cell Cardiol.

28: 1375
-
1383, 1996.

Geeves, M.A., S.S. Lehrer. Dynamics of the muscle thin filament regulatory switch: the size of the
cooperative unit.
Biophys. J.

67:273
-
282, 1994.

Godt, R.E., D.W. Maughan. Influence of osmotic compre
ssion on calcium activation and tension in
skinned muscle fibers of the rabbit.
Pfluegers Arch.

391: 334
-
337, 1981.

Harrison, S.M., C. Lamont, D.J. Miller. Hysteresis and the length dependence of calcium sensitivity in
chemically skinned rat cardiac musc
le.
J Physiol

401:115
-
143, 1988.

Herzig, J.W. A model of stretch activation based on stiffness measurements in glycerol extracted insect
fibrillar muscle. In:
Insect Flight Muscle
. R.T. Tregear, editor. North
-
Holland, Amsterdam.
209
-
219, 1977.

Hill, A.

V. The heat of shortening and the dynamic constant of muscle.
Proc. Roy. Soc.
126B: 136
-
195,
1938.

Hill, T.L. Cooperativity Theory in Biochemistry. Springer
-
Verlag, New York, 1985.

Hill, T.L. and L. Stein. Critical behavior of two
-
state, steady
-
state Is
ing system, according to the Bragg
-
Williams approximation.
J. Chem. Phys.
69: 1139
-
1150, 1978.

Hofmann, P.A. and F. Fuchs. Effect of length and cross
-
bridge attachment on Ca
2+

binding to cardiac
troponin C.
Am J Physiol

253: C90
-
C96, 1987a.

Hofmann, P.A
. and F. Fuchs. Evidence for a force
-
dependent component of calcium binding to cardiac
troponin C.
Am J Physiol

253: C541
-
C546, 1987b.

Homsher, E., J. Lacktis, M.Regnier. Strain
-
dependent modulation of phosphate transients in rabbit
skeletal muscle fibe
rs. Biophys. J. 72:1780
-
1791, 1997.

Huxley, A. F. and R. M. Simmons. Proposed mechanism of force generation in striated muscle.
Nature
233: 533
-
538, 1971.

Huxley, A.F. Muscle structure and theories of contraction.
Prog. Biophys. biophys. Chem.

7: 255
-
31
8,
1957.

Julian, F. J., and M. R. Sollins. Variation of muscle stiffness with force at increasing speeds of
shortening.
J. Gen Physiol.

66:287
-
302, 1975.

Julian, F.J., K.R. Sollins, and M.R. Sollins. A model for muscle contraction in which cross
-
bridge
a
ttachment and force generation are distinct.
Cold Springs Harbor Symp. Quant. Biol
. 37: 685
-
688, 1973.

Kaufman, K.R., K.
-
N. An, W.J. Litchy, E.Y.S. Chao. Physiological prediction of muscle forces


I.
Theoretical formulation.
Neuroscience

40:781
-
792, 1
991.

Kawai, M. and P. Brandt. Sinusoidal analysis: a high resolution method for correlating biochemical
reactions with physiological processes in activated skeletal muscles of rabbit, frog and crayfish.
J. Mus. Res. & Cell Motil.

1: 279
-
304, 1980.

Kawai,
M., H.R. Halvorson. Role of MgATP and MgADP in the cross
-
bridge kinetics inj chemically
skinned rabbit psoas fibers.
Biophys. J.

55: 595
-
603, 1989.

Kawai, M., J.S. Wray, K. Guth. Effect of ionic strength on crossbridge kinetics as studied by sinusoidal
an
alysis, ATP hydrolysis rate and X
-
ray diffraction techniques in chemically skinned rabbit
psoas fibers.
J. Musc. Res. & Cell Motil.

11: 392
-
402, 1990.

Kawai, M., Y. Saeki, Y. Zhao. Crossbridge scheme and the kinetic constants of elementary steps
deduced
from chemically skinned papillary and trabecular muscles of the ferret.
Circ. Res.

73:
35
-
50, 1993.


33

Krylow, A.M., W.Z. Rymer. Role of intrinsic muscle properties in producing smooth movements.
IEEE Trans Biomed Eng

44:165
-
176, 1997.

Lehrer, S. S. The re
gulatory switch of the muscle thin filament Ca
2+

or myosin heads?
J. of Muscle Res.
And Cell Mot.

15: 232
-
236, 1994.

Lymn, R.W. and E. W. Taylor. Mechanisms of adenosine triphosphate hydrolysis by actomyosin.
Biochem.

10: 4617
-
4624, 1971.

Maughan, D., J.
Moore, J. Vigoreaux, B. Barnes, L.A. Mulieri. Work production and work absorption
in muscle strips from vertebrate cardiac and insect flight muscle fibers. In:
Mechanisms of Work
Production and Work Absorption in Muscle
. Edited by Sugi and Pollack. Plen
um Press, New
York, 1998, pp471
-

480.

McKillop, D. F. A. and M. A. Geeves. Regulation of the interaction between actin and myosin
subfragment 1: evidence for three states of the thin filament.
Biophys. J.

65: 693
-
701, 1993.

McMahon, T. A., Muscles, Refle
xes, and Locomotion
.
Princeton University Press, Princeton, New
Jersey, 1984.

Moore, J.R., J.O. Vigoreaux, D.W. Maughan. The Drosophila projectin mutant, bentD, has reduced
stretch activation and altered indirect flight muscle kinetics.
J. Musc. Res. Cel
l Motil.

20: 797
-
806, 1999.

Moss, R.L. Ca
2+

regulation of mechanical properties of striated muscle.
Circ Res

70: 865
-
884, 1992.

Murase, M., H. Tanaka, K. Nishiyama, and H. Shimizu. A three
-
state model for oscillation in muscle:
sinusoidal analysis.
J.

Muscle Res. Cell Motil.

7: 2
-
10, 1986.

Murphy, K.P., Y. Zhao, M. Kawai. Molecular forces involved in force generation during skeletal muscle
contraction
. J. Exp. Biol.

199: 2565
-
2571, 1996.

Murray, J. M. and A. Weber. Cooperativity of the calcium switch

of regulated actomyosin system.
Mol. Cell. Biochem. 35:11
-
15, 1980.

Pate, E. and R. Cooke. A model of crossbridge action: the effects of ATP, ADP, and Pi.
J. musc. Res.
Cell Motil.

10: 181
-
196, 1989.

Razumova, M. V., A. E. Bukatina, K. B. Campbell. St
iffness
-
distortion sarcomere model for muscle
simulation.
J. Appl. Physiol.

87(5): 1861
-
1876, 1999.

Razumova, M. V., A. E. Bukatina, K. B. Campbell. Different myofilament nearest neighbor interactions
have distinctive effects on contractile behavior.
Biop
hys. J.

78: 3120
-
3137, 2000.

Saeki, Y., M. Kawai, and Y. Zhao. Comparison of crossbridge dynamics between intact and skinned
myocardium from ferret right ventricles.
Circ. Res.
68: 772
-
781, 1991.

Schoechting, J.F., M. Flanders. Evaluating an integrated m
usculoskeletal model of the human arm.
J.
Biomech. Eng.

119:93
-
102, 1997.

Schultz, A.B., J. A. Faulkner, V.A. Kadhiresan. A simple hill element
-
nonlinear spring model of muscle
contraction biomechanics.
J. Appl. Physiol.

702:803
-
812, 1991.

Shue G.
-
H., a
nd P. E. Cargo. Muscle
-
tendon model with length history
-
dependent activation
-
velocity
coupling.
Ann Biomed Eng

26:369
-
380,1998.

Solaro, R. J., H.M. Rarick. Troponin and tropomyosin: Proteins that switch on and tune in the activity
of cardiac myofilaments
.
Circ. Res

83:471
-
480, 1998.

Squire, J.M.
The Structural Basis of Muscular Contraction
. Plenum. New York. 1981.

Squire, J. M. and E. P. Morris. A new look at thin filament regulation in vertebrate skeletal muscle.
FASEB J.

12, 761
-
771, 1998.


34

Taylor,

TW, and Y Goto, H Suga. On the solutions of Huxley
-
type models in cardiac muscle fiber
contractions.
J. Theor. Biol.

165:409
-
416, 1993

Thomas, N, and R.A. Thornhill. A theory of tension fluctuation due to muscle cross
-
bridges.
Proc. R.
Soc. Lond. B

259:

235
-
242, 1995.

Thomas, N, and R.A. Thornhill. Stretch activation and nonlinear elasticity of muscle cross
-
bridges.
Biophys. J.

70: 2807
-
2818, 1996.

Thorson, J. and D.C. White. Role of cross
-
bridge distortion in the small
-
signal mechanical dynamics of
in
sect and rabbit striated muscle.
J. Physiol.

343: 59
-
84, 1983.

Tobacman, L.S. Thin filament mediated regulation of cardiac contraction.
Annu Rev Physiol
. 58:447
-
481, 1996.

Trombitas, K. and A. Tigyi
-
Sebes. Cross
-
bridge interaction with oppositely polari
zed actin filaments in
double
-
overlap zones of insect flight muscle.
Nature

309: 168
-
170, 1984.

Trombitas, K. and A. Tigyi
-
Sebes. How actin filament polarity affects crossbridge force in doubly
-
overlapped insect muscle.
J. musc. Res. Cell Motil.

6: 447
-
459, 1985.

van den Bogert, A. J., K. G. M. Gerristen and G. K. Cole. Human muscle modeling from a user’s
perspective.
Proc 9
th

Conf. Canadian Soc. Biomech
. 22, 1996.

Wang, D., E. Pate, R. Cooke, and R. Yount. Synthesis of non
-
nucleotide ATP analogues and

characterization of their chemomechanical interaction with muscle fibers.
J. musc. Res. Cell
Motil.

14: 489
-

497, 1993.

Wang, G., M. Kawai. Effects of MgATP and MgADP on the cross
-
bridge kinetics of rabbit soleus
slow
-
twitch muscle fibers
. Biophys. J.

71: 1450
-
1461, 1996.

Wexler, A.S.,J. Ding, S.A. Binder
-
Macleod. A mathematical model that predicts skeletal muscle force
.
IEEE Trans Biomed Eng

44(5):337
-
348, 1997.

White, D.C.S. Dynamics of contraction in insect flight muscle.
37
th

Cold Springs Harbor
Symp. Quant.
Biol
. 37: 201
-
213, 1973.

Williams, T.L., G. Bowtell, N.A. Curtin. Predicting force generation by lamprey muscle during applied
sinusoidal movement using a simple dynamic model.
J. Exp. Biol.

201: 869
-
875, 1998.

Winters, J. M. An improved mu
scle
-
reflex actuators for use in large
-
scale neuromusculoskeletal models.
Ann. Biomed. Eng.
23: 359
-
374, 1995.

Winters, J.M. and L. Stark. Muscle models: what is gained and what is lost by varying model
complexity.
Biol. Cybern
. 55:403
-
420, 1987.

Wu, J.Z
., W. Herzog. Modelling concentric contraction of muscle using an improved cross
-
bridge
model.
J. Biomechancis

32:837
-
848, 1999.

Zahalak, G. I. A comparison of the mechanical behavior of the cat soleus muscle with a distribution
-
moment model.
J. Biomech
. Eng.

108: 131
-
140, 1986.

Zahalak, G. I. A distribution
-
moment approximation for kinetic theories of muscular contraction.
Mathematical Biosciences
, 55: 89
-
114, 1981.

Zahalak, G. I., An overview of muscle modeling. In:
Neural Prostheses
, R.B. Stein, P.
H. Peckham, and
D. B. Popovich, eds. Oxford Univ Press, New York, pp 17
-
57, 1992.

Zajac F.E. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor
control.
Critical Reviews in Biomedical Engineering
17: 359
-

411, 1989.

Z
hao, Y., M. Kawai. Kinetics and thermodynamic studies of the cross
-
bridge cycle in rabbit psoas
muscle fibers.
Biophys. J.

67: 1655
-
1668, 1994.


35