## Abstract

Analysis of the hystereses in the force-length relationship at constant Ca^{2+} concentration and in the force-calcium relationship at constant sarcomere length (SL) provides insight into the mechanisms that control cross-bridge (XB) recruitment. The hystereses are related here to two mechanisms that regulate the number of strong XBs: the cooperativity, whereby the number of strong XBs determines calcium affinity, and the mechanical feedback, whereby the shortening velocity determines the duration for which the XBs are in the strong state. The study simulates the phenomena and defines the role of these feedbacks. The model that couples calcium kinetics with XB cycling was built on Simulink software (Matlab). Counterclockwise (CCW) hysteresis, wherein the force response lags behind the SL oscillations, at a constant calcium level, is obtained in the force-length plane when neglecting the mechanical feedback and accounting only for the cooperativity mechanism. Conversely, the force response precedes the SL oscillations, yielding a clockwise (CW) hysteresis when only the mechanical feedback is allowed to exist. In agreement with experimental observations, either CW or CCW hysteresis is obtained when both feedbacks coexist: CCW hystereses are obtained at low frequencies (<3 Hz), and the direction is reversed to CW at higher frequencies (>3 Hz). The cooperativity dominates at low frequencies and allows the muscle to adapt XB recruitment to slow changes in the loading conditions. The changeover frequency from CCW to CW hysteresis defines the velocity limit above which the muscle absorbs rather than generates energy. The hysteresis in the force-calcium relation is conveniently explained by the same cooperativity mechanism. We propose that a single cooperativity mechanism that depends on the number of strong XBs can explain the hystereses in the force-length as well as in the force-calcium relationships.

- excitation-contraction coupling
- regulated actin
- work loop
- force-length relationship

hystereses were noted to exist in the force-calcium relationship at a constant sarcomere length (SL) (13, 28) and in the force-length relationship (FLR) at constant free Ca^{2+} concentration ([Ca]) (26). Hysteresis was first observed in the force-calcium relationship in the skinned skeletal fiber (28); the force was higher during stepping down in [Ca] than during stepping up in the free calcium, at constant SL. The dependence of calcium affinity on the number of strong, force-generating cross-bridges (XBs) was suggested as the underlying mechanism (28). At the same SL and free calcium, calcium affinity is higher and XB recruitment is larger when passing through a higher force (i.e., higher calcium) rather than through a lower force. Similar hysteresis in the force-calcium relation was observed in the rat cardiac trabeculae (13). It was therefore instructive to explore the response of the force to large length oscillations, at constant free calcium, and, indeed, we have recently identified the hystereses phenomena (26) in the FLR at constant free calcium, in the intact tetanized rat cardiac trabeculae (Fig. 1). The shapes and the directions of these hystereses depended on the frequency of the SL oscillations and on the SL (26). The force response lagged the SL oscillations at low frequencies (Fig. 1*A*), leading to hystereses in the counterclockwise (CCW) direction in the force-length plane; at a given SL the force was higher during shortening than during lengthening (Fig. 1*B*). At 4 Hz and higher frequencies, the force response preceded the SL oscillations (Fig. 1*C*), and clockwise (CW) loops were observed in the force-length plane (Fig. 1*D*). The changeover frequency from CCW to CW occurred between 2 and 4 Hz (26). These different hystereses are attributed to changes in the number of strong force-generating XBs during the length oscillations. Obviously, the analysis of these hystereses should give better insight into the mechanisms underlying the regulation of XB recruitment.

Our model of the sarcomere control of contraction (19–22) couples the kinetics of calcium binding to the troponin regulatory proteins with the kinetics of XB cycling and includes two feedback mechanisms: a positive (20) and a negative one (21, 22). The positive feedback (20), denoted as the cooperativity mechanism, suggests that the number of force-generating XBs determines the affinity of troponin C for calcium. The cooperativity regulates XB recruitment and provides the explanation for the cardiac steep FLR (20) and the related Frank-Starling law at the whole heart level (19). This mechanism also provides the adaptive control mechanism whereby the loading conditions (preload and afterload) modulate the muscle energy consumption (22). The negative feedback (21), denoted as the mechanical feedback, suggests that the sarcomere shortening velocity determines the rate of XB turnover from the strong to the weak, non-force-generating state. Inclusion of mechanical feedback in our model (Fig. 2*A*) leads to the analytic solution for the experimentally established Hill's equation for the force-velocity relationship (21). The mechanical feedback also explains the linear relationship between the energy consumption and the muscle mechanical energy (22).

The combination of the two feedback loops (19, 22) allows us to describe a wide spectrum of phenomena at the cardiac fiber level as well as at the global left ventricular (LV) level (19), including the Fenn effect (22), the regulation of the cardiac relaxation (21), the length-dependent calcium sensitivity, and the changes in the free calcium transient after mechanical perturbations (21).

