INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION MODEL PROTOCOLS AND RESULTS REFERENCES APPENDIX

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

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION MODEL PROTOCOLS AND RESULTS REFERENCES APPENDIX

Glossary

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

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, A_i(t) is the number of parallel cross bridges in the ith of n attached state, and x_i(t) is the average distortion among cross bridges in the ith state.

Both the A_i(t) and the x_i(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 eA_i(t) products, which physically represent the stiffness associated with each of the i populations within the n states, eA_i(t) = E_i(t). With this, Eq. 1 may be rewritten in a stiffness-distortion format

(2)

Importantly, the population stiffness E_i(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 E_i(t) and x_i(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

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 (A₁ and A₂) and detached (D) states. Attached cross bridges may be in the pre-power-stroke (A₁) and post-power-stroke (A₂) states. The power stroke, representing the A₁ to A₂ 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, R_off(t) = R_T  D(t)  A₁(t)  A₂(t), and R_T is a constant representing the total number of cross bridges for a given filament overlap. Changes in filament overlap bring about changes in R_T 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 A₁(t) and A₂(t) result in direct changes in E₁(t) and E₂(t) according to arguments transforming Eq. 1 to Eq. 2 and, thus, A₁(t) and A₂(t) will be referred to as stiffness variables. Equations 3-5 provide the basis for dynamic changes in stiffness.

Cross-Bridge Distortion

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 A₂ cross bridges, is x₀. This x₀ distortion is lost from a single cross bridge both when it, in the A₂ state, returns to the A₁ state by the reverse power stroke (regulated by h') and when it detaches during the A₂ to D transition (regulated by g). In addition, changes in sarcomere length (SL) produce sliding movement between thick and thin filaments, causing x₂(t) to vary from x₀ (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 x₂(t) may be derived by performing a distortional balance over all parallel cross bridges in the A₂ state. The collective distortion among the parallel A₂ 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 X₂(t) is used to designate the collective distortion and the x₂(t) designates the average distortion among the A₂(t) cross bridges. X₂(t) at some t + t can be written as






where













It is assumed that cross bridges enter the A₂ state from the A₁ state via the power stroke with x₀ distortion, regardless of whatever distortion they may have possessed in the A₁ state before the transition to A₂.

From the above, it can be written

(7)

Rearranging gives

(8)

Taking the limit as t  0, yields

(9)

Now, ₂(t) may be eliminated by noting that differentiation of Eq. 6 yields

(10)

Equating Eqs. 9 and 10, making appropriate substitutions for ₂(t) from Eq. 5, and solving for ₂(t), gives the desired differential equation

(11)

In like manner, it can be shown that

(12)

where cross bridges enter the A₁ state from D with no distortion and also return from A₂ 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

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.

(13)

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 A₁ 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)

(15)

Model Summary

(16)

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 (k_on, k_off, f_r, f '₀, h₀, h'₀, g₀); an index for distortion dependence (); an index for cooperativity (); and a value for isometric distortion (x₀).

	MODEL PROTOCOLS AND RESULTS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION MODEL PROTOCOLS AND RESULTS REFERENCES APPENDIX

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

(17)

where e₁(t) is the variation in E₁ induced by L(t), and x₁(t), e₂(t), and x₂(t) are corresponding variations induced in the respective variables, x₁, E₂, and x₂. 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 x⁰₂ = x₀ and x⁰₁ = 0 yield, for the isometric force

(18)

and, for the linear small-amplitude variation in F(t) around F₀

(19)

Note that e₁(t) does not enter into F(t), because the e₁(t) contribution is multiplied by x⁰₁, which is zero.

Frequency-dependent complex stiffness. The concept of stiffness embodied in the time-domain variables E₁(t) and E₂(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 e₂(t), resulting from the combination of the three cross-bridge state Eqs. 3-5. The k_i in Eq. 23 are combinations of the rate coefficients f, f ', h, h', and g; these relations are given in the APPENDIX.

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.

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.

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), x₁( j), x₂( j), and e₂( j), assigning L( j) = 1/j, and then taking the inverse Fourier transform to solve for F(t), x₁(t), x₂(t), and e₂(t).

1

1

Force-Velocity Relationship

Model-predicted force-velocity relationships during shortening and stretch were made by using the same set of parameters (Table 1) as in the small-amplitude analysis. Results are presented in Fig. 7 as normalized for V_max and F₀. These results are examined in two regions: force-shortening velocity (V/V_max > 0), and force-stretching velocity (V/V_max < 0)

View larger version (14K): [in this window] [in a new window]   Fig. 7.   Model-predicted force-velocity relationship (solid line) and its comparison with theoretical Hill-equation (dashed line, parameter of the Hill curve k = 0.23). Note good agreement during shortening (V/Vmax) > 0, where V is shortening velocity and Vmax is maximal V) and deviation of the 2 curves during stretching (V/Vmax < 0).

In accord with experimental data obtained from a wide variety of muscles, the force-shortening velocity portion of the model-predicted curve was closely approximated by the hyperbolic curve of the classic Hill force-velocity equation

(24)

with  = 0.15 for the closely fitting dashed curve in Fig. 7. Typical values of  for muscle of many types range between 0.15 and 0.25 (25). Thus the value of  that gave a close approximation to the model-predicted force-velocity relationship was within the accepted range. Given the close fit, similar to the fits obtained with experimental data from muscle and the acceptable  value, it can be concluded that the model generated satisfactorily realistic force-velocity behavior.

The model-predicted force-stretching velocity region also shares important characteristics in common with muscle, in demonstrating an initial shallow dependence of velocity on force for small clamps above F₀ and then increasingly larger velocities for larger force clamps. This behavior is unlike that in the force-shortening velocity region and cannot be approximated by the classic Hill force-velocity curve. Similar to actual muscle, "yielding" (19) was observed in the model predictions as force approached an apparent asymptotic value somewhere in the neighborhood of 2 F₀. A clear asymptote could not be defined by our methods, because the numerical solution procedure became unreliable as the asymptote was approached. Nonetheless, evidence for yielding can clearly be seen in Fig. 7, in those parts of the curve that were obtained at the highest force levels allowed by our numerical procedures.

We examined features of the model responsible for the force-velocity predictions. Whereas distortion-dependent g played no role in the small-amplitude sinusoidal and step-response model predictions (these small-amplitude behaviors were strongly impacted by the cooperative effects in f, as shown in Fig. 5), distortion-dependent g played an important role in the force-velocity relationship. We sought to determine the relative contributions of cooperativity in f and distortion dependence in g in the force-velocity relationship.

When f and g were frozen at their values in the isometric condition, the model predicted a linear force-velocity relationship, as given by the straight line in Fig. 8. In this case, the maximal velocity of shortening was much less (65%) than what was obtained with the full model. When only f was allowed to vary by expressing cooperative effects, the dotted curve below the straight line in Fig. 8 was obtained. From this it can be seen that these cooperative effects in f, by themselves, tend to reduce the velocity at any given load. This effect is most pronounced at middle loads, and there is no effect on the V_max at zero load. Furthermore, these cooperative effects in f did not result in a detectable yielding behavior in the force-stretching velocity region, as velocity retained only a shallow dependence on F during stretching. In contrast, when f was frozen at its isometric value and only g was allowed to change with distortion, the dashed curve in Fig. 8 was obtained. From this, it is seen that the distortion dependence in g had its greatest effect during shortening at low loads where it brought about a considerable increase in shortening velocity, including a substantial increase in V_max. In addition, distortion dependence in g caused a pronounced yielding effect with an apparent force asymptote at a load only slightly higher than isometric force.

View larger version (20K): [in this window] [in a new window]   Fig. 8.   Model-predicted force-velocity relationship (thick solid line) and its comparison with predictions when there was no variable-dependent rate coefficients (thin solid line); cooperative effect in f only (dashed line); and distortion-dependence in g only (dotted line).

When both cooperative effects in f and distortion dependence in g were taken together, the resultant curve was bounded by the effects of f alone and of g alone. At the lowest loads during shortening and the highest loads during stretching, the effects of distortion dependence in g predominate to elevate V_max and cause yielding. At intermediate loads (0.5-1.5 F₀), the effects of cooperativity in f predominated. The transition between dominance by distortion dependence in g and cooperative effects in f resulted in a curve that was much more curvilinear than would be obtained by either effect alone.

To test the validity of our assumption that distortion dependence in g was the most important distortion-dependent effect among the rate coefficients, the remaining rate coefficients were assigned distortion dependence according to

(25)

(26)

(27)

When these, all with the same distortion-dependent coefficient , were incorporated into the model, and force-velocity was predicted, it was found that, relative to distortion dependence only in g, distortion dependence in all rate coefficients caused an increase in V_max, a greater curvilinearity during shortening, and a slightly greater tendency for yielding during stretching. These effects, due to distortion dependence in all coefficients, were a matter of degree only, and it was concluded that, for all intents and purposes, the effects of distortion dependence in g alone sufficed to enable the model to predict the important behaviors due to distortion-dependent rate coefficient effects.

In addition to the influence on the force-velocity relationship of varying rate coefficients with cooperative effects and distortion dependence, we investigated the sensitivity of this relationship to fixed quantities of selected rate coefficients. In Fig. 9, we compared the sensitivity to the forward rate coefficient governing the power stroke, h, to that governing detachment, g. The reference curve of Fig. 7 (calculated with the parameters of Table 1: h = 8 s¹, g₀ = 4 s¹) is reproduced in Fig. 9, A and B, as the solid curve. In Fig. 9A, results from variation in h were obtained by setting h to be variously 12 and 4 s¹ and by calculating new force-velocity curves. These were then graphed by using V_max and F₀ of the reference curve as normalizing factors. Changes in h had a large effect on curve orientation and isometric force but very little effect on V_max. A decrease in h shifted the curve down so as to decrease F₀; an increase in h shifted the curve up so as to increase F₀. Even though changes in h did not change V_max appreciably, these brought about marked changes in shortening velocity for all F > 0. This is not unexpected as one would anticipate shortening velocity to be influenced by the speed of the power stroke (likened to the cross-bridge rowing motion that propels filaments past one another).

View larger version (14K): [in this window] [in a new window]   Fig. 9.   Sensitivity of force-velocity curve to changes in rate coefficients h (A) and g (B). Solid lines in A and B are the reference force-velocity curve using parameters in Fig. 6 and in Table 1. Velocity and force are scaled to their maximum values calculated by using reference set of parameters: h0 = 8 s1, g0 = 4 s1. A: dashed line h0 = 12 s1, dotted line h0 = 6 s1. B: dashed line g0 = 2 s1, dotted line g0 = 6 s1.

In Fig. 9B, results from variation in isometric g were obtained according to the following: g₀ was made to be variously 6 and 2 s¹, distortion-dependent g was calculated according to Eq. 13, and the resulting force-velocity curves were normalized as was done in Fig. 9A. An increase in isometric g caused a decrease in isometric force and an increase in V_max, whereas a decrease in isometric g had opposite effects. Thus isometric force is enhanced by changes in kinetic coefficients that promote the formation of the A₂ state (increased h and/or decreased g) while V_max is enhanced strongly by increased A₂ breakage due to increased g but only weakly by increased h. The sensitivity of model-predicted force-velocity relations to model parameters suggests that the model may be used to describe and interpret force-velocity data obtained from muscles of very different characteristics.

Large-Amplitude Steps and Ramps

Large-amplitude step responses differ from small-amplitude step responses in the rate of recovery from the distortional response, in the relative amplitude and time course of the delayed transient, and in asymmetry of responses to stretch relative to those of release (Fig. 10). For release responses, the classic four phases of the response may be identified as in Ref. 17. Force response to large-amplitude ramp stretches differs from force responses to ramp shortening (Fig. 11) in exhibiting yielding during ramp stretches, which limits the magnitude of force (see especially curve a in Fig. 11). Yielding in the response to positive ramps is as would be expected from the yielding observed in the force-stretching velocity relationship in Figs. 7 and 8. Thus large-amplitude behaviors, obtained from the model by numerical integration, exhibit features different from, but consistent with, those derived analytically for small-amplitude steps and with those derived from steady-state solutions for force-velocity relationships.

View larger version (13K): [in this window] [in a new window]   Fig. 10.   Time course of force response to step changes in length with the amplitude 0.4 * x0, 0.2 * x0, 0.1 * x0 (A) and 0.5 * x0, 0.3 * x0, 0.1 * x0 (B). Nos. 1-4 correspond to phases of response, as originally specified by Huxley and Simmons (18).

View larger version (20K): [in this window] [in a new window]   Fig. 11.   Force response to large-amplitude ramp stretches. Amplitude of length changes ±x0, duration: a = 0.1 s, b = 0.3 s, c = 0.5 s.

These realistic large-amplitude step and ramp responses, in conjunction with realistic model-predicted force-velocity relationships, lead us to conclude that the model may, in general, be used with confidence in the prediction of muscle behavior in response to large-amplitude perturbations of all kinds.

Changes in Dynamic Profile to Represent Muscles of Different Types

A particularly dramatic change in the model's complex-stiffness profile is achieved by adding a length-sensing feature to the attachment rate constant f, as shown in Fig. 12. Length sensing in f was introduced according to

(28)

where f_r is the f_r of Eq. 15, f_m is a constant, SL₀ is a reference sarcomere length, and  is the parameter that scales the influence of length sensing in f. Similar to the findings of Thorson and White (34), increased length sensing in f moves the zero-frequency stiffness from the origin to increasingly rightward locations on the positive real axis (Fig. 12). Additionally, the rightward movement along the positive real axis reduces the phase angle and frequency range of the positive-phase low-frequency region of the polar plot while increasing the phase angle and frequency range of the negative-phase region (Fig. 12, locus b). With sufficient length sensing in f, the low-frequency positive-phase region disappears, and the polar plot begins with a negative-phase leftward sweep through the fourth quadrant before eventually passing into the first quadrant at higher frequencies (Fig. 12, locus c). Polar plots like locus c are characteristic of cardiac muscle and asynchronous insect flight muscle (2, 7, 23, 27-29, 34). Changes in the overall complex stiffness due to length sensing are entirely the result of corresponding changes induced in the e₂ component, and there is no effect on the x₁ and x₂ components. These effects on the e₂ component exhibit in the time domain as slowed, enhanced, and sustained delayed tension; an effect known as "stretch activation."

View larger version (15K): [in this window] [in a new window]   Fig. 12.   Effect of length sensing on model-generated complex stiffness. Length-dependent coefficient  was varied as follows:  = 0, a;  = 3, b;  = 5, c. , indicate 1-Hz and  indicate 10-Hz positions on various loci. Za, Zb, Zc: zero-frequency real-axis intercept for loci a, b, and c, respectively. Note that the position of zero frequency stiffness on real axis shifts increasingly rightward as  is increased. Furthermore, negative-phase region of locus is enhanced. Through these length-sensing features, the model's complex stiffness changes from one characteristic of skeletal muscle (locus a) to one characteristic of cardiac and insect flight muscle (locus c).

Thus the model contains features that allow flexibility in reproducing dynamic behavior of constantly activated muscles with very different types of contractile behavior, from skeletal muscle of different speeds to muscles such as cardiac muscle and asynchronous insect flight muscle, which exhibit pronounced stretch activation (27).

Summary

Furthermore, the model is structured so that it may easily be extended to applications where time-varying calcium activation and changes in myofilament overlap are important. The former can be achieved by considering k_on and k_off as function of calcium concentration as in Ref. 7. Supplying the appropriate time-varying myoplasmic free calcium to appropriate expressions for k_on and k_off will result in a simulated twitch contraction or brief tetanus, as desired. Myofilament overlap effects are achieved by supplying appropriate length-dependent R_T as required for Eqs. 3-5. Also, complications concerning thin-filament regulatory processes, such as cooperative feedback effects as treated in Ref. 6, may be accommodated through an appropriate functional dependence of the rate coefficients k_on and k_off on any of the several model variables. Whereas these complications involving contractile processes may be easily introduced through the present structure of the model, complications such as noncontractile series elastance and spatial inhomogeneity will require more elaborate model structures in which the proposed sarcomere model would be one component.

In conclusion, we present a model of the muscle sarcomere that is sufficiently simple that it may be used in application settings where simulation of natural muscle behavior is desired for purposes of control, engineering design, and greater understanding.