Vol. 87, Issue 5, 1861-1876, November 1999
Stiffness-distortion sarcomere model for muscle simulation
Maria V.
Razumova1,3,
Anna E.
Bukatina1,4, and
Kenneth B.
Campbell1,2
Departments of 1 Veterinary and Comparative Anatomy,
Pharmacology and Physiology and 2 Biological Systems
Engineering, Washington State University, Pullman, Washington 99164;
3 Department of Physics, Division of Biophysics, Moscow
State University, Moscow; and 4 Institute
of Theoretical and Experimental Biophysics, Russian Academy of
Sciences, Puschino, Russia
 |
ABSTRACT |
A relatively simple method is presented for incorporating
cross-bridge mechanisms into a muscle model. The method is based on
representing force in a half sarcomere as the product of the stiffness
of all parallel cross bridges and their average distortion. Differential equations for sarcomeric stiffness are derived from a
three-state kinetic scheme for the cross-bridge cycle. Differential equations for average distortion are derived from a distortional balance that accounts for distortion entering and leaving due to
cross-bridge cycling and for distortion imposed by shearing motion
between thick and thin filaments. The distortion equations are unique
and enable sarcomere mechanodynamics to be described by only a few
ordinary differential equations. Model predictions of small-amplitude
step and sinusoidal responses agreed well with previously described
experimental results and allowed unique interpretations to be made of
various response components. Similarly good results were obtained for
model reproductions of force-velocity and large-amplitude step and ramp
responses. The model allowed reasonable predictions of contractile
behavior by taking into account what is understood to be basic muscle
contractile mechanisms.
mathematical model; muscle mechanics; muscle cross bridge; muscle
contraction
| |
INTRODUCTION |
THOSE WHO HAVE INTEGRATED muscle models into
engineering applications have followed various approaches with more or
less fidelity in representing contractile processes. As recounted in
several reviews (36, 41, 42), serious attempts to capture the essence of muscle behavior have followed one of two options: Hill-based models
[based on the experimental work of A. V. Hill (15)] or Huxley-based
models [based on the theoretical work of A. F. Huxley (18)]. The pros
and cons of Hill-based vs. Huxley-based muscle models have been
enumerated repeatedly (30, 37, 41, 42). Briefly, Hill-based models
incorporate phenomenological descriptions relating muscle force to
shortening velocity and, in addition, often include the relationship
between muscle tension and length. In contrast, Huxley-based models
attempt to derive overall muscle behavior from basic contractile
processes arising from sliding filament and actin-myosin interaction
mechanisms. Because of their purported computational simplicity,
Hill-based models have been the most popular (30, 37). However, the
desirability of incorporating basic contractile mechanisms into a
muscle model for application purposes, as with the Huxley-based
approach, is well recognized (42), and various authors have given
different means for reducing the computational complexity of
Huxley-based models (12, 33, 38-40). Despite the elegance of these
simplification attempts, a recent comparison between a Hill-based and a
simplified Huxley-based model concluded that present theories of
cross-bridge dynamics are not a suitable basis for development of
mathematical models for force production in human skeletal muscle (8,
35) and that Hill-based models remain the preferred approach.
Considering the immensity of our present understanding of contractile
processes, it seems that the problem associated with integrating that
knowledge into models for application purposes is more related to
building a suitable bridge between the molecular and the macroscopic
than to a deficiency in the understanding of molecular mechanisms.
Accordingly, the primary objective of this work was to simulate dynamic
behavior of muscle with mathematical equations that were inspired by
muscle contractile processes. In this sense, the model has its origins
in contractile mechanisms but ultimately exists as a descriptor and
predictor of general dynamic muscle phenomenology. This method of
dynamic physiological modeling is an advance over previous modeling,
based solely on limited phenomenological observations. Among the many
uses for the resulting model, one will be to reproduce in vivo muscle
function in animate movement simulations and to simulate muscle
function in actuators in the control of robots or limb prostheses.
| |
MODEL DESCRIPTION |
Glossary
f, f ', h, h', g, kon, koff
| A1, A2, D, R |
Cross-bridge states and number of cross bridges in attached
(A1, A2), detached (D), and regulatory (R)
states, respectively
|
| f, f ', h, h', g,
kon, koff |
Rate coefficients
|
| fr, f '0,
h0, h'0,
g0 |
Reference values of rate coefficients in isometric condition
|
| F(t) |
Total force from all cross bridges
|
| F |
Constant muscle force
|
| Fi(t) |
Force from cross bridges in i state
|
| F0 |
Isometric force
|
| e |
Stiffness of 1 cross bridge
|
| Ei = eAi |
Stiffness of all cross bridges in the i state
|
| xi |
Average distortion among cross bridges in i state
|
| Xi |
Total distortion among cross bridges in i state
|
| x0 |
Isometric value of average distortion among cross bridges in
A2 state
|
| RT |
Total number of cross bridges
|
| SL |
Sarcomere length
|
| SL0, L0 |
Reference or initial muscle length
|
 SL, L |
Changes in length around reference values
|
| ei(t) |
Variation in Ei induced by
L(t)
|
 xi(t) |
Variation in xi induced by
L(t)
|
| E0i |
Isometric values for total stiffness in i state
|
k  |
Probability of transition between  and  states
|
| E |
Activation energy for state transition
|
| k |
Boltzman's constant
|
| T |
Absolute temperature
|
 |