This study simulates the hystereses phenomena by relating to the suggested feedback mechanisms and aims to establish the dominant role of each feedback loop. The analysis highlights the effects of length perturbations and slow stepwise changes in the free calcium on the regulation of XB recruitment and force development.

## METHODS

### The Physiological Model

The physiological model and its underling assumptions were derived in our previous studies (19–22). Consequently, only the assumptions relevant to the present analysis of the effects of length oscillations on the force at a constant [Ca] (tetanic contractions), and to changes in the free Ca^{2+} at constant SL (skinned fiber), are introduced here for convenience and clarity.

The regulatory unit includes a single troponin-tropomyosin complex with 13 neighboring actin monomers and the adjacent heads of the myosin.

The XB cycle between weak and strong conformations (9) relates to nucleotide binding and release. A single ATP molecule is hydrolyzed when the XB turns from weak to strong conformation. Force is produced only by XBs at the strong conformation.

The XBs operate asynchronously of one another as assumed in other modeling work (8, 16). The active force is proportional to the number of XBs in the strong conformation, in the single overlap region. The heads of the myosin and the troponin complexes are distributed uniformly along the thick and thin filaments. The number of heads in the single-overlap region depends linearly on the length of the single-overlap region, but the fraction of XBs in the strong conformation depends on calcium kinetics and XB cycling.

The XB behaves as a Newtonian viscous element (14), as can be derived from Huxley model (16). The average force per XB (*F*_{XB}) during steady shortening is a linear function of the shortening velocity (*V*) and is given by (1)

where *F̄* is the force per XB at isometric regime (unitary force) and η is the viscous coefficient determined by the maximal shortening velocity (*V*_{U}), i.e., η = *F̄*/*V*_{U} (21, 22). We assume that *Eq. 1* applies also for sarcomere lengthening and that sarcomere lengthening increases the average force per XB (see discussion).

Calcium binding to troponin modulates XB transition from the weak to the strong conformation (5). The major role of the bound calcium is regulation of ATPase activity (5). ATP hydrolysis and phosphate release are required for the transition of the XBs from the weak to the strong conformation (9). The activity of the ATPase is inhibited in the absence of calcium (5). Thus only an insignificant number of XB transitions from the weak to the strong conformation occur without calcium binding to the neighboring troponin site.

Calcium may dissociate from troponin after the XB turns to the strong conformation. The strong XB can generate force without having calcium bound to the adjacent troponin (20, 27).

Corresponding to the mechanical feedback concept (21, 23), the rate of XB weakening (*g*) is determined by the sarcomere shortening velocity (Fig. 2*C*): (2)

where *g*_{0} represents the rate of XB weakening in isometric regime and *g*_{1} is the mechanical feedback coefficient. We assume that *Eq. 2* applies also for sarcomere lengthening and that sarcomere lengthening decreases the rate of XB weakening.

XB troponin-tropomyosin-actin interactions along the sarcomere affect the kinetics of calcium binding to the low-affinity sites on the troponin along the sarcomere (11, 12, 15). Calcium affinity *K*_{[Ca]} is determined by the number of strong XBs (*N*_{XB}) in the single-overlap region, a mechanism denoted as the cooperativity mechanism (20, 22). Only one cooperativity mechanism was included in the model. The apparent average calcium affinity was considered here, without simulating the interactions between adjacent regulatory units. The function used to describe this cooperativity is approximated at best, given our lack of knowledge. We speculate that calcium affinity is a sigmoidal function of *N*_{XB} (Fig. 2*B*), i.e., (3)

where *K*_{0}, *K*_{1}, *K*_{0.5}, and *n* are constants. *K*_{0} is calcium affinity at rest and *K*_{0} + *K*_{1} gives the maximal calcium affinity. *K*_{0} and *K*_{1} were calculated (Table 1) based on the assumption that the affinity at full activation is ∼20-fold larger than the affinity at rest (12, 19). *K*_{0.5} and *N* (*N* = 3.5) were fitted to produce the observed large phase delay (of 40°) at low frequencies. Larger *N* allows greater phase delay.

The generated force and the regulatory unit distribution are determined by two kinetics: *1*) calcium binding and dissociation from troponin, and *2*) XB cycling between the weak and the strong conformations (Fig. 2*A*). The regulatory units (*assumption 1*) are thus distributed between four different states (Fig. 2*A*). *State R* represents the rest state, in which the XBs are in the weak conformation and no calcium is bound to the troponin. Calcium binding to troponin leads to *state T*. *State A* denotes a regulatory unit “activated” by bound calcium, but the adjacent XBs are still in the weak conformation. XB turnover from the weak to the strong conformation leads to *state T*, in which calcium is bound to the troponin and the XBs are in the strong conformation. Calcium dissociation from *state T* leads to *state U*, in which the XBs are still in the strong conformation but without bound calcium.

The system variables *R*, *A*, *T*, and *U* denote the probability of having regulatory units in each state, respectively. The total density of the regulatory units, i.e., the number of regulatory units per filament unit length, is constant (4)

The system (Fig. 2) is described by three first-order, nonlinear differential equations (5) (6) (7)

where *k*_{l} and *k*_{−ℓ} represent the rate kinetics of calcium binding and dissociation from troponin, respectively. *k*_{l} is assumed to be constant (25, 29) and diffusion limited, whereas the dissociation rate, *k*_{−ℓ}, is determined by the cooperativity mechanism (*Eq. 3*), i.e. (8)

An increase in the number of strong XBs decreases the rate of calcium dissociation through the cooperativity mechanism (*Eqs. 3* and *8*). ATP hydrolysis occurs when the regulatory unit turns from *state A* to *state T* and *f* represents the kinetic rate of ATP hydrolysis (*assumption 2*). Therefore, the sarcomere energy consumption (*E*) is given by the number of regulatory units, which turn from *state A* to *T* during the twitch (9)

where *E*_{ATP} is the free energy liberated from the hydrolysis of a single ATP molecule, and TD is the entire twitch duration.

The force is determined by the number of strong regulatory units at *states T* and *U*, in the single overlap region (*L*_{s}) between the thin and thick filaments, and the average force per XB (*assumption 4*). The single overlap length is proportional to half of the SL, because the sarcomere includes two single-overlap regions (10)

where the constant *L*_{0} is determined by the actin and myosin filament lengths (20).

*Assumptions 2* and *3* lead to the number of force-generating XBs (*N*_{XB}) (11)

where the constant *N*_{C} is determined by the density of sarcomeres per muscle fiber cross-sectional area, the fiber cross-sectional area, and the number of regulatory units per half sarcomere.

The force is a product of the number of strong XBs (*N*_{XB}) and the average force per XB (*assumption 4*): (12)

where *F*_{C} = *N*_{C}*F̄* is a constant and *V*_{U} is the maximal shorting velocity (*V*_{U} = *F̄*/η).

