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
Comments 0
Log in to post a comment