Stiffness divided by kT
|
 |
Velocity of filament sliding
|
 S |
Fraction of all cross bridges that are in the force-bearing state
|
| v |
Parameter of cooperativity among cross bridges
|
| j |
|
 |
Frequency of sinusoidal length changing
|
| Vmax |
Maximum shortening velocity
|
 |
Parameter in Hill equation for force-velocity relationship
|
Sarcomeric Stiffness and Distortion
The foundations of a muscle model reside in the representation of its
elemental contractile unit, the sarcomere. Our starting point for a
sarcomere model was that it may be represented in terms of its
stiffness and an average distortion among the stiffness elements. This
idea was taken from Berman et al. (2), where it was employed as a
deductive consequence of Huxley's assumption that muscle cross bridges
act as independent, parallel force generators (17). These ideas were
successfully used by us in a previous analysis (7) of the similarities
between the mechanical properties of the heart and heart muscle.
Briefly, muscle sarcomeres (Fig. 1) are
composed of overlapping thick and thin filaments and are bounded by
adjacent z lines. Thin filaments project from a z line toward the
middle of the sarcomere. Thick filaments are centered between z lines.
The two halves of a sarcomere are mirror images of one another 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
move preferentially toward either z line. Force is generated as a
result of interactions between thick and thin filaments through myosin
cross bridges. Myosin cross bridges extend out of the thick filament
and may be attached or not attached to their actin-binding sites on the
thin filament. Because cross bridges along a thick filament are
arranged in parallel with one another and because thick filaments
within a sarcomere also are in parallel with each other, cross bridges
within a half sarcomere are all in parallel. Thus the force generated
in a half sarcomere equals the sum of forces generated by attached
cross bridges in that half sarcomere.
View larger version (25K):
[in this window]
[in a new window]
|
Fig. 1.
Diagrammatic representation of sarcomere as the elemental contractile
unit, with myosin cross bridges as elemental force generators. Cross
bridges may be attached and force bearing or detached. Half sarcomere
on left and right must have equal force generation. All
force generators in a half sarcomere are in parallel.
|
|
Each cross bridge is considered to be linearly elastic. Force is
generated when parallel, linearly elastic cross bridges in attached
states are distorted. Cross bridges in detached states (i.e., those not
bound to the thin filament) cannot generate force. Attached cross
bridges (i.e., those capable of being distorted) may be in n
different states. Thus the force generated by all cross bridges within
a half sarcomere is given by
|
(1)
|
where e is the stiffness of a single cross bridge,
Ai(t) is the number of parallel cross
bridges in the ith of n attached state, and
xi(t) is the average distortion
among cross bridges in the ith state.
Both the Ai(t) and the
xi(t) are dynamic variables
that undergo changes with time to bring about dynamic variation in
F(t). Because e is a constant and equal for all cross bridges
in attached states, it may be brought inside the summation symbol to
form eAi(t) products, which physically
represent the stiffness associated with each of the i
populations within the n states,
eAi(t) = Ei(t). With this, Eq. 1
may be rewritten in a stiffness-distortion format
|
(2)
|
Importantly, the population stiffness
Ei(t) (unlike the single cross-bridge
stiffness e, which is a fixed value) is a dynamic variable that is
subject to both transient and steady-state variations. The dynamic
nature of this stiffness is a subject of detailed analysis in this work.
A dynamic sarcomere model consists of mathematical descriptions of
processes responsible for time variations in
Ei(t) and
xi(t). These processes may be
divided into several categories: 1) availability of activator
calcium; 2) changes in thick- and thin-filament overlap;
3) dynamics of myofilament regulatory proteins; and 4)
cross-bridge kinetics. Of these, categories 3 and 4 represent a dynamic core that is modulated by processes in categories
1 and 2. For this reason, this report will focus only
on processes within the dynamic core, i.e., myofilament regulation and
cross-bridge kinetics during constant calcium activation. It will be
shown that a sarcomere stiffness-distortion model based on myofilament
regulation and cross-bridge kinetics generates the dynamic behaviors
observed in the constantly calcium-activated muscle in response to
1) small-amplitude sinusoidal and step changes in muscle
length, 2) perturbations used in protocols to generate the
force-velocity relationship, and 3) large-amplitude step and ramp changes in muscle length.
Myofilament Regulation and Cross-Bridge Kinetics
We look for a blueprint for stiffness dynamics in some elementary
aspects of muscle contraction kinetics. Consider the myofilament regulation and cross-bridge cycling scheme in Fig.
2, consisting of an actin thin filament, a
regulatory tropomyosin-troponin complex, and a myosin cross bridge. The
thin-filament regulatory unit may be in the "off" position
(Roff) or the "on" position. Switching between off
and on positions is governed by the "on" rate coefficient (kon) and the "off" rate coefficient
(koff). These regulatory on and off rate
coefficients vary depending on the concentration of available activator
calcium, as given in Ref. 7. Here, we consider that activator calcium
concentration is constant.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 2.
Illustrative representation of myofilament regulation and cross-bridge
cycling giving rise to stiffness and distortion dynamics of sarcomere.
Actin thin-filament possessing the myosin-binding site is represented
by chain of 3 circles (no stoichiometric relations are implied by 3 circles). Thin-filament regulatory tropomyosin-troponin complex,
controlling myosin cross-bridge access to the thin-filament binding
site, is represented by the bar spanning the binding site. Myosin cross
bridge is represented by the globular head with coiled tail.
Thin-filament regulatory unit may be in 1 of 2 steric configurations:
the "off" position (Roff) or the "on"
position. The on state may have different status depending on whether
cross bridges are attached or not including "on" and cross
bridges detached (D) and "on" and cross bridges attached
(A1, A2). Attached cross bridges may be
pre-power stroke (A1) or post-power stroke
(A2). Cross bridge attachment, power stroke, and
detachment occur cyclically. Steps in the cycle are governed by rate
coefficients f, f ', h, h',
g, where the forward attachment is governed by f, the
forward power stroke is governed by h, and cross-bridge
detachment is governed by g. Primes designate reverse
reactions.
|
|
Cross bridges can be in attached (A1 and A2)
and detached (D) states. Attached cross bridges may be in the
pre-power-stroke (A1) and post-power-stroke
(A2) states. The power stroke, representing the
A1 to A2 transition, is the point in the cycle
of chemical-to-mechanical energy transduction. Cross-bridge attachment
can occur only when the regulatory unit is in the on configuration.
Attachment, power stroke, and detachment occur cyclically. Steps in the
cycle are governed by rate coefficients f, f',
h, h', and g, where the forward attachment is
governed by f, the forward power stroke is governed by
h, and cross-bridge detachment is governed by g. Primes
designate reverse reactions. Isometric force is generated as cross
bridges go through the power stroke.
Ordinary differential equations for the cross-bridge kinetic scheme of
Fig. 2 may be written from inspection as
|
(3)
|
|
(4)
|
|
(5)
|
where overdots indicate first time derivative,
Roff(t) = RT 
D(t) A1(t) A2(t), and RT is a constant
representing the total number of cross bridges for a given filament
overlap. Changes in filament overlap bring about changes in
RT and thus influence cross-bridge recruitment
through this mechanism. Such length-dependent changes will not be
investigated here but may be easily introduced into the model during
future studies.
Changes in A1(t) and
A2(t) result in direct changes in
E1(t) and E2(t)
according to arguments transforming Eq. 1 to Eq. 2 and,
thus, A1(t) and A2(t)
will be referred to as stiffness variables. Equations 3-5
provide the basis for dynamic changes in stiffness.
Cross-Bridge Distortion
Whereas the preceding cross-bridge kinetic scheme and resulting
differential equations are more or less similar to a great many
previously published schemes (see Refs. 9, 10, and 16 for equivalent
three-state cross-bridge cycle and Refs. 14, 22, 24, and 31 for similar
myofilament regulation schemes), there is no satisfactory treatment
available for dynamic changes in average cross-bridge distortion.
Because our expression for force (Eq. 2) requires an assessment
of average distortion, we provide an original treatment for average
distortion in the following.
There are two mechanisms for inducing distortion in cross bridges (Fig.
3). During isometric conditions, the power
stroke (i.e., the forward step in Fig. 2 regulated by the factor
h and shown schematically in Fig. 3A) induces a
distortion in the cross-bridge equivalent to the stretching of a spring
that, on the average among A2 cross bridges, is
x0. This x0 distortion is
lost from a single cross bridge both when it, in the A2
state, returns to the A1 state by the reverse power stroke
(regulated by h') and when it detaches during the
A2 to D transition (regulated by g). In addition,
changes in sarcomere length (SL) produce sliding movement
between thick and thin filaments, causing
x2(t) to vary from
x0 (Fig. 3B).
View larger version (49K):
[in this window]
[in a new window]
|
Fig. 3.
Illustration of 2 mechanisms for inducing distortion in cross bridges.
A: an attached cross bridge in the pre-power-stroke state
(A1) experiences no elastic distortion during isometric
conditions and does not generate force. Power stroke that converts
A1 state to A2 state is a chemical-mechanical
energy transduction step in which elastic energy is stored by
distorting cross bridge. Thus transformation from A1 to
A2 state is associated with genesis of an increment of
distortion, x0, in the A2 state.
B: when there is sliding motion between the thick and thin
filaments as during changes in sarcomere length, cross bridges in both
the A1 and A2 states undergo additional
distortion. Amount of distortion due to this shearing motion is
governed by balance between the velocity of sliding and the rate at
which elements within states are formed and dissipated (see Eqs.
11, 12 in text). Distortion induced by filament sliding is
transient and is dissipated by distorted cross bridges leaving the
population while being replaced by new cross bridges without
shear-induced distortion.
|
|
The differential equations for x2(t)
may be derived by performing a distortional balance over all parallel
cross bridges in the A2 state. The collective distortion
among the parallel A2 cross bridges differs from the force
associated with these cross bridges only in not including the stiffness
coefficient e. This collective distortion is given by
|
(6)
|
where the X2(t) is used to designate
the collective distortion and the x2(t)
designates the average distortion among the A2(t) cross bridges.
X2(t) at some t + t
can be written as
where
It is assumed that cross bridges enter the A2
state from the A1 state via the power stroke with
x0 distortion, regardless of whatever distortion
they may have possessed in the A1 state before the
transition to A2.
From the above, it can be written
![<IT>X</IT><SUB>2</SUB>(<IT>t</IT> + &Dgr;<IT>t</IT>) = <IT>X</IT><SUB>2</SUB>(<IT>t</IT>) + &Dgr;S<IT>L</IT> A<SUB>2</SUB>(<IT>t</IT>) + [<IT>h</IT> A<SUB>1</SUB>(<IT>t</IT>)]&Dgr;<IT>t</IT>](/content/vol87/issue5/fulltext/1861/img027.gif)
|
(7)
|
Rearranging gives
|
(8)
|
Taking the limit as t 
0, yields
![<A><AC>X</AC><AC>˙</AC></A><SUB>2</SUB>(<IT>t</IT>) = A<SUB>2</SUB>(<IT>t</IT>)S<A><AC>L</AC><AC>˙</AC></A>(<IT>t</IT>) + [<IT>h</IT> A<SUB>1</SUB>(<IT>t</IT>)]<IT>x</IT><SUB>0</SUB> − (<IT>g</IT> + <IT>h</IT>′)A<SUB>2</SUB>(<IT>t</IT>) <IT>x</IT><SUB>2</SUB>(<IT>t</IT>)](/content/vol87/issue5/fulltext/1861/img031.gif)
|
(9)
|
Now, 
2(t) may be
eliminated by noting that differentiation of Eq. 6 yields
|
(10)
|
Equating Eqs. 9 and 10, making appropriate
substitutions for 
2(t) from
Eq. 5, and solving for
2(t), gives the desired differential equation
![<IT><A><AC>x</AC><AC>˙</AC></A></IT><SUB>2</SUB>(<IT>t</IT>) = −<IT>h</IT> <FR><NU>A<SUB>1</SUB>(<IT>t</IT>)</NU><DE>A<SUB>2</SUB>(<IT>t</IT>)</DE></FR> [<IT>x</IT><SUB>2</SUB>(<IT>t</IT>) − <IT>x</IT><SUB>0</SUB>] + <A><AC>S</AC><AC>˙</AC></A><IT>L</IT>(<IT>t</IT>)](/content/vol87/issue5/fulltext/1861/img035.gif)
|
(11)
|
In like manner, it can be shown that
|
(12)
|
where cross bridges enter the A1 state from D
with no distortion and also return from A2 via the reverse
power stroke with no distortion.
Equations 11-12, describing the dynamic changes for
average distortion of the post-power-stroke and pre-power-stroke
states, are the key developments that allow a simplified model of
sarcomeric dynamics. We know of only one other approach (34) that
approximates that given here. However, the results from derivations in
Ref. 34 apply only to small-amplitude changes in SL, whereas
those given in Eqs. 11 and 12 apply generally to all
length perturbations of any amplitude.
Variable-Dependent Rate Coefficients
From the differential Eqs. 3-5 and 11-12, it
appears that SL(t) drives the system only through
its effect on cross-bridge distortion, with no impact at all on
stiffness variables, A1(t) and
A2(t). In fact, the stiffness variables may
respond to SL(t) through the indirect effect of
SL(t) on the rate coefficients. We introduce two of
the several potential mechanisms for this into our model, but in so
doing we encounter some theoretical difficulties.
Rate coefficient distortion dependence.
A central tenet in the cross-bridge theory of muscle contraction is
that rate coefficients governing state transitions depend on the
mechanical distortions (strains) experienced by the cross bridge (13,
17, 18, 26). Previous treatments of these distortion-dependent effects
have taken into account thermal effects causing a distribution of
distortions within a population and the impact of strain energy on
state transitions (see APPENDIX). One consequence of these
considerations is that the differential equation describing
distortion-dependent kinetic events becomes a partial differential
equation in the distortional variable. The resulting mathematical
problem takes on a fair measure of difficulty.
There have been various approaches for simplifying the problem of
distributed distortion among cross bridges. As outlined by Taylor et
al. (33), 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 cross-bridge
state were uniform. More direct have been approaches to formally
simplify the problem. Perhaps the first serious attempt was by
Deshcherevskii (11), who considered that above one value of distortion,
cross bridges were pulling in the direction of shortening, whereas for
distortion below this value, cross bridges 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 distortion 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. (12) where the distortion coordinate was discretized, and all
cross bridges within a discrete distortion range x ± x were grouped into a bin. All cross bridges 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 cross bridges entering and leaving each bin
during filament sliding. In accord with the approach of Deshcherevskii
(11), the number of bins employed by Dijkstra et al. (12) could be
reduced to as few as two, and reasonable results were obtained. A third
approach is to consider just small distortional changes such that
spatial average variations in the detachment rate coefficient are not
large and can be ignored (34). This changed the role of the
distortional variable in the myofilament system by removing its
function as a random variable, which allowed representing the system as
a system of ordinary differential equations. Another approach is that
of Zahalak (39, 40) where a distribution function for attached cross
bridges 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 quantities 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-crossbridge state needs to be considered.
Because the basic premise of our model was that muscle force equaled
sarcomeric stiffness times average cross-bridge distortion (Eq. 2), our challenge was to find some reasonable relationship between
population average rate coefficient and average distortion. We
recognized that, 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 (see APPENDIX). However, as discussed by Thorson and White (34),
variation in average distortion from its reference value would be
associated with variation in the average rate coefficient from its
reference value. To express this, it was convenient to adopt
parametrically succinct functional forms that satisfied the
experimental fact that the numbers of attached cross bridges decline
during constant shortening (20). This meant that, during shortening,
rate coefficients leading away from an attached state must increase
relative to those leading into that state. Additional considerations
resulted from adopting parameter values (Table
1), such that cross bridges in the
A1 state were short lived relative to cross bridges in the
A2 state (this is somewhat, but not totally, analogous with the commonly used concepts of weak-binding/strong-binding states). The
consequence was that, for a given
L(t),
x1 was never forced far from its isometric value of
0, whereas x2 was forced relatively far from its
isometric value of x0. Thus those rate
coefficients leading away from the A1 state were much less
affected by distortion than those leading away from the A2
state.
Given the above, attention was focused on distortion dependence in
g alone, and the following arbitrary expression was chosen as a
parametrically succinct candidate for introducing distortion-dependent effects into the model
|
(13)
|
In Eq. 13, the zero subscript refers to the value at
the isometric condition, and 
is a parameter that grades
distortional dependence, = 0 meant no distortional dependence.
Accordingly, when > 0, g increased whether
x2 was forced below x0, as
during shortening, or above x0, as during
stretching. Retrospectively, we found that best results were obtained
when took on a larger value when x2 > x0 than when x2 < x0. These different values for in these two
distortion ranges are given in Table 1.
To summarize, the stiffness-distortion approach to muscle modeling
leads to development of ordinary differential equation expressions for
average distortion of attached cross-bridge states. This average
distortion may then be used to introduce distortion dependence into
population-average kinetic rate coefficients. Our approach is only an
approximation because 1) no account is taken of the role of
distortion dependence in the derivation from the macroscopic balances
leading to Eqs. 3-5 and 11, 12 and these ordinary
differential equations are global approximations of the partial
differential equations that properly describe the detailed behavioral
features of the system; and 2) lacking a strict mathematical relationship between population-average quantities, the underlying thermodynamic constraints on kinetic relationships cannot 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.
Cooperative effects and rate coefficients.
A second mechanism for variable-dependent rate coefficients is through
cooperative effects (3). These cooperative effects may be of many
kinds, but we focus on just one kind, i.e., that which could occur via
force-bearing cross bridges. We introduced force-bearing cooperativity
into the model by letting force in one cross bridge enhance the
attachment rate coefficient f of a neighboring cross bridge.
This enhancement occurs by modification of the activation energy for
transition from detached D to attached A1 state (Eq. A1 in APPENDIX); an attached neighbor reduces the
activation energy for that transition. Each cross-bridge attachment
site has two adjacent neighbors, and the equation for the attachment rate coefficient must account for the impact of force-bearing cross bridges at both neighboring sites. The neighboring effect may be
represented in f by considering
|
(14)
|
The r subscript on the fr reference
value does not refer to the isometric value, as in distortion-dependent
g, but to the value when no neighbors are in the force-bearing state.
Using the forgoing considerations, it may be derived that
![<IT>f</IT> = <IT>f</IT><SUB>r</SUB> {1 + &lgr;<SUP>A<SUB>1</SUB></SUP> [e<SUP><IT>x</IT><SUB>1</SUB>/<IT>x</IT><SUB>0</SUB> (&ngr;−1)</SUP> − 1] + &lgr;<SUP>A<SUB>2</SUB></SUP> [e<SUP><IT>x</IT><SUB>2</SUB>/<IT>x</IT><SUB>0</SUB> (&ngr;−1)</SUP> − 1]}<SUP>2</SUP>](/content/vol87/issue5/fulltext/1861/img039.gif)
|
(15)
|
where A1 and
A2 are the probabilities
(= A1/RT and
A2/RT, respectively) of finding a
neighbor in one of the attached states, the exponential terms arise
from activation energy considerations, and is a number (
1)
expressing the influence of force-bearing in one cross bridge to affect
the attachment of its neighbor. When there is no cooperative influence,
= 1. In this way, a single index, , may be used to grade the
influence of nearest neighbor cooperative effects.
Because of the appearance of x in Eq. 15, the same
caveats and cautions as given above for distortion dependence must be
taken into consideration in this representation of cooperative effects.
Model Summary
In summary, the model consists of three differential equations
(Eqs. 3-5), describing dynamic changes in stiffness
variables, and two differential equations (Eqs. 11, 12),
describing dynamic changes in distortion variables. It also contains
equations for distortion-dependent rate coefficients (Eq. 13)
and cooperative effects (Eq. 19). Taken together, these
equations constitute the stiffness-distortion model. This is a
nonlinear set of equations of fifth dynamic order in the state
variables D, A1, A2,
x1, and x2. The forcing
function is SL(t), which appears explicitly (as the
first time derivative) only in the distortion equations (Eqs. 11,
12) but enters implicitly through the variable-dependent rate coefficients, Eqs. 13 and 19. The system output is
F(t), which, as derived from Eq. 2, is
|
(16)
|
where E1(t) = eA1(t) and E2(t) = eA2(t).
The model differential equations are nonlinear in that they contain
products and ratios of state variables at several locations. In
addition, the equations are nonlinear because the rate coefficients, which appear in the equations as multipliers of the state variables, depend on the distortion variables as a consequence of
distortion-dependent effects and on variables of the cycle state
because of cooperative effects. In its full nonlinear form, the model
contains 10 parameters including 7 reference values for the rate
coefficients (kon, koff, fr, f '0,
h0, h'0,
g0); an index for distortion dependence (); an
index for cooperativity (); and a value for isometric distortion (x0).
| |
MODEL PROTOCOLS AND RESULTS |
The purpose of model protocols was to elicit dynamic features of the
model that could be compared to experimentally determined dynamic
features of constantly activated muscle. The initial focus was on
history-dependent dynamic features arising from differential equation
relations among variables. This was because such history-dependent dynamics are important features of constantly calcium-activated muscle,
and these cannot be reproduced by Hill-based models. Second, the
differential-equation-based model was examined for its ability to
reproduce hyperbolic force-velocity behavior. Such special-case behavior is the essence of Hill-based models. Third, the model was
asked to create large-amplitude behavior in which features of both
history-dependent and force-velocity behaviors could be observed.
Finally, to demonstrate its versatility, the model was modified to
change its dynamic signature from one that is characteristic of
skeletal muscle to one that is characteristic of cardiac muscle.
Linear Dynamic Response
The standard figures of merit for evaluating dynamic behavior of
constantly activated muscle are the frequency-dependent complex stiffness determined by force responses to small-amplitude sinusoidal length changes (21, 23) and the force response to step length changes
(1, 18). Accordingly, model-predicted small-amplitude behaviors were
obtained by linearization procedures in which small deviations around a
steady-state isometric value were considered, and then truncation was
performed to eliminate second-order terms. Thus, for small length
variation, L(t), around some baseline length,
L0, Eq. 16 may be expanded to
give
![= (E<SUP>0</SUP><SUB>1</SUB><IT>x</IT><SUP>0</SUP><SUB>1</SUB> + E<SUP>0</SUP><SUB>2</SUB><IT>x</IT><SUP>0</SUP><SUB>2</SUB>) + [E<SUP>0</SUP><SUB>1</SUB>&dgr;<IT>x</IT><SUB>1</SUB>(<IT>t</IT>) + e<SUB>1</SUB>(<IT>t</IT>) <IT>x</IT><SUP>0</SUP><SUB>1</SUB> + E<SUP>0</SUP><SUB>2</SUB> &dgr;<IT>x</IT><SUB>2</SUB>(<IT>t</IT>) + e<SUB>2</SUB>(<IT>t</IT>) <IT>x</IT><SUP>0</SUP><SUB>2</SUB>] + [e<SUB>1</SUB>(<IT>t</IT>)&dgr;<IT>x</IT><SUB>1</SUB>(<IT>t</IT>) + e<SUB>2</SUB>(<IT>t</IT>)&dgr;<IT>x</IT><SUB>2</SUB>(<IT>t</IT>)]](/content/vol87/issue5/fulltext/1861/img042.gif)
|
(17)
|
where e1(t) is the
variation in E1 induced by L(t), and
x1(t),
e2(t), and
x2(t) are corresponding variations
induced in the respective variables, x1,
E2, and x2. Superscript 0 refers to
the value at the reference isometric condition. The three terms in
square brackets on the right-hand side represent zero-order, first-order, and second-order effects, respectively. Truncation of
Eq. 17 to eliminate second-order effects and substitution
of the reference isometric values
x02 = x0 and x01 = 0 yield, for the isometric
force
|
(18)
|
and, for the linear small-amplitude variation in
F(t) around F0
|
(19)
|
Note that e1(t) does not enter into
F(t), because the e1(t)
contribution is multiplied by
x01, which is zero.
Frequency-dependent complex stiffness.
The concept of stiffness embodied in the time-domain variables
E1(t) and E2(t) is now
generalized to the complex quantity F/L in the frequency
domain. Taking the Fourier transform of Eq. 19 and dividing
through by the Fourier transform of L(t) gives the incremental complex stiffness as
|
(20)
|
where
|
(21)
|
|
(22)
|
|
(23)
|
Equations 21 and 22 were derived from
linearized versions of Eqs. 11 and 12, and Eq.
23 was derived from the linearized, single third-order differential
equation in e2(t), resulting from the combination of the three cross-bridge state Eqs. 3-5. The
ki in Eq. 23 are combinations of
the rate coefficients f, f ', h, h',
and g; these relations are given in the APPENDIX.
Distortion-dependent features of g do not enter into the
ki in this small-amplitude expression,
because distortion dependence was formulated as affecting g
only as a result of deviation from isometric (steady-state) values.
These deviations fall out in the linearization procedure. Thus, the
ki depend just on the reference values of
g. In contrast, steady-state values depend strongly on the
cooperative effects in f. Consequently, the
ki do depend on the cooperative effects in
f as they appear in the truncated series expansion. Through
these effects, the distortional variables x1 and
x2 appear in the numerator of Eq. 23 as a
result of the dependence of f on x1 and
x2 in Eq. 15.
Model-predicted complex stiffness consists of three components
corresponding to the three terms on the right-hand side of Eq. 20. It is important to note that the physical units for each of the
three terms derive from the stiffness variable as it was used in
Eq. 2. However, the dynamic features in two of these three complex-stiffness components derive from dimensionless dynamics of
distortional variables. We refer to each of these components as the
x1, x2, and
e2, respectively, because: 1) the frequency dependence of the first component derives from the dynamic response of
x1 to changes in length; 2) the frequency
dependence of the second component derives from the dynamic response of
x2 to changes in length; and 3) the
frequency dependence of the third component derives from the dynamic
response of e2 to changes in length.
Utilizing a model parameter set (Table 1) that yielded negative-phase
complex stiffness at frequencies within the 1- to 10-Hz range, the
polar plot of the model-predicted overall complex stiffness and of each
of its component parts was determined as shown in Fig.
4. The model-predicted frequency variation
in the overall complex stiffness (Fig. 4, A) is reminiscent
of experimentally measured complex stiffness of skeletal muscle (4, 19)
in exhibiting a looping locus that may be divided into three regions: a
low-frequency positive-phase region (<1 Hz); a midfrequency negative-phase region (for frequencies in a range roughly between 1 and
10 Hz); and a high-frequency positive-phase region (>10 Hz). These
various regions reflect changing contributions of the three components,
i.e., the x1, x2, and
e2 components from Eq. 20. For instance, compare
the relative magnitudes of the complex vector for each of the three
components at 1 Hz (Fig. 4, A, B, and C).
From these relative magnitudes, it can be concluded that the
low-frequency positive-phase region of the overall complex stiffness is
dominated by the e2 and x2 components
(each with large vectors at 1 Hz), with the x1
component (with an extremely small vector at 1 Hz) playing a negligible
role. The midfrequency negative-phase region is due to the predominance
of the imaginary part of the e2 component; it being the
only component with a negative imaginary part and a locus that enters
the fourth quadrant. The transition from negative phase to positive
phase as frequencies increase above 10 Hz is due to the increasing
importance of dynamics associated with the x1
component.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 4.
Polar plots of model-generated overall complex stiffness (A)
and its 3 dynamic components (B, C, and D)
over frequency range 0-50 Hz. Distortion dynamics contribute to
complex stiffness through the x1 and
x2 components. The e2 component
expresses dynamics of cross-bridge recruitment. Frequencies of 1 Hz and
10 Hz ( ) are marked on polar plots. Magnitudes of complex vector at
1 Hz in each of 3 component polar plots are indicated. Vector in the
overall complex stiffness is the sum of vectors in the 3 components.
Contribution of each component to overall complex stiffness changes
with changes in the respective component vectors at each frequency.
, Frequency of sinusoidal length changing; Im, imaginary; Re,
real.
|
|
The genesis of complex stiffness in terms of relative contributions of
the three above components is unique to the present model and contrasts
with previous reproductions of experimentally determined stiffness
dynamics using sums of exponentials (1) or sums of first-order filters
(43). At least two important interpretive points are to be made from
the present model results. The first of these concerns the mechanisms
responsible for the negative-phase region of the overall complex
stiffness. The negative-phase frequencies are frequencies of positive
work (32, 34) in which the subject (muscle or model) is capable of
doing work on its (actual or simulated) mechanical environment. At all
frequencies of positive-phase complex stiffness, the mechanical
environment performs work on the subject. The negative-phase
frequencies in the overall model (i.e., the frequencies of positive
work) arise because of the dominance of e2 dynamics at
these frequencies. Importantly, these e2 dynamics are
representative of the entire myofilament-regulation/cross-bridge-cycle
system. This may be appreciated by recognizing that the roots of the
characteristic equation of Eq. 23 represent the essential
dynamic features of the e2 component response and depend on
combinations of the ki. In turn, the
ki depend on products and sums of all the rate coefficients (APPENDIX) in such a way that no single
rate coefficient dominates in determining any individual eigenvalue.
Thus myofilament regulatory and cross-bridge cycling processes combine
in a complicated fashion to determine overall model dynamics at
frequencies in the negative-phase, or positive-work, region of the
complex stiffness.
This demonstrates that the power stroke has no more obvious influence
on the dynamics of positive work than any other step in the cycle.
Instead, positive work is the result of a recruitment phenomenon where
the relative rates of cross-bridge association exceed those of
dissociation, and the A2 state (isometric force-bearing state) is transiently populated to a greater extent than during the
initial isometric conditions (which is clearly demonstrated in the
delayed tension response to a step shown in Figs.
5D and 6C). The kinetics of the
recruitment phenomenon responsible for negative-phase frequencies of
the complex stiffness should not be confused with the kinetics of the
power stroke.
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 5.
Changes in model-generated complex-stiffness polar plot (A)
plus its e2 component (B) and the step response
(C) plus its e2 component (D) with
varying degrees of cooperativity. In all panels, thick lines are
results with strong cooperative effect (v = 3); thin lines
are results with moderate cooperative effect (v = 2.7); and
dashed lines (which reduce to the dot at the origin in B) are
results with no cooperative effect (v = 1). Complex stiffness
at 1 and 10 Hz are indicated by  and , respectively. Among the 3 components of the overall dynamic response (A, B,
C, in Figs. 5 and 7), cooperativity impacts just dynamic
changes in the e2 component. Thus e2 component
alone is responsible for changes in the overall response with changes
in cooperativity. See Glossary for symbol definitions.
|
|
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 6.
Model-generated small-amplitude step response (A) and its
dynamic components (B, C, and D). This
time-domain step response and its components were calculated directly
from frequency-domain complex stiffness and its respective components
given in Fig. 4. Overall response is the sum of the 3 components.
|
|
Whereas the dynamic features of the recruitment phenomenon are
determined by dynamic features of the entire system, the origin and
magnitude of the recruitment response are due entirely to the
cooperative effect in f. This can be appreciated by examining the influence of changes in , the coefficient expressing strength of
the cooperative effect, on complex stiffness. In Fig. 5A, the polar plot of the overall complex stiffness has been normalized to the
steady-state stiffness, and changes for different values of are
given that correspond to no cooperative effect ( = 1), modest
cooperative effect ( = 2.7), and strong cooperative effect ( = 3). The polar plots for the associated e2
component are given in Fig. 5B. Clearly, the magnitude of the
e2 component increases with increasing and, since this
is the only component of the overall response that is affected by ,
these changing magnitudes of the e2 component are entirely
responsible for the looping behavior of the overall response that leads
to negative-phase (and positive-work) behaviors.
A second important interpretive point arises from the increasing
importance and eventual dominance of dynamics associated with the
x1 component at frequencies approaching and
surpassing 10 Hz. This implicates an important mechanical role for the
pre-power-stroke state at frequencies that are within the range of
frequencies of physiological importance. Consequently, the
pre-power-stroke state is not weakly bound in the usual sense (4, 5)
and may not be lumped with detached states as in two-state cycle
schemes. As others have found (34), a three-state cross-bridge cycle appears to be a minimum configuration for reproducing the full frequency spectrum of the characteristic frequency-dependent complex stiffness.
Step response.
The time domain step response (Fig. 5C and Fig. 6A)
and respective component responses (Fig. 6, A-C)
corresponding to the polar plots in Fig. 4 were obtained by operations
on Eqs. 25-27 to allow solving for the respective
time-domain dynamic variables. This was done by rearranging Eqs.
25-27 to solve explicitly for F( j),
x1( j),
x2( j), and
e2( j), assigning
L( j) = 1/j, and then taking the
inverse Fourier transform to solve for F(t),
x1(t), x2(t), and
e2(t).
The small-amplitude model-predicted step response had important
features in common with experimentally measured step responses, as
described in Refs. 1, 18. Both model-predicted and experimentally measured responses possess a leading edge coincident with the leading
edge of the length step, a rapid recovery phase, and a delayed tension
transient. In Fig. 5C, the dashed line is the sum of just the
x1 and x2 component responses.
From this, it can be seen that the leading edge and rapid recovery
phases are due to the x1 and x2
component responses, i.e., they are due to distortional responses. The
individual x1 and x2 component
responses are qualitatively similar but quantitatively different. Their
leading edges differ by the multiplying factors
A01 and A02 of Eq. 20 and, from Eqs. 25, 26, their time courses of
recovery are controlled by rate coefficients regulating transitions
away from the respective states. The rate coefficient for
x1 recovery, h + f ' (= 408
s
1 from Table 1), is much larger than the rate
coefficient for x2 recovery, g + h' (= 10 s1 from Table 1). Thus, in line with
the earlier argument about relative disturbances in
x1 and x2 values at different
velocities, length-change-induced distortion of the pre-power-stroke
state recovers much faster than that for the post-pow
TOP
ABSTRACT
INTRODUCTION
![(no. A<SUB>2</SUB> cross bridges existent at <IT>t</IT>)]](/content/vol87/issue5/fulltext/1861/img015.gif)
![(distortion of newly formed cross bridges)]](/content/vol87/issue5/fulltext/1861/img020.gif)
![= [<IT>h</IT> A<SUB>1</SUB>(<IT>t</IT>)] &Dgr;<IT>t</IT>(<IT>x</IT><SUB>0</SUB> + &Dgr;S<IT>L</IT>/2)](/content/vol87/issue5/fulltext/1861/img021.gif)
![(average distortion of these lost A<SUB>2</SUB> cross bridges)]](/content/vol87/issue5/fulltext/1861/img025.gif)
![= [(<IT>g</IT> + <IT>h</IT>′) A<SUB>2</SUB>(<IT>t</IT>)] &Dgr;<IT>t x</IT><SUB>2</SUB>(<IT>t</IT>)](/content/vol87/issue5/fulltext/1861/img026.gif)
![· <FENCE><IT>x</IT><SUB>0</SUB> + <FR><NU>&Dgr;S<IT>L</IT></NU><DE>2</DE></FR></FENCE> − [(<IT>g</IT> + <IT>h</IT>′)A<SUB>2</SUB>(<IT>t</IT>)]&Dgr;<IT>t x</IT><SUB>2</SUB>(<IT>t</IT>)](/content/vol87/issue5/fulltext/1861/img028.gif)
![+ [<IT>h</IT> A<SUB>1</SUB>(<IT>t</IT>)] <FENCE><IT>x</IT><SUB>0</SUB> + <FR><NU>&Dgr;S<IT>L</IT></NU><DE>2</DE></FR></FENCE> − [(<IT>g</IT> + <IT>h</IT>′)A<SUB>2</SUB>(<IT>t</IT>)]<IT>x</IT><SUB>2</SUB>(<IT>t</IT>)](/content/vol87/issue5/fulltext/1861/img030.gif)
![F(<IT>t</IT>) = F<SUB>0</SUB> + &Dgr;F(<IT>t</IT>) = [E<SUP>0</SUP><SUB>1</SUB> + e<SUB>1</SUB>(<IT>t</IT>)][<IT>x</IT><SUP>0</SUP><SUB>1</SUB> + &dgr;<IT>x</IT><SUB>1</SUB>(<IT>t</IT>)] + [E<SUP>0</SUP><SUB>2</SUB> + e<SUB>2</SUB>(<IT>t</IT>)][<IT>x</IT><SUP>0</SUP><SUB>2</SUB> + &dgr;<IT>x</IT><SUB>2</SUB>(<IT>t</IT>)]](/content/vol87/issue5/fulltext/1861/img041.gif)