The sarcomere shortening velocity is defined as (13)

The mathematical model was built on SIMULINK software (MATLAB). To analyze the force response to length oscillations at a constant [Ca] (0.6 μM), the SL oscillations were imposed after the tetanic contraction had reached a steady force level. Tetanus was simulated by assuming exponential rise in the free calcium to a steady level (Fig. 3). The parameters used for the simulation are summarized in Table 1. The same program and equations were also used for the simulation of the force-calcium hysteresis (13, 28) at constant SL, in the skinned fibers.

## RESULTS

The CCW and CW hystereses in the force-length plane were simulated by relating to the two feedback loops. The individual roles of the feedback loops were examined by studying the effect of each feedback loop separately. This was accomplished by eliminating one feedback loop by setting the feedback coefficients *g*_{1} = 0 or *K*_{1} = 0, and observing the effect of the remaining feedback on the phase of the force response. Thereafter, the force response to SL oscillation was studied with both feedback loops. The hysteresis in the force-calcium plane was simulated by relating only to the cooperativity mechanism.

### Assuming No Feedback Mechanisms and Unitary Force per XB

Figure 4*A* shows the force response for a simple forward system, which does not include any feedback loop, i.e., *g*_{1} = *K*_{1} = 0, and the average force per XB is constant (η = 0). As shown, the force follows the SL changes without any phase shift. Identical results are obtained for all the oscillation frequencies. In this simple case, the rate of XB weakening and the rate of calcium dissociation are constant (*g* = *g*_{0}; *k*_{−l} = *k*_{l}/*K*_{0}). Under these conditions, all of the parameters in the state equations (*Eqs. 5*–*7*) are constant and length independent. Consequently, length perturbations do not affect the distribution of the state variables and the state variables reach steady-state levels (*A* = *A*_{0} *T* = *T*_{0} *U* = *U*_{0}). With the use of *Eqs. 10* and *12*, the force is determined only by the change in the SL (14)

and there is no phase shift between the force and the length perturbations, for all the oscillation frequencies.

### Role of Mechanical Feedback

When neglecting only the cooperativity feedback (*K*_{1} = 0, in *Eq. 3*) and the XB viscous properties (η = 0), the force response to SL oscillation is determined only by the mechanical feedback. In this case, the force precedes the SL and a CW hysteresis is obtained in the force-length plane (Fig. 4*B*). CW hystereses are obtained for all tested frequencies. At any given SL, the force is lower during shortening than during lengthening, because the XB weakening rate (*g*) is faster during shortening.

An approximate analytic expression for the force response to sinusoidal length oscillations, when only the mechanical feedback exists, can be derived by reducing the four-state model (Fig. 2*A*) to two states: *W* (weak) and *S* (strong) (*W* = *R* + *A* and *S* = *T* + *U*). These two states describe the role of the mechanical feedback at full activation. Thus taking (15)

under full activation (*k*_{l}·Ca ≫ *k*_{−l}), *Eqs. 5*–*7* reduce to (16)

For sinusoidal length oscillations, SL = SL_{0} + *l*, the change in the number of strong XBs is given by (17)

and the analytic solution of *Eq. 17* yields (18)

where

The first term (*Eq. 18*) describes the transient response, which rapidly decays with a time constant of 20 ms [1/(*f* + *g*_{0})], based on Table 1.

For *t* ≫ 1/(*f* + *g*_{0}), the second term reduces to (19)

Hence the probability for having strong XBs in *state T* precedes the sinusoidal length changes. The force in this simplified two-state model is given by (20)

Utilizing *Eq. 19* for *S(t)* and assuming small length oscillations (ℓ ≪ SL_{0}) yields (21) where (22)

The force precedes the length oscillations for all oscillation frequencies because φ is positive. It is emphasized that this phenomenon does not mean that the consequence (force) precedes the cause (length perturbations), because the cause is essentially the imposed velocity of length oscillations and not the length per se. The mechanical feedback is determined by the sarcomere shortening velocity. The velocity of oscillations, the derivative of the length changes, precedes the length oscillations. Note that when we impose sinusoidal length oscillations we actually determine the oscillation velocities, at a given initial SL. At a high frequency (>3–4 Hz), the force follows the velocity, and therefore precedes the length oscillations.

At high frequencies *Eq. 21* reduces to (23)

Note that the magnitude of the force oscillation increases with the increase in the oscillation frequency (*w*/2π).

Obviously, all four states exist (Fig. 2*A*) at partial activation, but the above mathematical reduction to the two-state model does not eliminate the generality of the approximated analytic result. The transitions between *R* and *A*, or *T* and *U*, which relate to the kinetics of calcium binding to troponin, are constant and length independent when only the mechanical feedback exists. Therefore, the force response to length perturbations is determined by two factors when only the mechanical feedback exists: *1*) the effect of length changes on the single overlap length (*Eq. 10*) and *2*) the effect of shortening velocity on the rate of XB turnover (weakening) from the strong to the weak conformation. As shown in Fig. 4*A*, the first factor (changes in the overlap length) causes no phase shift, while the role of the mechanical feedback is rather evident (*Eqs. 21* and *22*) by this simplified two-state model analysis. Obviously, the same conclusion, that the force precedes the SL oscillations when only the mechanical feedback exists, holds for the four-state model too.

### Role of Cooperativity Mechanism

The effect of the cooperativity mechanism on the direction (phase) of the hysteresis was studied by neglecting the mechanical feedback loop (*g*_{1} = 0) and the XB viscous property (η = 0). As shown in Fig. 4*C*, the force lagged the SL oscillations and a CCW hysteresis was obtained in the force-length plane. The force lagged the SL oscillation and the phase delay was negative for all the tested oscillation frequencies.

An approximate expression for the phase delay is obtained by reducing the four-state model to a two-state model, *B* and *O*, which denote the regulatory units with (*B*) and without (*O*) calcium bound to troponin, respectively (*B* = *A* + *T* and *O* = *R* + *U*) (24)

*B* in this two-state model represents regulatory units associated with strong XBs. This two-state model relates only to the kinetics of calcium binding to troponin and to the cooperativity mechanism. When under constant [Ca], *Eqs. 5*–*7* are reduced to (25)

The rate of calcium dissociation (*k*_{−ℓ}) is determined by the cooperativity mechanism. No analytic expression can be derived for a sigmoidal function of calcium affinity (*Eq. 3*). Therefore, calcium affinity is assumed to be a linear function of the number of strong XBs (26)

where *k*_{0} and *k*_{1} are constant and *k*_{1}*N*_{XB} > *k*_{0} at the simulated physiological range.

The changes in *B* during sinusoidal length oscillations, SL = SL_{0} + *l*·sin(*w*·*t*), were derived by substituting *Eqs. 8*, *11*, and *26* in *Eq. 25*: (27)

Rearranging *Eq. 27* and utilizing the approximation 1/(1 + *x*) ≅ 1 − *x* for small oscillation amplitudes [*x* = ℓ/(SL_{0} − *L*_{0}) ≪ 1] yields (28)

The solution of *Eq. 28* yields the steady-state number of regulatory units at *state B*. For *t* > 1/*k*_{l}[Ca] or *t* ≫ 30 ms based on Table 1: (29)

where *B*_{SS}, the steady-state number of regulatory units associated with strong XBs when there are no oscillations, is given by (30) and the phase shift (θ) is (31)

The number of strong regulatory units (B) lags the sinusoidal length changes (φ < 0). Consequently, the force lags the sinusoidal length oscillations for all oscillation frequencies.

The amplitude of the oscillation decreases with the increase in the oscillation frequency (*Eq. 29*) and the cooperativity mechanism acts as a low-pass filter. The amplitude is inversely related to the oscillation frequency for oscillation frequencies faster than the rate of calcium binding to troponin. At low frequencies (*w* ≪ *k*_{l}[Ca]), the amplitude is apparently constant (32)

Note that the amplitude of the oscillation is inversely related to [Ca] and the SL.

### Dominant Domains of Two Feedback Loops

Consistent with the experimental data (26), either CW or CCW hysteresis is obtained in the force-SL plane when both feedback loops are allowed to exist, as shown in Fig. 5. The direction of the hysteresis depends on the oscillation frequency. A CCW hysteresis is obtained at slow oscillation frequencies (<3 Hz), where the force lags the length changes and the phase is negative (Fig. 5). Hence, the cooperativity mechanism dominates at low frequencies, as it is the only mechanism in the present model that can produce CCW hysteresis. At high oscillation frequencies (>3 Hz), the mechanical feedback dominates and the phase is positive (Fig. 5), leading to CW hysteresis, consistent with the experimental data (26). Note that the changeover frequency (3 Hz) is a function of the two feedback loops and is smaller than the rate of XB turnover from the weak to the strong stated (*f*).

### Hysteresis in Force-Calcium Plane

Stepwise changes in [Ca] have led to CCW hystereses in the force-calcium plane at constant SLs, in the skinned fibers (13, 28). The hysteresis was length dependent and diminished at longer SL (13). Note that only the cooperativity existed, because *V* = 0 in these isometric contractions. These experiments were simulated here using our model and employing only the cooperativity mechanism. The rate kinetics of XB cycling are summarized in Table 1. The coefficients *K*_{0} and *K*_{1} were reduced here, compared with the intact fiber values, because calcium affinity is lower in these skinned fiber and the maximal force is obtained at [Ca] of 10 μM (28), which is ∼10-fold larger than that required to reach maximal force in the intact fiber.

Figure 6 simulates the experimental studies in the skinned fibers at two different SLs: SL = 1.7 μm and SL = 2.2 μm. The stepwise increase in the free calcium, from “relaxing” concentration (pCa = 6.3) to partial activation (pCa = 5.6) and then to a complete activating solution (pCa = 5.0), was followed by a stepwise decrease from full activation to the same partial activation concentration (pCa = 5.6), in both SLs. As shown here, two different force levels were obtained at the same partial activation, at the short SL. The force depends on the history of contraction, and is higher after passing through a greater force level (full activation). However, this hysteresis diminishes at longer SL and identical force levels were obtained at partial activation, at the longer SL.

## DISCUSSION

The FLR is an important constitutive feature of the cardiac muscle (1). As shown here, a unique FLR can be obtained by a simple model of the sarcomere that does not include feedback loops, where calcium affinity and the rates of XB cycling are taken as constants (Fig. 4*A*). In this case, the model's state variables *A*, *T*, and *U* are constant and independent of the history of contraction (no memory!), for a constant [Ca]. The SL changes affect the force only by changing the single overlap length, and a unique FLR is obtained.

Hystereses are observed in the force response to length oscillations when the feedback loops that regulate calcium affinity and XB weakening rate are included in the analysis. The hysteresis phenomenon is characteristic of control systems with memory (6), and the memory in the sarcomere resides in the state variables *R*, *A*, *T*, and *U*. The length perturbations affect the distribution between these states, through the two feedback loops: the cooperativity mechanism and the mechanical feedback, which have opposite effects on the direction (phase) of the hysteresis. The advantage of this theoretical study is that it allows us to identify the working mechanisms of sarcomere control by separately testing the role of each feedback loop. Obviously, it is difficult, if not impossible, in experimental approaches to isolate and study the effects of any particular mechanism by suppressing the effects of all other mechanisms.

The study reveals that the force lags the SL oscillations when only the cooperativity mechanism is allowed to exist, leading to CCW hystereses in the force-length plane, where the force is greater during shortening than during lengthening at the same given SL. The mechanical feedback, which determines the rate of XB weakening, generates only CW hystereses. When both feedbacks act simultaneously in the system, the direction (phase) of the hysteresis depends on the frequency of length oscillations due to the opposite effects of the two feedback loops. The cooperativity mechanism acts like a low-pass filter: the amplitude of the oscillatory force response derived from the cooperativity feedback is maximal at low frequency and decreases with the increase in the oscillation frequency (*Eq. 29*). The mechanical feedback and the XB viscous property act as a high-pass filter (*Eq.* 23). The amplitude of CW oscillations derived from the mechanical feedback (*Eq. 23*) and the effect of the XB viscous property (*Eq. 1*) decrease at slower oscillation frequencies. Consequently, at slow oscillation frequencies (bellow 3Hz) the cooperativity mechanism dominates and CCW hystereses are obtained, whereas at higher frequencies the increased velocity leads to CW hystereses. The results are consistent with the experimental data (26) presented in Fig. 1.

The existence of a CCW hysteresis in force-calcium relation at constant SL (13, 28) as well as the CCW hysteresis obtained in the FLR at constant free calcium and low frequencies (26) can be explained by the same cooperativity mechanism, whereby calcium affinity depends on the number of strong XBs (28); longer SL or higher [Ca] induces more strong XBs and higher calcium affinity. Therefore, a larger force is obtained when reaching a given length and a [Ca] from a longer SL or a higher [Ca]. Note that only CCW hystereses are observed in the force-calcium plane, at constant SL, because under these conditions, the shortening velocity is zero and only the cooperativity mechanism prevails.

The hystereses between the SL (input) and the force (output) in the present model result from a phase shift between the state variables and the SL. However, the analysis of the quasi-steady-state hysteresis in the force-calcium plane reveals that the cooperativity mechanism turns the system to a nonlinear system because the rate of calcium dissociation, *K*_{−l} in *Eqs. 5*–*7*, is not constant but depends on state variables *T* and *U*, through the cooperativity mechanism.

It is evident that only the suggested cooperativity mechanism, which is determined by the number of strong XBs, can explain the existence of several force levels at a given SL and a given [Ca]. Other conceivable mechanisms that are governed by the SL (1) or by [Ca] (11) yield the same steady-state force at a given SL and free calcium values without a hysteresis. The steady-state force in the four-state model (Fig. 2) is given by (20) (33)

Assuming that *k*_{l}[Ca] ≫ *g*_{0} for the skinned fibers reduces *Eq. 33* to (34) where

and *C*_{I} = *g*_{0}/(*f* + *g*_{0}) are constants at a given SL. Calcium affinity, however, is a function of the number of strong XBs (*Eq. 3*). Substituting *Eq. 3* and rearranging *Eq. 34* yields the following polynomial (35)

where *a*_{0}, *a*_{1}, *a*_{n}, and *a*_{n+1} are constants that are determined by the cooperativity coefficients (*K*_{0}, *K*_{1}, *K*_{0.5}, and *n* from *Eq. 3*), the SL and [Ca] (36)

The steady-state force can be derived from the polynomial (*Eq. 35*), which is one order larger (*n* + 1) than the order of the equation for the cooperativity mechanism (*Eq. 3*). The analysis yields that more than one steady force level can be obtained at the same SL and same free calcium, if the order of the cooperativity mechanism (*n* in *Eq. 3*) is ≥1. These different force levels can be reached from different initial conditions. This does not imply that several force levels always exist. For a given set of parameters, SL and [Ca] (*Eq. 36*), only one (physical) solution exists and can be reached, whereas for the other set of parameters, SL and [Ca], two distinct physical solutions exist (13, 28). Thus the same system can explain the observed significant hysteresis at short SL and the insignificant or no hysteresis at longer SL. This is simulated in Fig. 6, and is consistent with the experimental data (13).

The term hysteresis is used in various fields of engineering and science and various models were suggested to simulate the hysteresis phenomenon. Sain et al. (31) described four phenomenological models of hysteresis between input and output signals and suggested that the basis for the hysteresis phenomena is the existence of two or more switching thresholds. Cook (6) has defined hysteresis as a nonlinearity, where the output can get two or more values for a given input so that the output value depends on the history of the input as well as on its current value. This general definition does not require any threshold point. The present mechanistic model does not include any threshold but still generates hystereses in the force-length and the force-calcium planes. Similarly, viscoelastic substances produce CW hysteresis in the force-length plane without having any threshold point. A physiological textbook (30) defines hysteresis as “a system which fails to follow identical paths on application and withdrawal of a forcing agent exhibits hysteresis.” Our observations (26) and the present simulations are consistent with these (6, 30) classic definitions.

The strength of the present mechanistic four-state model over other phenomenological models (31) is based on two significant features. First, the present model compresses a broad spectrum of physiological data (21–23) into a set of theoretical assumptions and rules that include only two kinetics and two feedback loops. The model provides an explanation to basic key features of the cardiac muscle, including the force-velocity relationship (21) and the regulation of energy consumption (22). Second, the same set of relatively few parameters (Table 1) describes both the CW and CCW hystereses, whereas phenomenological models (31) require different sets of parameters to describe the various observations.

Small-amplitude sinusoidal length oscillations (18, 32, 34), of the order of 0.25% of the muscle length, have been used to study the muscle dynamic stiffness in the presence of saturated [Ca] or barium. These studies (18, 32, 34) were used mainly to characterize the single XB dynamics and the kinetics of XB turnover between the various biochemical conformations. Campbell et al. (4) have presented a model that couples the actin on-off kinetics with a three-state model of XB kinetics to explain the effect of small-amplitude oscillations on the dynamic stiffness. Campbell et al. have shown that XB kinetics alone cannot explain the observed negative phase in the dynamic stiffness, where the stiffness lags behind the length changes, and suggest that such behavior should be attributed to some cooperativity mechanism. The present study focuses on the regulation of force development (and not the stiffness) and relates to an experimental study that utilized large amplitudes (26) of length oscillations. These large length perturbations were well above the single XB stroke step (of ∼5 nm) and therefore the observations and the simulated hystereses relate to the regulation of XB recruitment and the number of strong XBs. Campbell et al. (4) have used the linearization technique to describe the effect of small length oscillation (<1%). However, the system is highly nonlinear and the observations described here were obtained at 10-fold larger amplitude than those used (4) for the measurement of the dynamic stiffness. The present model simulates the large changes in the SL and free calcium in this nonlinear system, without linearization or limitation to small changes.

In addition, Campbell et al. (4) have suggested that the SL has no effect on the phase of the stiffness. However, the SL showed significant effect when large length oscillations were studied (26). Maximal phase shift was observed at small SL (1.75 μm) and the phase of the CCW hysteresis decreased at longer SLs. The present model clearly explains the effect of the SL. The effect of the cooperativity mechanism decreases at longer SLs (*Eq. 32*), whereas the effects of the mechanical feedback (*Eq. 19*) or the viscous property of the XBs (8) depend on the amplitude of length oscillation but not on the absolute SL value. Hence the cooperativity dominates at short SLs and its effect decreases at longer SLs.

The CCW hysteresis described here differs from the so-called “work loop” obtained at physiological twitch contractions (3, 7). Different phase shifts during twitch contractions, which were observed between the force and the imposed muscle length oscillations, were also denoted as hystereses (10, 24). However, these hystereses at the force-length plane resulted mainly from changes in the activation level during the contractions. The term work loops (3, 7) is more precise for these observations during twitch contractions, as these loops describe the work done during the twitch contraction. The hystereses studied here in the force-length plane were obtained at constant activation level, utilizing tetanic stimulation. The term hysteresis is used here for the length oscillations when, for a given length (input), there are two or more finite force values (output) while all the other controllable contractile parameters, including the activation level, are constant.

The mechanical feedback (21) and the XB viscous property (8) produce only CW hysteresis in the present model. The sarcomere shortening velocity affects the number of strong XBs, through the mechanical feedback (*Eq. 2*), and the average force per XB (*Eq. 1*). During sarcomere shortening (*V* > 0), the rate of XB weakening increases and the force decreases. The opposite holds during lengthening (*V* < 0), and the force is consequently larger during lengthening than during shortening at the same SL, again due to the mechanical feedback. The XB viscous property that determines the average force per XB provides a similar effect: the force per XB is larger during lengthening than during shortening. The differentiation between the role of the mechanical feedback and the XB viscous property is of significant importance. Measurements of the dynamic stiffness, which is proportional to the number of strong XBs, allow distinguishing between the two. This required superimposing small amplitude (∼3–5 nm) and high-frequency oscillations (200 Hz) over the large length oscillations used for quantifying the feedback mechanisms. We have found (26) that the dynamic stiffness, as measured from the amplitude of the high-frequency oscillations (200 Hz), is modulated by the low-frequency oscillations; similar phase delays appeared in the force-length and stiffness-length planes (Figs. 2–4 in Ref. 26). CCW hystereses were obtained in the force-length and in the stiffness-length planes at low frequencies, and CW hystereses were observed in both planes at high frequencies. These observations suggest that the feedback mechanisms rather than the viscous properties play the key role in determining the hystereses phenomena.

The area within the force-length loop represents the total work done. The CCW direction implies that the muscle does work upon its load while a CW direction infers that work is actually done on the muscle. CCW hysteresis can be explained only by the existence of more XBs during the shortening phase than during the lengthening phase of the SL oscillations. Therefore, CCW hysteresis must be attributed to a cooperativity mechanism that regulates XB recruitment. The XB viscous property and the mechanical feedback produce CW hystereses that resemble the common appearance of CW hystereses in various physical systems, e.g., viscoelastic and magnetic systems, where the energy dissipates as heat during the cycle. By comparison, the muscle behaves as a shock absorber and absorbs external energy at high frequencies, where the mechanical feedback dominates.

The present model suggests the existence of one dominant cooperativity mechanism whereby calcium affinity is determined by the number of strong XBs (*assumption 8*). Other types of cooperativity that depend on the SL (1) or [Ca] have also been suggested (11). However, length-dependent mechanisms (1) cannot explain a hysteresis at constant SL (13, 28) and calcium-dependent mechanisms (11) cannot explain the hysteresis in the FLR at constant free calcium. Only the suggested cooperativity mechanism (20) can explain the hystereses in the FLR as well as the force-calcium relationships.

While the existence of more than one cooperativity mechanism cannot be ruled out, the dominance of the cooperativity mechanism suggested here is supported by additional observations: it explains the “length-dependent calcium sensitivity” (20), the steep cardiac FLR, the effect of loading conditions and length perturbations on calcium transient (20, 21) and the adaptive control of energy consumption (22) whereby the loading conditions (preload and afterload) affect XB recruitment and energy consumption. However, despite the ability of the present model to describe a wide spectrum of phenomena at the cardiac fiber level (19–21, 23) as well as the global LV level (20), it does have some limitations.

It is assumed here that the average force per XB is a linear function of the shortening velocity, at steady state, but the shortening velocity is not at steady state during sinusoidal length oscillations. However, it can be assumed that for most practical slow physiological contractions and the slow oscillations studied here, the average force per XB is in a quasi-steady state because the dynamics of myosin head conformational changes are by far faster than the length perturbations. The force per XB is determined by conformational changes in the myosin head and has a characteristic time response of <0.25 ms, as observed in the T_{1}-T_{2} studies (17). The average force per XB reaches the steady level within few milliseconds (8, 17). The model simulates the slow transient changes in the force response to slow length oscillations. These changes are determined by the slow kinetics involved in the regulation of XB recruitment and XB weakening (rate constants of 10–40 s^{−1}, relating to *g*_{0} or *f*, respectively). Thus the fast kinetic that determines the average force per XB is practically at steady state during the slow sinusoidal length perturbations.

The exact form of the cooperativity function is still unknown. It was assumed that calcium affinity depends on the number of XBs (*Eq. 3*). Yet it is still unclear whether the affinity depends on the number of XBs or on the force. It is impossible to determine this point based on our present experimental results, and a different experiment has to be designed to solve this enigma. However, the slow transient changes in [Ca] in response to fast length perturbation (2) suggest that calcium affinity is indeed determined by the number of strong XBs and not by the force or the number of attached XBs.

The model includes 12 parameters (Table 1), and their exact values are not known. Also, [Ca] during the steady tetanus was not measured (26) and the exact function of the cooperativity mechanism is unknown. These severe limitations dictate that we can only explain the phenomena in a qualitative manner. Clearly, the ability to simulate the phenomena is the first step in the pursuit of better understanding. The aims of the present study were to define the effects of the suggested feedback loops on the direction of the hystereses and to examine whether we can explain the experimentally observed changeover from CCW to CW hysteresis. Future studies, wherein the cooperativity function will be better defined and include the measurement of [Ca], are clearly indicated.

The present simulation assumes that [Ca] is constant during the imposed length oscillation. However, the simulated changes in calcium affinity are expected to produce oscillations in [Ca]. This phenomenon was ignored for the simplicity of the simulation.

Our model assumes that the rate of XB weakening is a linear function of the velocity, for both sarcomere shortening and lengthening. The mechanical feedback plays the key role during sarcomere lengthening and lengthening decreases the rate of XB weakening. This hypothesis is consistent with our observations (26) that the dynamic stiffness increases with sarcomere lengthening.

The above-mentioned limitations do not impose a major constraint on our hypotheses, although they present new questions and challenges for further research, e.g., what is the real nature of the cooperativity mechanism? Does the mechanical feedback apply to sarcomere lengthening too, with the same function and same coefficient?

In conclusion, the study identifies the role of the two feedback loops that modulate sarcomere contraction. The cooperativity mechanism that regulates calcium affinity and XB recruitment produces only CCW hystereses. The mechanical feedback that regulates the rate of XB weakening generates only CW hysteresis. The phase of the actual hysteresis depends on the frequency of the oscillations and the interplay between these two feedback loops. At low frequencies, the cooperativity mechanism dominates. As the frequency increases, the shortening velocity increases and the mechanical feedback becomes prominent. The changeover frequency form CCW to CW hystereses defines the maximal rate at which the muscle can adapt to the changes in the prevailing loading conditions by modulating XB recruitment. Whereas the existence of other cooperativity mechanisms cannot be completely ruled out, it seems that the suggested cooperativity mechanism, whereby calcium affinity is a function of the number of strong XBs, can explain the existence of hystereses in both the force-calcium and force-length planes. The analysis of these phenomena strengthens our insight on the intercellular control of contraction, and leads to a better understanding of the cardiac function.

## GRANTS

This study was partially supported by a grant from the United States–Israel Binational Foundation and by the Fund for Promotion of Research at the Technion (to A. Landesberg).

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society