## Abstract

To integrate myocardial contractile processes into left ventricular (LV) function, a mathematical model was built. Muscle fiber force was set equal to the product of stiffness and elastic distortion of stiffness elements, i.e., force-bearing cross bridges (XB). Stiffness dynamics arose from recruitment of XB according to the kinetics of myofilament activation and fiber-length changes. Elastic distortion dynamics arose from XB cycling and the rate-of-change of fiber length. Muscle fiber stiffness and distortion dynamics were transformed into LV chamber elastance and volumetric distortion dynamics. LV pressure equaled the product of chamber elastance and volumetric distortion, just as muscle-fiber force equaled the product of muscle-fiber stiffness and lineal elastic distortion. Model validation was in terms of its ability to reproduce cycle-time-dependent LV pressure response, ΔP(*t*), to incremental step-like volume changes, ΔV, in the isolated rat heart. All ΔP(*t*), regardless of the time in the cycle at which ΔP(*t*) was elicited, consisted of three phases: *phase 1*, concurrent with the leading edge of ΔV; *phase 2*, a brief transient recovery from phase 1; and *phase 3*, sustained for the duration of systole. Each phase varied with the time in the cycle at which ΔP(*t*) was elicited. When the model was fit to the data, cooperative activation was required to sustain systole for longer periods than was possible with Ca^{2+} activation alone. The model successfully reproduced all major features of the measured ΔP(*t*) responses, and thus serves as a credible indicator of the role of underlying contractile processes in LV function.

- myosin cross bridge
- rat
- Recruitment-Distortion-Instantaneous-Elastance model

cardiovascular scientists have long been interested in connecting heart function with the regulated action of molecular motors that gives rise to that function. Despite an effort in the 1960s and 70s to build a quantitative bridge between whole-heart function and muscle-contraction indices (6), a satisfactory connection between molecular motors and heart function remains an unrealized goal (26). One approach for making this connection uses complimentary experiments in preparations, representing several organizational levels in the cardiac hierarchy (56), i.e., isolated myofilament, myocyte, muscle fiber, and whole heart. Thus, the figment of a myofilament effect may be traced through the several levels of organization until finally it is seen at the level of heart function. Another approach is an international effort to synthesize a comprehensive heart model out of all of the heart's detailed and distributed molecular, cellular, tissue, valve, and fluid complexities (27). When completed, and when a sufficiently powerful computer is found to run such a model, the functional consequences of change in regulation of a molecular motor may be predicted. Alternative to the comprehensive model is the effort to connect regulated motor kinetics with heart function through lumped approximations and well-considered simplifications in reduced-order cardiac models (9, 11, 32, 37). We argue that an appropriate reduced-order cardiac model will be of value both in tying results together from the experimental approach, where different preparations are used, and in providing a simple expression of the underlying complexity to guide the comprehensive modeling approach. The work described in this report represents further effort toward realization of a satisfactory reduced-order heart model.

Left ventricular (LV) pressure and volume are the fundamental variables of heart function, and some variant of the pressure-volume relationship has been the mainstay in representing heart function since the earliest descriptions (23, 46). Whereas early work made use of end-state or cycle-average forms of the pressure-volume variables, Suga, Sagawa, and colleagues (40, 41, 51, 52) emphasized the time course of pressure-volume ratio over the heart cycle. Their results took the form of linear pressure-volume isochrones that changed progressively during the heart beat. These experimental findings legitimized representing the LV in terms of a time-varying elastance, a concept employed previously by early cardiovascular modelers (20, 57). Subsequently, time-varying elastance has been a cornerstone for LV representation in a large body of cardiovascular modeling and simulation.

However, by itself, time-varying elastance does not well reproduce instantaneous pressure-volume-flow events throughout ejection and relaxation. Thus many attempts were made to add resistance and deactivation elements to the time-varying elastance core of the model to improve predictive performance (5, 12, 15, 29, 44, 45, 52, 53, 55, 58). Although some improvements were realized, all such models generated predictions with ever-increasing errors as the cycle advanced into late systole and relaxation (7, 12). Generally, the predictive power of time-varying elastance-based models is limited to the data from which the model parameters are estimated. Thus these models take on an ad hoc character, which severely limits their general applicability.

To improve on time-varying, elastance-based LV simulations, information is needed that reveals how underlying causes of time-variation in LV properties (i.e., waxing and waning myocardial Ca^{2+} activation) may be distinguished from the dynamic effects of volume change in influencing pressure generation. To this end, Burkhoff (9) and Landesberg (32) have constructed reduced-order models that promise to couple the dynamics of LV pressure-volume relations with the kinetics of underlying contractile processes. In fact, several groups are now using measurements of global LV Ca^{2+} transients and the Burkhoff model (3, 9, 35, 38, 42, 50) to link Ca^{2+}-activated myocardial contraction with chamber pressure-volume characteristics.

However, existing model formulations do not consider the dynamics of volume-driven LV pressure responses interacting with the time-varying effect imposed by waxing and waning myocardial Ca^{2+}. To illustrate, Hunter et. al. (28, 29) described the pressure response to step-like volume changes in the isolated dog heart and how these pressure responses varied when volume changes were delivered at different times in the heart cycle. Any single response evolved according to dynamic processes initiated by the volume step and to the modulation of these dynamic processes by time-varying Ca^{2+} activation. This interaction between dynamic volume effects and Ca^{2+} effects was different depending on when in the cycle the volume change was delivered. None of the existing models can account for how time variation in Ca^{2+} activation interacts with volume-determined pressure effects to establish LV pressure at all instants in the cycle. There continues to be a need for a more comprehensive reduced-order cardiac model; one that is based on a complete accounting for the interaction between volume and Ca^{2+} activation as seen in the cycle-time-dependent pressure response.

In this study, we use the volume-perturbation methods of Hunter et al. (28, 29) to obtain a comprehensive dynamic signature of the isolated rat heart in terms of a family of cycle-time-dependent pressure responses, ΔP(*t*), to step-like incremental volume changes. We then build a model of LV pressure-volume dynamics, starting with myofilament contraction processes and the apparent physical properties of the chamber. We then use the family of ΔP(*t*) responses to elucidate parametric forms of kinetic coefficients. The result is a model that satisfactorily reproduces all important aspects of the ΔP(*t*) responses. The proposed model adds a dynamic dimension to our understanding of LV pressure-volume relationships and provides new perspective for prediction and interpretation.

## MODELING METHODS

We built a reduced-order LV model in two stages: a priori construction of essential model structure and a posteriori assignment of parametric forms to nonlinear kinetic coefficients based on experimental results.

### A Priori Model Construction

#### Stiffness is the unifying physical concept that connects molecular motors to LV chamber.

The myosin cross bridge (XB) is the molecular motor that drives LV function. The elastic segment of the molecule is stretched as part of the chemomechanical transduction that occurs when the head of the molecule undergoes conformations change as it binds to an actin attachment site. As a result, the XB generates a force, f(*t*), that equals the elastic distortion, *x*(*t*), times the molecular stiffness, ε_{0}. Thus, (1) A collection of parallel XB in a muscle sarcomere, or in a muscle fiber, produces a force, F(*t*), that is equal to the number, *N*(*t*), of parallel, force-bearing XBs times the average force generated by a single XB. Thus, (2) where *x*(*t*) is now the average distortion among the *N*(*t*) XBs. The overall stiffness of this fiber, ε(*t*), is given by (3) This allows Eq. 2 to be rewritten as (4) Thus muscle-fiber force equals a XB-based fiber stiffness times a XB-referenced lineal elastic distortion.

The issue now is how muscle-fiber physical properties of stiffness and elastic distortion become expressed in the physical properties of the LV chamber. A thick-walled sphere is the simplest geometric approximation of the LV of small mammals (Fig. 1). We consider that midwall, muscle fibers are circumferentially arranged and serially connected along the midwall circumference. The centrality of these circumferential midwall fibers in the overall fiber arrangement allows them to be taken as representative of all muscle fibers in the heart. This is a reasonable assumption providing that the distribution of material and physiological properties satisfy the condition of essential homogeneity. Support for this assumption comes from reports of uniform sarcomere shortening patterns throughout the myocardium (1, 22, 25, 39).

Because wall volume and wall mass are related by a density factor that is ∼1 g/ml, wall mass is equal to wall volume, V_{w}. Then the midmass circumference, *L*(*t*), may be taken as the length of a circumferentially-oriented midwall fiber, and this length is related to chamber volume, V(*t*), by (5)

The value of *L* at a baseline volume, V_{BL}, is designated *L*_{BL}. For small volume variations, ΔV, around V_{BL}, the resulting incremental change in *L*, Δ*L*, is given by the proportionate relation (6) where the multiplying factor on the right-hand side is the partial derivative of *L* with respect to V in Eq. 5 evaluated at V_{BL}.

Importantly, the transformation given by Eq. 6 applies to all conversions of lineal dimensions along the circumference into analogous volumetric dimensions. For instance, the elastic distortion in a muscle fiber corresponding to *x*(*t*) in Eq. 4 may be converted into its volumetric equivalent, υ(*t*), in the LV chamber by the equation (7)

The circumferential midwall fiber may be visualized as a torus with small cross-sectional area. The result is that the force produced by the fiber, F(*t*), is best expressed in units of force per unit cross-sectional area. Thus F(*t*) is the stress at the midwall and, as a result of half-shell force balance, (8)

In response to the Δ*L*, the fiber will generate an incremental change in F, ΔF, which results in an incremental change in P, ΔP. This is given by the proportionate relation (9) where the constant multiplying factor within the brackets on the right-hand side is the partial derivative of P with respect to F in Eq. 8 evaluated at V_{BL}.

By means of the transformations of Eqs. 6 and 9, it can be shown that the relation between mid-wall muscle stiffness, ε = ΔF/Δ*L*, and LV chamber elastance, E = ΔP/ΔV, is (10) where the multiplying term on the right-hand side is again a constant.

These equations demonstrate that for small volume perturbations around a baseline volume, muscle fiber stiffness is linearly related to LV chamber elastance, and elastic distortion within the muscle fiber is linearly related to volumetric distortion, υ(*t*), of the LV chamber.

The nontrivial outcome of these linear transformations is that it allows writing an expression for LV chamber pressure in terms of chamber elastance that is analogous to Eq. 4 for muscle fiber force in terms of muscle stiffness, i.e., (11)

#### Kinetics of myofilament activation give rise to differential equation for muscle-fiber stiffness.

In Eq. 4, ε(*t*) changes as XB are recruited into and out of the *N*(*t*) pool (see Eq. 3). This recruitment is a consequence of both myofilament activation and fiber-length-induced effects. Myofilament activation occurs as a result of Ca^{2+} binding to troponin C (Tn-C) of the regulatory complex on the thin filament and as a result of cooperative activation from the formation of strong-binding XB. We equate the number of force-bearing XB with the number of activated binding sites on the thin filament. A kinetic scheme for the activation state of thin-filament XB binding sites is given in Fig. 2. These sites are either not permissive for XB attachment (*N*_{off}, not activated) or permissive for XB attachment (*N*, activated). A thin filament site may have Ca^{2+} bound (indicated by * in Fig. 2, *left*) or not (indicated by ° in Fig. 2, *right*). If Ca^{2+} is bound, thin-filament site activation takes place readily according to the activation constant *k*_{on}. Activation may also occur as a result of cooperative XB action according to the coefficient *k*_{xb}(*t*). Loss of activation occurs according to the kinetic constant *k*_{off}.

A reduced-order equation for the four-state kinetic scheme in Fig. 2 may be written by adopting a Ca^{2+}-binding equilibrium assumption and applying a conservation constraint. Approximate equilibrium in Ca^{2+} binding to Tn-C exists if this binding is fast relative to the other kinetic steps. This equilibrium is characterized by the equilibrium coefficient, *K*_{eq}. Thus, Ca^{2+} bound to Tn-C in the presence of a myoplasmic activator calcium signal, Ca(*t*), is given by *K*_{eq}Ca(*t*). With this equilibrium assumption, the Ca^{2+}-bound and non-Ca^{2+}-bound states need not be treated independently. Accordingly, the four states in Fig. 2 may be reduced to two states: *1*) sites that are not activated (represented by *N*_{off} on the top of Fig. 2; *N*_{off} = *N*_{off}* + *N*_{off}°); and *2*) sites that are activated and to which XB have attached (represented by *N* on the bottom of Fig. 2; *N* = *N** + *N*°). In both states, the Ca^{2+}-bound and non-Ca^{2+}-bound conditions (*N** and *N*°, respectively) are related to the sum of Ca^{2+}-bound and nonbound conditions, *N*, by (12) The equilibrium condition allows economy in expression by defining the net outcome of the Ca^{2+} activation processes in terms of a time-varying coefficient, α(*t*): (13) Consequently, all the effects of Ca^{2+} activation appear in the model in the form of the α(*t*) coefficient. Further economy is realized by formulating a single coefficient governing all of activation, *on*(*t*), to be the sum of that due to Ca^{2+} activation plus that due to XB cooperative activation (14) Similarly, a coefficient governing deactivation, *off*(*t*), may be formulated as (15)

Simplifying from the conservation constraint comes from noting that the total number of thin-filament sites, *N*_{T}, equals the sum of sites in all states, (16)

Fiber length, *L*(*t*), now enters the recruitment equation when it is noted that on the ascending limb of the force-length curve, the total number of thin-filament sites that can participate in the contraction process, *N*_{T}, increases with fiber length. This relation may be approximated by (17) where *b* is the slope of the tangent line to the *N*_{T} vs. *L* relationship at a reference fiber length, and *L*_{0} is the extrapolated length-axis intercept of this tangent line.

Assuming a stoichiometry of one XB per activated thin filament site, noting the relation between ε(*t*) and *N*(*t*) in Eq. 3, and employing the above simplifications results in a simple one-equation formulation for myofiber stiffness in terms of XB recruitment kinetics. (18)

Importantly, ε(*t*) is an instantaneous stiffness, i.e., the muscle-fiber stiffness encountered when a quick stretch is applied. The differential equation character of Eq. 18 imparts a dynamic dimension to ε(*t*) as well as showing how ε(*t*) changes come about in response to changes in both fiber length and myofilament activation.

#### Kinetics of XB detachment gives rise to differential equation for muscle-fiber elastic distortion.

An equation for the effects of muscle fiber length on average distortion of the elastic XB, *x*(*t*), is based on XB cycling kinetics and has been developed in previous publications (11, 14, 16–18). Briefly, as XB attach to the thin filament during isometric contraction, they undergo a power stroke that imposes *x*_{0} distortion on the elastic segment of the XB, and consequently force develops. Changes in muscle length cause filament sliding, which leads to further elastic distortion of these attached, force-generating XB over and above that imposed by the XB power stroke. However, because XB lose their distortion when they detach and because when they reattach they do so with *x*_{0} distortion, the differential equation describing changes in XB distortion with changes in muscle length is, to an approximation, (19) where *x*_{0} is XB distortion during isometric conditions, and *g* is the XB detachment rate constant. Thus, *L*(*t*) affects *x*(*t*) only through its first-time derivative. Without a rate-of-change in *L*(*t*), *x*(*t*) takes on its isometric value *x*_{0} regardless of the fiber length.

#### Dynamic relations between LV pressure and volume can be derived from the dynamic relations between muscle-fiber force and length.

The muscle-fiber differential equations (Eqs. 18 and 19) were adopted to describe the dynamic behavior of LV elastance and LV volumetric distortion because of the proportionate relation between: *1*) Δ*L* and ΔV, Eq. 6; *2*) ΔP and ΔF, Eq. 9; *3*) υ(*t*) and *x*(*t*), Eq. 7; and *4*) ε(*t*) and E(*t*), Eq. 10. Thus, (20) (21) By analogy with XB recruitment, *on*(*t*) and *off*(*t*) represent kinetic coefficients governing the formation and dissolution of the elastance state, E(*t*), from elements in a nonelastance state; β[V(*t*) − V_{0}] represents the ascending limb of the Frank-Starling mechanism; *g* represents a coefficient by which volume-induced changes in υ(*t*) are dissipated; and υ_{0} represents a baseline volumetric distortion imposed during the formation of the elastance state during an isovolumic beat. Time variation in *on*(*t*) and *off*(*t*) accrues primarily from a coefficient, α(*t*), which is related to the time course of Ca^{2+} binding to the thin filament (Eq. 13). By analogy with ε(*t*), E(*t*) is an instantaneous elastance, i.e., the elastance exhibited during an instantaneous volume change.

Instantaneous LV elastance, E(*t*), appears to be composed of elastance elements that have been recruited into an elastance state, E, from a nonelastance state (E_{0}, Fig. 3). There is a total number of potential elastance elements, E_{T}, which is the sum of elements in E and E_{0} states, i.e., E_{T} = E_{0} + E. Also, E_{T} is a function of chamber volume just as the total number of potential force-bearing XB is a function of muscle length. In general this function is nonlinear, but in the case of incremental volume variations around a baseline volume, V_{BL}, the nonlinear character of the function cannot be discerned. Rather, the slope of this nonlinear function, β, around V_{BL} is observed, and it enters the relationship as (22) where V_{0} is an extrapolated volume axis intercept of the tangent line to the nonlinear function at V_{BL}.

Two additional note-worthy items: *1*) the time dependence of the [*off*(*t*) + *on*(*t*)] coefficient signifies that the LV model is time-varying in accord with the common understanding of the time-varying nature of the LV; and *2*) the two forcing inputs on the far right-hand-side of Eq. 20, i.e., *on*(*t*) and V(*t*), combine as a product rather than as a sum and, thus, neither *on*(*t*) nor [V(*t*) − V_{0}] will cause E(*t*) if the other is zero.

### A Posteriori Parameter Assignment

The model described above was developed prior to ΔP(*t*) data collection. Whereas this a priori construct served as the structural foundation of the model, key nonlinear features were introduced into model kinetic coefficients a posteriori to adequately reproduce three ΔP(*t*) response behaviors. These behavioral features are described in detail in the results section and, in brief, are: *1*) ΔP(*t*) to ΔV delivered prior to beat onset exhibited a sustained time course that extended well into late systole (Figs. 6–⇓8); *2*) ΔP(*t*) to ΔV delivered in early systole was exaggerated relative to ΔP(*t*) to ΔV delivered prior to beat onset (Fig. 8); and *3*) ΔP(*t*) to positive and negative ΔV were detectably asymmetric about the time axis (Fig. 9). The cogent nonlinearities to reproduce these features were as follows.

#### Cooperative activation accounted for extended ΔP(t) time-course.

Cooperativity is a mechanism whereby the intensity of contraction facilitates the development of greater contraction intensity. In the current model, we placed all cooperative effects in the activation coefficient, *k*_{xb}(*t*). If not well formulated in a mathematical sense, cooperativity can drive model contraction process toward progressively greater contraction, and this will not allow relaxation to occur. After attempting several versions, we found that the following was the best among a limited set of competing formulations for representing cooperative recruitment: (23) where *k*_{a} is a scaling factor, E_{50} is the value of E(*t*) at which *k*_{xb}(*t*) achieves half its potential maximum, and *n* is an exponent that shapes the time course of *k*_{xb}(*t*). Thus, *k*_{xb}(*t*) fluctuates up and down as E(*t*) rises and falls, but it does so with a distinctly different wave shape than that of E(*t*). Importantly, when *k*_{xb}(*t*) is formulated according to Eq. 23 (and *k*_{a}, E_{50}, and *n* are appropriately chosen), E(*t*), as an outcome of the solution of Eq. 20, asymptotically approaches zero when α(*t*) goes to zero. Because *k*_{a}, E_{50}, and *n* are all fixed values, assignment of these allow *k*_{xb}(*t*) to be estimated, providing a value of E(*t*) or a suitable surrogate is available.

#### Stretch-activation accounts for exaggerated ΔP(t) responses to early-cycle ΔV.

It is well known that cardiac muscle stretch contributes to its activation (47–49). For the midwall fiber, sudden stretch has two dimensions: *1*) a permanent increase in fiber length, *L*(*t*), and *2*) a transient increase in XB distortion, *x*(*t*). At the level of the LV, sudden volume inflation produces a permanent increase in V(*t*) and a transient increase in υ(*t*). One mechanism of stretch activation produces a sustained effect and works through cooperative activation as follows. The increase in V(*t*) causes an increase in E_{T} (Eq. 22), resulting in more E(*t*) (Eq. 20), which increases *k*_{xb}(*t*) (Eq. 23), leading to the formation of even more E(*t*). A second mechanism produces a transient effect that is linked to the distortion variable, υ(*t*). We represented these distortion effects by a fractional distortion, φ(*t*) where (24) During isovolumic beating, υ(*t*) never deviates from υ_{0}, and φ(*t*) equals zero throughout the beat. During rapid volume inflation, φ(*t*) increases precipitously with volume inflation and then dies away as υ(*t*) − υ_{0} approaches zero, which occurs in accord with Eq. 21.

Thus a transient effect of stretch on activation was obtained when *on*(*t*) was defined as (25) where ρ is a coefficient that grades the dependence of the recruitment coefficient, *on*(*t*), on the distortion fraction. The effect of ρφ(*t*) is to cause transient stretch activation when φ(*t*) > 0 and to cause transient shortening deactivation when φ(*t*) < 0.

#### Distortion-induced loss of elastance elements accounts for diminished ΔP(t) responses to late-cycle ΔV and for response asymmetry.

Just as force-bearing XB are lost from the force-generating pool when XB strain is increased due to sudden sarcomere stretch (49), we formulated *off*(*t*) such that elastance elements were lost from the elastance pool when volumetric distortion of these elements was increased due to sudden volumetric inflation. The effect of φ(*t*) on *off*(*t*) was given by (26) where γ is a coefficient that grades the dependence of *off*(*t*) on φ(*t*). To restrict this effect to stretch only, γ possessed value only when φ(*t*) > 0, and was set equal to zero when φ(*t*) < 0. Following a positive ΔV, this formulation had the effect to increase the rate at which elastance elements left the E(*t*) population. In and of itself, this effect would be transient, and post-ΔV kinetics would quickly reset to isometric levels as υ(*t*) returned to υ_{0}. To capture the asymmetry in the late-systolic ΔP(*t*) responses to positive and negative ΔV, the extra-displaced elements following a positive ΔV were removed from the *E*_{T} pool, i.e., they were taken out of the cyclic kinetic scheme in Fig. 3 such that they could not return to participate in pressure generation at a later time in the cycle. This was done by defining the extra-displaced elements, E_{L}(*t*) to be lost according to: Then, E_{T}(*t*) was redefined as E_{T}(t) = β[V(*t*) − V_{0}] − E_{L}(*t*). Thus, E_{L}(*t*) elements were lost from the potential recruitment pool.

#### Parameter summary.

There are 11 constant-valued model parameters, including four that carry units of s^{−1} (*k*_{on}, *k*_{off}, *k*_{a}, *g*); two that carry units of volume, μl (V_{0}, υ_{0}); two that carry units of elastance, mmHg/μl (β, E_{50}); and three that are dimensionless (*n*, ρ, γ). There is one time-varying parameter, α(*t*), that represents the net outcome of waxing and waning Ca^{2+} activation during the cycle, and it is also dimensionless.

## EXPERIMENTAL METHODS

Hearts for these studies were obtained from rats according to a protocol approved by the Washington State University Institutional Animal Care and Use Committee. All animals in this study received humane care in compliance with the animal use principles of the American Physiological Society and the Principles of Laboratory Animal Care formulated by the National Society of Medical Research and the National Institutes of Health's *Guidelines for the Care and Use of Laboratory Animals* (NIH Publication No. 85–23, Revised 1985).

Hearts were isolated from anesthetized (im 50 mg/kg ketamine, 5 mg/kg xylazine, 1 mg/kg acepromazine), young adult (2–6 mo), male rats (body wt 350–450 g) according to a protocol approved by the Washington State University Institutional Animal Care and Use Committee. Upon excision of the heart, the aorta was quickly cannulated, and the heart was immediately perfused with Krebs-Henselite solution (NaCl, 127 mM; KCl, 4.7 mM; NaH_{2}PO_{4}, 0.36 mM; CaCl_{2}, 1.25 mM; MgCl_{2}, 1.08 mM; NaHCO_{3}, 21 mM; glucose, 5 mM; pyruvate, 5 mM; insulin, 5 U/l; and bovine serum albumin, 0.08 g/l) containing high-dissolved O_{2} (Po_{2} > 600 mmHg) at a constant perfusion pressure of 100 mmHg.

Hearts were mounted onto an experimental system consisting of a constant-pressure perfusion subsystem, an environmental control subsystem, a pacing control subsystem, and a volume-servo subsystem. A latex balloon, attached to the obturator of the volume-servo subsystem and sized so as to not develop measurable pressure when inflated to 400 μl, was inserted into the LV chamber through the mitral orifice. The mitral annulus was secured to the obturator. Mounted and perfused hearts were then submerged in perfusate in a temperature-controlled environmental chamber that kept the epicardial surface wet and allowed for field stimulation. Electrical pacing was initiated by means of a computer-controlled stimulator and field electrodes placed in the bath of the environmental chamber on either side of the heart. The balloon was inflated to a baseline volume, V_{BL}, which was defined as the volume necessary to achieve an end-diastolic pressure of ∼5 mmHg. LV pressure, P(*t*), was measured with a 3-Fr catheter-tip Millar pressure transducer that was passed through a port in the volume-servo system, through the lumen of the obturator, and positioned in the balloon within the LV chamber.

The LV balloon was connected to a computer-controlled, volume-servo system with a metal-bellows piston in a servo chamber. Movement of the piston displaced volume out of, or into, the LV balloon according to the servo command. This allowed dynamic control and perturbation of LV volume with the simultaneous measurement of LV pressure.

Experiments were conducted on three groups of hearts; the temperature of each group was held at 35 (*n* = 9), 30 (*n* = 9), or 25°C (*n* = 10). Temperature was controlled by heating or cooling the perfusate and the water circulating through the jackets of the experimental subsystems. Pacing periods at each of the three temperatures were chosen just long enough to allow a well-defined diastolic period, wherein pressure persisted at the diastolic (passive) level without increasing or decreasing for ∼15% of the cycle period. According to this criterion, pace period was set at 500 ms in experiments conducted at 35 °C, at 800 ms in experiments conducted at 30 °C, and at 1,200 ms in experiments conducted at 25 °C.

Once stable isovolumic beating had been achieved at V_{BL}, ventricular volume was changed according to a five-volume, single-beat Frank-Starling protocol, which allowed evaluation of both systolic and diastolic pressure-volume relationships in each heart. In this protocol, specified volume changes, ΔV, were imposed 40 ms prior to the start of a test beat, and the test beat took place isovolumically at the altered volume. The smallest and largest volumes in the Frank-Starling protocol were 95% and 105% of V_{BL}, respectively, with intermediate volumes separated by 2.5% V_{BL} increments. After the Frank-Starling protocol, a volume-perturbation protocol was administered.

The volume perturbation protocol was designed to elicit a family of cycle-time-dependent LV pressure responses, ΔP(*t*), to incremental step-like volume changes, ΔV. This protocol was as follows. In an otherwise isovolumically beating heart, the volume was commanded at a specific cycle time to change by 2.5% of V_{BL} over a 20-ms period. In a sequence of ten records, the command for ΔV(*t*) was delivered at cycle times that advanced by 5% of the beat interval on each successive record. Thus in an 800-ms beat period, records were obtained with a ΔV(*t*) at each of 40-, 80-, 120-, 160-, 200-, 240-, 280-, 320-, 360-, and 400-ms cycle times into the beat period. In one subset of 10 records, ΔV*(t*) was equal to + 2.5% V_{BL}, and in a second subset of 10 records, ΔV(*t*) was equal to −2.5% V_{BL}. These two subsets of records [10 records each of ΔP(*t*) responses to both positive and negative ΔV(*t*)] were combined into one data set of 20 records for each heart.

At the end of the experiment, all atrial and great vessel tissue was trimmed from the heart and the heart was blotted dry and weighed. The right ventricular free wall was cut away from the ventricular septum, and the septum plus the LV free wall was weighed. A summary of average body and heart size for rats in each temperature group is given in Table 1.

## DATA HANDLING

### Processing Raw Data

For experiments conducted at 30°C, pressure records from 9 hearts (all with V_{BL} between 140 and 180 μl, all generating isovolumic pressures at V_{BL} between 150 and 200 mmHg, and all having approximately equivalent systolic and diastolic intervals) were analyzed. Because data to be analyzed from the incremental volume-perturbation protocol constituted small signals that were obtained by subtracting one large signal (pressure measured during an unperturbed beat) from another large signal (pressure measured during a perturbed beat), ensemble averaging techniques were used to minimize the influence of measurement noise and the influence of variation that arose because one beat in any string of beats was not an exact copy of another. Ensemble averaging over records obtained in all nine hearts during the Frank-Starling protocol was as follows. For the five isovolumic pressure records obtained at the five volumes in the Frank-Starling protocol, diastolic pressure in each record was determined by taking the average pressure during the last 0.1 s of the diastolic period. This diastolic pressure was assumed to be due entirely to passive LV properties, allowing the calculation of passive LV elastance, E_{p}, according to the slope of the regression of diastolic pressure on volume. Developed pressure in each beat was determined by subtracting diastolic pressure from the measured pressure. These developed pressure records, designated P_{iso}(*t*), were then ensemble averaged over the records from all nine hearts (Fig. 4*A*). An incremental effect on P_{iso}(*t*) by an incremental volume change prior to beat onset was calculated by subtracting the ensemble-averaged P_{iso}(*t*) at V_{BL} from ensemble-averaged P_{iso}(*t*) obtained at V_{BL} + ΔV, where ΔV equaled variously: 2.5%, 5.0%, −2.5%, and −5.0% of V_{BL}. This incremental effect was designated ΔP_{iso}(*t*) (Fig. 4*B*).

Ensemble averaging over records obtained in all nine hearts during the volume-perturbation protocol was as follows. In the records obtained from each of the nine hearts, the specific perturbed and unperturbed beats were identified from which the pressure response to a ΔV(*t*) at a specific time in the beat period was to be calculated. The pressure measurements taken from the perturbed and unperturbed beats were ensemble averaged over the records from all nine hearts. Such ensemble averaging was done for each of the 20 elicited ΔV responses.

Interest was in just the developed pressure, P(*t*), during a volume perturbation, which was determined by subtracting passive pressure [P_{p} = E_{p}(V_{BL} + ΔV)] from the measured pressure, P_{m}(*t*), i.e., P(*t*) = P_{m}(*t*) − P_{p}. The reference for all pressure responses to incremental volume changes was the isovolumic developed pressure, P_{iso}(*t*), at the baseline volume, V_{BL}. This P_{iso}(*t*) reference was taken from the pressure measurement in the isovolumic beat just prior to the beat receiving volume perturbation. The pressure response, ΔP(*t*) was then calculated as ΔP(*t*) = P(*t*) − P_{iso}(*t*). An example of such a calculation of ΔP(*t*) to a ΔV delivered near the time of peak P_{iso}(*t*) is given in Fig. 5.

For experiments conducted at 25 and 35°C, ensemble average records for each temperature were obtained from data collected in 10 and 9 hearts, respectively. Ensemble averages were taken of data collected from hearts only if these hearts demonstrated uniformity in V_{BL} and P_{iso}(*t*), which was the case for all hearts that generated the data reported here.

### Fitting Model to Data and Parameter Estimation

It was necessary to estimate model parameters prior to model evaluation in terms of its ability to reproduce the cycle-time-dependent ΔP(*t*) responses to ΔV. The Ca^{2+} activation coefficient, α(*t*) was estimated from given values of the constant-value parameters, and an isovolumic pressure measurement, P_{iso}(*t*), was obtained from an isovolumic beat at a baseline volume, V_{BL}. In an isovolumic beat, the volumetric-distortion variable, υ(*t*), never deviates from υ_{0}. Thus, from Eq. 11, (27) With υ_{0} a constant, P_{iso}(*t*) serves as a surrogate for E(*t*). After substitution of the P_{iso}(*t*) surrogate for E(*t*) into the first-order differential equation for E(*t*) (Eq. 20), the equation can be rearranged to solve for α(*t*) in terms of known quantities as (28)

The parameter *n* in the cooperative activation equation (Eq. 23) was set equal to 2. The remaining 10 constant-valued parameters (*k*_{on}, *k*_{off}, *k*_{a}, β, E_{50}, V_{0}, υ_{0}, ρ, γ, *g*) were estimated by conducting nonlinear regression procedures as follows. Best guess estimates of the 10 free parameters were chosen. The first time derivative of the measured P_{iso}(*t*) at V_{BL} was calculated. Then, α(*t*) was calculated according to Eq. 28. A volume-perturbation paradigm was defined computationally to simulate the volume perturbation protocol used in the experiments. The calculated α(*t*) and computationally defined ΔV(*t*) were fed as inputs into the differential equations (Eqs. 20, 21), which were then numerically integrated (4th order Runge-Kutta) to yield predictions of the state variables: E(*t*) and υ(*t*). These state variables were then multiplied according to Eq. 11 to yield a prediction of P(*t*), from which the ΔP(*t*) responses to ΔV(*t*) were calculated. The 10 free parameters were then submitted to a formal optimization procedure (lsqnonlin, MATLAB) that minimized the sum of square errors between predicted and observed ΔP(*t*). To avoid local minima, the optimization procedure was started with widely different initial values of parameters. Repeated convergence to the same set of final values was taken as evidence that a global minimum had been found. Successful termination of the optimization procedure yielded the set of parameter values that best matched the model with the data.

### Model Reproduction of ΔP(t) Responses

As an outcome of parameter estimation, the best possible model reproduction of the ΔP(*t*) responses was obtained. Model reproductions were evaluated according to the quality by which the various phases of the ΔP(*t*) responses were replicated and the manner in which these phases changed over the course of the cycle and with time of ΔV perturbation. Because qualitative evaluation was more important than quantitative evaluation, the only quantitative evaluative parameters that were obtained were the slope and *r*^{2} of the linear regression of model-generated ΔP(*t*) responses upon the measured ΔP(*t*) responses.

## RESULTS

### Experimental Results

Descriptive features of pressure responses to incremental volume change will be given for data collected at the middle temperature, 30°C. These features will then be compared with equivalent features in data obtained from hearts at lower and higher temperatures: 25 and 35°C, respectively.

#### Pressure response to ΔV administered prior to beat onset.

The five P_{iso}(*t*) values recorded during the single-beat Frank-Starling protocol (Fig. 4) were used to evaluate the uncomplicated effects of volume on pressure around V_{BL}. Because volume change occurred prior to beat onset in these isovolumic beats, there was no effect of rate-of-change of volume (i.e., flow) on pressure development, as occurs when ΔV is administered during the course of the beat. Further, because volume changes occurred only 40 ms before beat onset, it is unlikely that the volume change caused a change in activating calcium transients, for the reason that such calcium-transient changes take place, if they take place at all, only after several beats following a volume change (35). Thus differences in P_{iso}(*t*) at the various volumes reflect just the effects of differences in volume on pressure development.

The effect of incremental volume changes (Fig. 4) was: *1*) to systematically change the amplitude of P_{iso}(*t*) in the direction of the volume change; *2*) to not detectably change the time of peak P_{iso}(*t*) (210 ms for all 5 isovolumic beats); *3*) to not appreciably change the shape of the ascending limb of P_{iso}(*t*) [the ascending limbs of the five P_{iso}(*t*) values are scaled versions of one another]; and *4*) to systematically slow the fall of P_{iso}(*t*) on its descending limb, as can be seen in the divergence of the descending limbs of the five P_{iso}(*t*) value*s* and the fact that the time to fall from peak P_{iso}(*t*) to half peak value was longest (200 ms) for P_{iso}(*t*) at the largest volume and shortest (189 ms) for P_{iso}(*t*) at the smallest volume.

Additional insight into ΔV effects on P_{iso}(*t*) is obtained from examining the slope of the ΔP_{iso}/ΔV isochrones at different times during the cycle, as determined by linear regression of P_{iso} vs. V at times equal to 115, 210, and 410 ms, respectively (Fig. 4*C*). All regressions yielded *r*^{2} > 0.99, indicating apparently linear P_{iso} vs. V isochronal relationships over this small volume range (±5% V_{BL}) at each of these times. Note that over a larger volume range, P_{iso} vs. V isochrones can be nonlinear, with the nonlinearity becoming very pronounced during late systole and relaxation (8, 19, 33, 36). The slope of the P_{iso} vs. V isochrones (ΔP_{iso}/ΔV = E_{τ}) increased progressively with pressure as pressure rose on the ascending limb; e.g., at 115 ms, P_{iso}(*t*) at V_{BL} was 98 mmHg and E_{τ} was 0.633 mmHg/μl, whereas at 210 ms [time of peak P_{iso}(*t*)], P_{iso}(*t*) at V_{BL} was 173.4 mmHg and E_{τ} was 1.019 mmHg/μl. However, E_{τ} did not change progressively with pressure as pressure fell on the descending limb; e.g., at 410 ms, P_{iso}(*t*) at V_{BL} was 97 mmHg (approximately the same as at 115 ms), but E_{τ} was 0.990 mmHg/μl, approximately the same as it was 200 ms earlier at the time of peak P_{iso}(*t*).

These time variations in P_{iso} vs. V isochrones translate into characteristic features of the ΔP_{iso}(*t*) wave shape. The linear isochrones meant that ΔP_{iso}(*t*) was essentially symmetric in its response to positive and negative ΔV over a range equal to ±5% of V_{BL}. This apparent symmetry in ΔP_{iso}(*t*) responses to positive and negative ΔV is seen in the similar shapes of waveforms *a* and *c*, and again in waveforms *b* and *d* (Fig. 4*B*). Thus, marked asymmetry in the ΔP(*t*) responses to positive and negative ΔV delivered during the course of the cycle will have to be explained by something other than ΔV effects by themselves. Second, the manner in which the time course of E_{τ} differed from the time course of P_{iso}(*t*) meant that the ΔP_{iso}(*t*) waveform differed from that of P_{iso}(*t*). This is seen when the amplitude-normalized waveforms of P_{iso}(*t*) and ΔP_{iso}(*t*) are compared (Fig. 6). With ΔP_{iso}(*t*) to a ΔV of +2.5% V_{BL} as the reference (waveform *a* in Fig. 4), the amplitude-normalized ΔP_{iso}(*t*) waveform follows (or, at most, slightly leads) the amplitude-normalized P_{iso}(*t*) waveform on the ascending limb, reflecting the common trajectories of the ascending limbs of P_{iso}(*t*)s in Fig. 4. However, once a peak value was reached, the ΔP_{iso}(*t*) wave shape persisted with near-maximum amplitude for ∼200 ms, while P_{iso}(*t*) was falling to approximately half its peak value, reflecting the divergence of trajectories of the descending limbs of P_{iso}(*t*) in Fig. 4. Thereafter, ΔP_{iso}(*t*) quickly falls to reach zero at about the same time P_{iso}(*t*) reaches zero. The divergence of the ΔP_{iso}(*t*) trajectory from P_{iso}(*t*) is another way of viewing the previously reported phenomenon whereby the characteristic time of relaxation increases following increased pressure development in the heart (54) just as it increases with increased force development in isolated cardiac muscle (31).

The time courses of both P_{iso}(*t*) and ΔP_{iso}(*t*) were used as reference trajectories to describe the time-courses of pressure responses to ΔV administered at different cycle times.

#### Pressure responses to ΔV administered at various cycle times.

An example of records from which a calculation was obtained of ΔP(*t*) response to a ΔV delivered near the time of peak P_{iso}(*t*) is given in Fig. 5. The complete family of 20 ΔP(*t*) values in response to ΔV administered at various times in the cycle is given in Fig. 7. In brief, the descriptive characteristics of ΔP(*t*) responses in beats in which ΔV was delivered at all times of the cycle are as follows:

The ΔP(*t*) response [different from the ΔP_{iso}(*t*) response] consisted of three distinct phases: *phase 1*, coincident with the leading edge of ΔV; *phase 2*, a rapid recovery from the *phase 1* excursion back towards zero; and *phase 3*, sustained aspect of the response lasting until the end of relaxation (Figs. 5, 7–9).

The peak of *phase 1* amplitude more or less tracked a P_{iso}(*t*) envelope, with the exception that peak *phase 1* amplitude of very early ΔP(*t*) responses reached above the envelope (Fig. 7).

*Phase 2* became progressively more pronounced and extended for longer time periods as ΔV was delivered later in the cycle (Figs. 8 and 9).

*Phase 3* was associated with the ΔP_{iso}(*t*) trajectory (Fig. 8).

*Phase 3* of early ΔP(*t*) responses began quickly with initial amplitudes above the ΔP_{iso}(*t*) trajectory, but eventually joined with that trajectory by midsystole and completed the response coincident with ΔP_{iso}(*t*). These early-cycle *phase 3* responses were said to be exaggerated (Fig. 8).

*Phase 3* of late ΔP(*t*) responses began slowly and rose to amplitudes below the ΔP_{iso}(*t*) trajectory. *Phase 3* of late ΔP(*t*) responses to positive ΔV never joined with the ΔP_{iso}(*t*) trajectory before the response terminated. Late-cycle *phase 3* responses were said to be diminished (Fig. 8).

ΔP(*t*) responses to positive and negative ΔV were asymmetric about the time axis in both *phase 2* and *phase 3*, with the asymmetry in the *phase 3* response becoming more pronounced for late-cycle ΔP(*t*) (Fig. 9). This contrasted with the symmetry in ΔP_{iso}(*t*) obtained from beats with positive and with negative ΔV (Fig. 5).

#### Temperature dependence of ΔP(t) responses to step-like ΔV.

The most pronounced difference among responses collected at the three different temperatures was the speed at which pressure rose and fell during isovolumic beating. To meet our requirement for approximately equal percentage of time in diastole and systole during the isovolumic beat at each temperature, the heart was paced with a beat period of 1,200 ms at 25°C, 800 ms at 30°C, and 500 ms at 35°C. To indicate the relative time course of pressure development during beats at each temperature, the ensemble average P_{iso}(*t*) at each temperature is given in the right column of Fig. 10. Despite the different pacing periods, peak P_{iso}(*t*) at V_{BL} was, on the average, the same in hearts at all three temperatures, i.e., ∼180 mmHg.

All features of the cycle-time-dependent ΔP(*t*) responses to step-like ΔV described for the responses obtained at 30°C were also present in responses obtained at 25 and 35°C (Fig. 10). Further, these ΔP(*t*) response features changed with cycle time in the same characteristic pattern at 25 and 35°C as they did at 30°C. This consistency in the ΔP(*t*) response patterns over a 10°C range of temperatures suggests that there was one mechanodynamic pattern for ΔP(*t*) responses. This means that a model that reproduced this pattern at one temperature would then apply, with proper parameter adjustment, to the LV over the entire range of temperatures.

### Model Results

There were two types of model-generated results: *1*) model fits to the ΔP(*t*) response data in which model parameters were estimated, and *2*) model predictions of isovolumic pressures, which represented prediction of measured data other than that to which the model was fit.

#### Model fits to ΔP(t) response data.

Following parameter estimation, the model produced the 20 ΔP(*t*) responses (Fig. 7, *right*) that are to be compared with the measured ΔP(*t*) responses (Fig. 7, *left*). The goodness of fit may be appreciated from a plot of model-generated vs. measured ΔP(*t*) responses (Fig. 11). Linear regression of this data resulted in a standard error of the estimate equal to 0.42 mmHg and *r*^{2} equal to 0.97. Thus errors were small, with the consequence that only 3% of the variance in the data was not accounted for by the model. Although equal to 0.97, the slope of the regression line was significantly different from 1 (*P* < 0.01). This meant that there was some, albeit small, bias in the ΔP(*t*) and parameter estimates. Parameter estimates from the best fit to the 30°C ΔP(*t*) response family are given in Table 2.

The qualitative correspondence between the detailed features of model-generated and measured ΔP(*t*) responses were just as important as the quantitative aspects. This correspondence is given in Figs. 7, 8, and 9, where the overlay of the model-generated on the measured ΔP(*t*) responses in Fig. 9 is so close that discriminating between them is difficult at the scale of the figure. All model-generated ΔP(*t*) responses exhibited the three-phase response characteristic of the measured ΔP(*t*), and these phases changed with time in the cycle in much the same manner as did the measured ΔP(*t*). For instance, the amplitude of *phase 1* of the model-generated ΔP(*t*) tracked the P_{iso}(*t*) envelope with some tendency [just was with the measured ΔP(*t*)] for early ΔP(*t*) to reach above this envelope (Fig. 7). Further similar to the measured data, *phase 2* was small and indistinct for early responses but became more prominent during later responses.

Model-generated *phase 1* was most sensitive to the isovolumic-distortion parameter (υ_{0}), and the rate of recovery of model-generated pressure during *phase 2* was most sensitive to the distortion rate constant (*g*). Thus *phases 1* and *2* were due to the distortion component of the model. Without the distortion equation in the model, the model did not predict a response with *phases 1* and *2*. The distortion response was initiated at the leading edge of ΔV and quickly died away.

Similar to *phase 3* of measured ΔP(*t*) responses, *phase 3* of model-generated ΔP(*t*) responses was closely associated with ΔP_{iso}(*t*) (Fig. 8). Further [like measured ΔP(*t*)], *phase 3* of early-cycle, model-generated ΔP(*t*) responses exhibited the exaggerated early-cycle response and began well above the ΔP_{iso}(*t*) trajectory, but eventually joined with that trajectory by midsystole and completed the response coincident with ΔP_{iso}(*t*) (note response *2* in Fig. 8). Also [like measured ΔP(*t*)], *phase 3* of late-cycle, model-generated ΔP(*t*) responses exhibited a diminished late-cycle response. *Phase 3* of late-cycle, model-generated ΔP(*t*) responses to positive ΔV began well below the ΔP_{iso}(*t*) trajectory (note arrows in responses *6* and *8* in Fig. 8) and never joined with the ΔP_{iso}(*t*) trajectory before the response terminated. However, *phase 3* of late model-generated ΔP(*t*) responses to negative ΔV began only slightly below the ΔP_{iso}(*t*) trajectory but approached and joined with the ΔP_{iso}(*t*) trajectory before the response terminated.

Additional similarities between model-generated and measured ΔP(*t*) responses were observed in the character of the asymmetry about the time axis in the responses to positive and negative ΔV (Fig. 9). Both model-generated and measured ΔP(*t*) showed little asymmetry in early-cycle responses but exhibited progressively greater asymmetry in late-cycle responses. A notable aspect of the response asymmetry in both measured and model-generated responses was that *phase 3* in the response to negative ΔV was larger in magnitude that that in response to positive ΔV (Fig. 9).

*Phase 3* was due entirely to the recruitment response. Its amplitude was sensitive to β and V_{0}; its rate of rise and fall were sensitive to *k*_{on} and *k*_{off} as well as E_{50}; and the extended time course of *phase 3* was sensitive to the cooperative activation parameters of *k*_{xb}(*t*), *k*_{a}, and E_{50}. The stretch activation parameter, ρ, was responsible for the exaggerated early-cycle *phase 3* response (when ρ was set to zero, there was no exaggerated early-cycle response). The asymmetry in late-cycle *phase 3* responses to positive vs. negative ΔV was due to the parameter γ (when γ was set to zero, late-cycle *phase 3* responses were symmetric about the ordinate).

There was one qualitative feature of the measured response that the model did not reproduce with fidelity. Note in Figs. 8 and 9 the pronounced nadir in the measured responses in the transition between *phases 2* and *3* followed by a rebound in *phase 3* during late ΔP(*t*) responses (especially in responses to negative ΔV). In contrast, the model generated only a shallow nadir and a weak rebound response. In terms of overall goodness of fit, this qualitative model discrepancy contributed to the <3% of the variance in the data that was not accounted for by the model. Whether this model discrepancy is important to an application will depend on the goals of the application.

#### Time course of model recruitment parameters.

Model construction was based on kinetic processes whereby elastance elements were recruited into the elastance state according to a recruitment kinetic coefficient, *on*(*t*). There were two mechanisms for recruitment: calcium activation as embodied in the coefficient *k*_{on}α(*t*), and cooperative cross bridge activation as embodied in the coefficient *k*_{xb}(*t*). The relative time course and magnitude of these mechanisms are important to interpreting model results and to the model's ultimate utility. As expected from a calcium-activating process, *k*_{on}α(*t*) rose early before the rise of P_{iso}(*t*), peaked well before the time of peak P_{iso}(*t*), and fell to zero while P_{iso}(*t*) continued (Fig. 12). And, as expected from a cross bridge-based activation process, *k*_{xb}(*t*) rose some time after *k*_{on}α(*t*) and persisted to become the primary mechanism for recruiting elastance elements and sustaining P_{iso}(*t*) during the late phases of the cycle (Fig. 12). What was not expected was that the peak amplitude of *k*_{xb}(*t*) was higher than that of *k*_{on}α(*t*). The implication was that cooperative activation, as a mechanism for recruitment, is at least as strong as is calcium activation. Be that as it may, recruitment is initiated by calcium activation, and no recruitment from cooperative mechanisms could occur without preceding recruitment by the Ca^{2+}-based process, *k*_{on}α(*t*).

#### Model predictions of isovolumic pressure.

Isovolumic pressures were not used in the data set to which the model was fit. Thus after parameters were estimated from fitting to the ΔP(*t*) response family, model reproduction of P_{iso}(*t*) and the resulting calculations of ΔP_{iso}(*t*) were true predictions. Such predictions are shown in Fig. 4, where it is demonstrated that the model generated apparently linear pressure-volume isochrones over ±5% volume variations, and the temporal variation in these model-predicted isochrones was much the same as that observed in the measured data (Fig. 4*C*).

Model-predicted ΔP_{iso}(*t*) wave shapes were images with respect to amplitude and time-course of the measured ΔP_{iso}(*t*) wave shapes (compare curves *a*–*d* in the right and left columns of Fig. 4*B*). Thus model-predicted ΔP_{iso}(*t*) extended with high amplitude well beyond the time of peak P_{iso}(*t*) in a manner similar to that observed in the measured ΔP_{iso}(*t*) (Fig. 6). The extended time course of model-predicted ΔP_{iso}(*t*) was due entirely to the cooperative recruitment parameter, *k*_{xb}(*t*), and the amplitude of this extended time course was very sensitive to the *k*_{xb}(*t*) scalar, *k*_{a}. In addition, amplitude of the model-predicted ΔP_{iso}(*t*) wave shape was sensitive to β and V_{0}, and the rate of rise and fall were sensitive to *k*_{on}, *k*_{off}, and E_{50}. These were the same parameter sensitivities as found for *phase 3* of the ΔP(*t*) response.

#### Temperature dependence of model parameters.

The model fit the families of ΔP(*t*) responses obtained at 25 and 35°C as well as it fit the ΔP(*t*) responses obtained at 30°C; *r*^{2} equaled 0.96 for fits at both 25 and 35°C. The estimated parameter values differed between fits at 25 and 35°C, and these differed from parameters estimated at 30°C (Table 2). As a test of model validity that went beyond how well the model fit the data, the estimated parameters were scrutinized as to whether the difference in estimated values at the different temperatures made sense. Among the four parameters that carried units of s^{−1} (*k*_{off}, *k*_{on}, *g*, and *k*_{a}), the expected systematic increase in value with increase in temperature was seen only in the distortion rate constant, *g*, which increased by a factor of 1.83 over the 10°C range. However, among these four parameters, only *g* was not multiplied by a time-varying factor, and therefore, traditional interpretations of temperature dependence of constant-valued kinetic coefficients may apply only to *g*. Thus the temperature changes in recruitment kinetics were all contained within the time course of the Ca^{2+} activation parameter α(*t*) and were not reflected in the scalar coefficients *k*_{off}, *k*_{on}, and *k*_{a}.

Interestingly, υ_{0} showed very little difference among the three groups of hearts. When expressed as a fraction of V_{BL}, υ_{0} was 0.20, 0.23, and 0.18 in hearts at 25, 30, and 35°C, respectively. Therefore, a sudden removal of 20% of V_{BL} would completely depressurize these hearts regardless of temperature. Because only elastic physical properties are involved in this depressurization, it is consistent with expectations that this parameter would not change with temperature.

One bothersome outcome of these temperature-dependent parameter estimates is the calculated differences in E_{T} [E_{T} = β(V_{BL} − V_{0})] among hearts at these three temperatures. From the theory on which the model was based, E_{T} depends on structure and could not change with temperature. However, calculated E_{T} at V_{BL} ranged from 46.4 mmHg/μl for hearts at 30°C to 7.6 and 7.3 mmHg/μl for hearts at 25 and 35°C, respectively. Considering that E_{T} represents the number of elastance elements that could potentially be recruited into pressure generating states, the values estimated at 25 and 35°C (potentially producing pressures of 263 and 282 mmHg, respectively) are much more believable than the value estimated at 30°C (potentially producing a pressure of 1,653 mmHg). We are unable to comment on this seeming model inadequacy, whereby unrealistic potential pressure-generating capability is predicted in some circumstances.

## DISCUSSION

We built a model of LV pressure-volume dynamics starting with myofilament contraction processes. This derivation led naturally to treating separately the dynamics of processes responsible for the recruitment of pressure-generating elements and the dynamics of processes responsible for elastic deformation (distortion) of these elements. We dubbed this the Recruitment-Distortion-Instantaneous-Elastance (RDIE) model. Experimental observations of cycle-time-dependent ΔP(*t*) to incremental step-like ΔV were then used to inform the parametric character of model kinetic coefficients. Important outcomes of these experiments are as follows.

### Pressure Responses to Volume Changes Delivered Prior to Beat Onset Suggested an Important Role for Cooperative Activation

The ΔP_{iso}(*t*) wave shape was markedly different from the P_{iso}(*t*) wave shape in that it extended with high amplitude to well beyond the time of peak P_{iso}(*t*) (Figs. 4, 6). To account for this, we note that both myocardial activation and LV volume determine LV pressure generation. The difference in wave shape between ΔP_{iso}(*t*) and P_{iso}(*t*) suggests that there is interaction between these two causative factors (activation and volume). We assumed that Ca^{2+} activation was the same in all isovolumic beats whether at V_{BL}, V_{BL} + ΔV, or V_{BL} − ΔV. If the effect of volume in determining pressure was linearly independent of the effect of activation, then P_{iso}(*t*) values at different volumes would simply be scaled versions of one another and the ΔP_{iso}(*t*) wave shape would be identical to that of P_{iso}(*t*) at V_{BL}. The marked difference in ΔP_{iso}(*t*) and P_{iso}(*t*) wave shapes means that there is some nonlinear interaction in the determination of P_{iso}(*t*). Furthermore, in as much as this prolonged high-amplitude wave form is an essential aspect of the *phase 3* response to ΔV delivered early in systole, but is a less pronounced part of responses to ΔV delivered late in systole, the nonlinear interaction is also cycle-time-dependent. A significant challenge was to represent this nonlinear interaction and its cycle time dependence.

We responded to this challenge by incorporating a cooperative activation term, *k*_{xb}(*t*), into the model. The formulation of *k*_{xb}(*t*) was such (Eq. 23) that the transition of elastance elements into the elastance state favored the further transition of more elements into that state, but this could not be sustained indefinitely in the absence of an initiating Ca^{2+} activation episode. However, the incorporation of *k*_{xb}(*t*) into the model enabled contraction to be prolonged for periods beyond what could be supported by Ca^{2+} activation alone. Once Ca^{2+}-based activation initiated the formation of pressure-bearing elastance elements, cooperative activation became ever more important until it was the sole activation mechanism during late systole (Fig. 12). Because these cooperative activation effects were most pronounced during late systole, they had a profound effect on the onset and rate of relaxation, causing beats that generated higher pressure to relax slower. This was responsible for the sustained high amplitude of ΔP_{iso}(*t*). Our findings confirm the claims of others (3, 34, 42, 49) that cooperative activation is important for sustaining contraction and the systolic event.

### ΔP(t) Responses to Step-Like ΔV Are Reminiscent of Force Response to Quick Stretch and Release of Constantly-Activated Muscle

Even though these ΔP(*t*) responses to step-like ΔV were elicited from a contractile system that was undergoing temporal variations in activation, and even though these responses were different at different times in the heart cycle, each ΔP(*t*) response exhibited features that were reminiscent of the force response to quick stretch and release in constantly activated cardiac muscle. Similar features include a component of the response that was coincident with the forcing perturbation (*phase 1*); a component of the response in which there was relatively rapid recovery from the initial *phase 1* disturbance (*phase 2*); and a component of the response that was sustained and, at later times in cycle, slowly developing (*phase 3*). These similarities beg an interpretation of the ΔP(*t*) response in terms of the accepted interpretations for the equivalent phases in the response of constantly activated cardiac muscle undergoing quick stretch (see Ref. 49 for a brief synopsis of the quick stretch response of cardiac muscle). *Phase 1* in the quick-stretch muscle response represents force as a result of stretching elastic, force-bearing (attached) myosin cross bridges (XB) beyond the stretch imposed in the formation of these XB in the isometric state. *Phase 2* represents the detachment of these overstretched XB and their replacement with cycling XB possessing only isometric strains. *Phase 3* represents length-induced recruitment of additional force-bearing XB. This reproduction in the LV pressure response to quick volume change of muscle force response to quick stretch and release is strong justification for the method of construction of the RDIE model and using contractile processes as the starting point.

By analogy with elastic XB in cardiac muscle, the RDIE model incorporated volumetric elastance elements. These elements undergo baseline volumetric distortion, υ_{0}, when formed, resulting in P_{iso}(*t*). The leading edge of a step-like ΔV causes additional volumetric distortion of these elements, resulting in *phase 1* of the ΔP(*t*) response. That *phase 1* amplitude tracked P_{iso}(*t*) (Fig. 7) is evidence for the instantaneous elastance construct. Thus the instantaneous LV elastance is proportional to P_{iso}(*t*) at all times of the cycle, just as muscle stiffness is proportional to muscle force.

By analogy with cycling XB in cardiac muscle, cycling elastance elements in our RDIE model recovered from the elastic distortion imposed during *phase 1* to produce *phase 2*, which was a consequence of the detachment of distorted elastance elements and their replacement with newly-formed elastance elements that had not experienced the distorting influence of ΔV. Cycling elastance elements are also consistent with our earlier demonstration that, for ΔV imposed at the time of peak P_{iso}(*t*) in isolated rabbit and ferret hearts, changes in *phases 1* and *2* with variations in ΔV amplitude and ΔV rate are consistent with the elastic distortion hypothesis (11, 13, 16, 17, 43).

Finally, by analogy with *phase 3* of the quick stretch response of cardiac muscle, *phase 3* of the ΔP(*t*) response in the RDIE model represented the ΔV-induced recruitment of either more (for +ΔV) or fewer (for −ΔV) pressure-generating elastance elements. We base this on the relationship of *phase 3* to ΔP_{iso}(*t*). Because *phase 3* of the ΔP(*t*) response follows ΔP_{iso}(*t*), there can be no other conclusion than that the LV, upon receiving a volume increase during systole, enhances its pressure-generating capacity just as it would have done had the beat begun at the greater volume. Thus ΔV during systole recruits pressure-generating capacity.

### Instantaneous LV Elastance Differs from Traditional Isochronal Elastance

Despite the conventional appearance of Eq. 11, both elastance, E(*t*), and volumetric distortion, υ(*t*), in that equation differ substantially from the traditional definitions of LV elastance and of the volume associated with that elastance to determine pressure. For instance, in the conventional time-varying elastance formulation (51), LV elastance (ΔP/ΔV) is calculated as the slope of a linear pressure vs. volume relationship at a given cycle time (isochrones), where pressure and volume variations are obtained from beats taking place at different volumes. Indeed, we followed this procedure to obtain the P_{iso} vs. V isochronal relationships in Fig. 4. Based on isochronal pressure-volume relationships and the isochronal elastance calculated from these isochrones, a large fraction of the chamber volume is apparently involved in the elastic deformation of wall material. As an example, the volume associated with isochronal elastance to predict pressure at the time of peak isovolumic pressure in Fig. 4 is ∼100% of the total LV volume (near zero volume axis intercept). However, if LV volume is removed suddenly at any time during the heart beat, the volume of elastic distortion is found to be only 18–20% of the chamber volume (from Tables 1 and 2) and less than this in other studies (5, 11, 13, 16, 43). Correspondingly, chamber elastance calculated according to the conventional isochronal formulation is several times smaller than elastance calculated from sudden volume withdrawals (e.g., isochronal elastance at peak systole calculated from the data in Fig. 4 is ∼1 mmHg/μl, whereas instantaneous elastance at peak systole calculated from the ΔP response in Fig. 5 is at least 5 mmHg/μl). Thus the E(*t*) in Eq. (1) needs to be designated the “instantaneous chamber elastance” to distinguish it from the “isochronal chamber elastance”.

A second difference stems from the fact that the conventional isochronal elastance waxes and wanes over the cycle as if it is a result of variations in myocardial activation/contraction, and it is not influenced by chamber volume. In contrast, the instantaneous elastance in Eq. 11 is driven to change both by activation and by chamber volume according to Eq. 21. A third difference is that, whereas the conventional isochronal elastance varies with time but does not have dynamic properties, instantaneous elastance in our formulation is a dynamic variable that evolves over time in response both to input functions (Ca^{2+} activation and volume) and to its own history.^{1} A fourth difference is that volumetric distortion, υ(*t*), is also a dynamic variable that changes according to its own history as well as in response to the rate of change of volume, but not according to volume itself. Repeating this important feature of the current model, volumetric distortion of LV elastance is not related to chamber volume per se, but to the rate of change of volume. These points of difference arise because we used the dynamic relationship between force and length of a cardiac myofiber as the centerpiece for the dynamic relationship between LV pressure and volume.

### Model Uniqueness

There are many models for representing LV function, so many that only a small fraction dealing with global LV pressure-volume relationships will be mentioned here. As described in the introduction, several LV models have been based on the time-varying isochronal elastance concept and its modification with resistance and deactivation elements (2, 5, 12, 15, 21, 24, 29, 44, 45, 52, 53, 55, 58). Others, recognizing the need to build a bridge between global LV properties and underlying contractile mechanisms, have complemented isochronal elastance with muscle contractile phenomena to improve on the model representation (4, 37). Still others have built their models from the ground up by starting with the kinetics of muscle contractile processes and extending them to their ultimate mechanical expression at the level of the LV (3, 9, 11, 32, 42, 50). While relying on contractile kinetics as the core, most of these models have also incorporated the muscle force-velocity relationship as a necessary adjunct to account for pressure behavior during LV ejection.

Choosing between models depends on the user's objectives and confidence in the model. For instance, those who want a simple model to convey basic notions of LV pressure-volume behavior, especially for didactic purposes, would be well served to use the classical time-varying, isochronal elastance model. Those who measure the LV global Ca^{2+} activation signal and, simultaneously, the LV pressure will find the model of Burkhoff and colleagues (3, 9, 42, 50) valuable for establishing the status of the many factors that convert the Ca^{2+} input signal into LV pressure. Further, just as the classical time-varying isochronal elastance model was useful to us as a beginning point in thinking about LV mechanics, the Burkhoff model (9) was valuable to us in that it formatted a kinetic pattern whereby force generation from XB could come from both Ca^{2+}-bound and non-Ca^{2+}-bound myofilament states. In summary, some of the previous models continue to be useful for analyzing and interpreting experimental data, and some are useful as background for constructing new models.

There are at least five features of the RDIE model that make it unique. One, the model is couched in terms of newly defined and identifiable global LV physical properties, i.e., instantaneous elastance and the volumetric distortion of that elastance. Two, the driving inputs for changes in instantaneous elastance are both activator Ca^{2+} and circumferential fiber length secondary to volume change. Third, the driving input for volumetric distortion change is rate of change of circumferential fiber length secondary to rate of change of volume. The distortion component of the RDIE model is unlike any part of previously proposed models and, in addition to being responsible for *phases 1* and *2* of the ΔP(*t*) response, acts to reduce pressure during ejection. Thus it takes the place of muscle force-velocity relationship in other models. Fourth, variable-dependent kinetic coefficients required to represent the global LV expressions of cooperative XB activation, stretch activation (shortening deactivation), and strain-dependent XB detachment take their form from global behaviors such as extended time-course of ΔP_{iso}(*t*), exaggerated *phase 3* response to early-cycle ΔV, diminished *phase 3* response to late-cycle ΔV, and asymmetry in ΔP(*t*) response to positive and negative ΔV. Five, the RDIE model is the only model that has been shown to reproduce the cycle-time-dependent ΔP(*t*) responses to step-like ΔV. In these ways the RDIE model is unique among the many previously proposed LV models.

Earlier, we described a different version of the recruitment-distortion LV model (18) that, in fact, does a credible job of reproducing the cycle-time-dependent ΔP(*t*) response to step-like ΔV. This earlier model can be considered a phenomenologic version of the current model. This is best seen by referring to the input term on the far right of Eq. 20, which is a product of an activation-dependent term, *on*(*t*), and a volume-dependent term, β[V(*t*) − V_{0}]. In the terminology of the previous approach, the effect on pressure by the *on*(*t*) input alone is through a property that was called “activance”, whereas the effect on pressure by V(*t*) alone is through the property that was called “elastance”. However, because *on*(*t*) and V(*t*) combine as a product and thus are not independent in their effects, there is an additional effect on pressure due to their interaction, and this is what was called “interactance” in the previous manuscript. The nonlinearities in the current model (cooperativity, stretch activation, and distortion-dependent loss) complicate this interpretation somewhat but not excessively. The previous approach worked well for analyzing LV responses to sinusoidal volume perturbations given throughout the cycle and remains useful for coupling experimental findings of myocardial stiffness in constantly-activated cardiac muscle with LV pressure-volume behavior. It does not work as well for the cycle-time-dependent step-response data of the current manuscript, where the current model is to be preferred. The current model also has the advantage of being tangibly linked, rather than phenomenologically linked, to underlying contractile mechanisms.

### Model Critique

All models are approximations of reality. Their value comes from the uses to which they may be put. These uses may be categorized as description, prediction, and explanation. Descriptive uses of models are for characterization of systems that give rise to observed behavior. The first criterion for characterizing an observed system is the reproduction, with acceptable fidelity, of those features of behavior that are important enough to be noted. By our judgment, the current RDIE model adequately reproduced the essential features of the cycle-time-dependent ΔP(*t*) responses to ΔV, including the three distinct ΔP(*t*) phases that changed in amplitude and in relationship to one another with cycle time, and also the extended time course of ΔP_{iso}(*t*). Therefore, an initial conclusion is that the model is descriptively valid and may be used to characterize the system (i.e., the LV) that gave rise to the observed behaviors.

However, the second criterion for descriptive validity is that model parameters be uniquely estimated. There were 10 fixed-value parameters estimated; eight of these were associated with the recruitment component (*k*_{on}, *k*_{off}, *k*_{a}, β, E_{50}, V_{0}, ρ, γ) and two were associated with the distortion component (*g* and υ_{0}). We encountered difficulty in the optimization algorithm satisfying convergence criteria even after many thousands of iterations. This meant that the square-error surface in parameter space on which the algorithm was attempting to find the least square error did not have a well-defined minimum; an indication that the parameters were not independent in their effects. It appeared that where parameters combined in product terms, such as β and V_{0} in the expression for E_{T} and *k*_{a} and E_{50} in the expression for *k*_{xb}(*t*), they were co-dependent. Further, because length-dependent recruitment and Ca^{2+}-dependent recruitment [arising from the time-varying parameter α(*t*)] appeared in the equations in product relations, clear separation of each of these recruitment effects could not be made based on the ΔP(*t*) response family alone. Thus parameters contained in the recruitment component of the model were not estimated uniquely. (Note the earlier concerns regarding the E_{T} calculation in hearts at different temperatures using the estimates of β and V_{0}.) Either reformulation of these expressions or acquisition of extra data containing information beyond what is contained in the ΔP(*t*) response family will be required to allow unique estimation of recruitment parameters.

In contrast to the difficulty encountered with unique estimation of recruitment parameters, the two distortion parameters, *g* and υ_{0}, were estimated without ambiguity. There were three reasons for this: *1*) the model was formulated such that the distortion component (Eq. 21) was independent of the recruitment component; *2*) within the distortion component, *g* and υ_{0} had very different roles in determining distortion behavior, and they did not combine in a nonlinear way; and *3*) distortion behavior, i.e., *phases 1* and *2* of the ΔP(*t*) responses, was a prominent part of the response that was largely separable from recruitment behavior and did not exhibit overt nonlinearities. Thus, in comparing *g* and υ_{0} estimated from hearts at different temperatures, *g* increased with increase in temperature, and υ_{0} did not, as expected.

To summarize, the RDIE model satisfied many requirements for descriptive validity but will need further refinement to allow unique estimation of parameters associated with its recruitment component.

Predictive validity depends on the model satisfactorily reproducing system behavior other than that used in data fitting/parameter estimation. Using parameters estimated from fitting to the ΔP(*t*) response family, the model successfully predicted isovolumic pressures and ΔP_{iso}(*t*) associated with beats at volumes above and below V_{BL}. Features of these model-predicted isovolumic pressures that corresponded to features of measured isovolumic pressures included: *1*) the dependence of isovolumic pressure amplitude on LV volume as expected from the Frank-Starling mechanism; *2*) slowed relaxation with increased pressure as a consequence of cooperative activation through *k*_{xb}(*t*); and *3*) reproduction of apparent isochronal-elastance lines at times corresponding to time at which rising pressure equaled half-maximal isovolumic pressure, time of maximal isovolumic pressure, and time at which falling pressure equaled half-maximal isovolumic pressure (Fig. 4). An important extension of the model's predictive capabilities will be its ability to reproduce realistic cardiac cycle events during simulated ejection into an arterial load. Such predictions will establish the potential utility of the RDIE model in situations where LV function is of primary interest.

Explanative uses of models enable causal explanation of system behavior. Thus the contribution of specific underlying mechanisms to global behavior can be identified, quantified, and explained through the use of the model. The RDIE model is strongly informed by and borrows its format from muscle contractile mechanisms. Thus the cycling of LV elastance elements between elastance and nonelastance states is analogous to the cycling of muscle XB between force-bearing and non-force-bearing states. The distortion of elastance elements with changing LV volume is analogous to the distortion of elastic XB with changing muscle length. Calcium activation of the global myocardium was treated no differently than if we were considering Ca^{2+} activation of a single sarcomere. The cooperative activation concept was taken from our earlier work on cooperative interactions within the myofilament system (10, 14). Thus, without attempting to make overly detailed connections, the RDIE model does build a bridge between underlying contractile mechanisms and global LV behavior, allowing explanations of LV function in terms of muscle contractile mechanisms.

Although many assumptions, approximations, and simplifications were needed to construct the model, the weakest assumption is that the LV myocardium is homogenous with respect to contractile activity and that circumferential chamber dimension changes are dominant. Good arguments for homogeneity in these isolated rat hearts may be made for early- and midsystole (1, 22, 25, 39). But as late relaxation periods are reached, heterogeneous contractile areas within the myocardium probably participate in global mechanical behaviors. For instance, the rebound in *phase 3* of late-cycle ΔP(*t*) responses that was poorly represented by the model could be due to a breakdown in uniformity in contractile activity within the myocardium, with resulting regionally distinct dimension changes. Such mechanisms are not a part of the RDIE model either explicitly or implicitly, and a model based on uniform contraction throughout the myocardium would never successfully capture these nonuniform events.

It is probable that the RDIE model will serve several purposes that call for application of an LV mechanodynamic model. But we believe that it will be most useful as a starting point for future model development. For instance, the RDIE model treats recruitment and distortion as completely independent processes, even though both processes must certainly share some common kinetic steps. Indeed, a single interconnected myofilament kinetic scheme accounts for both the kinetics of activation and the kinetics of XB cycling; some combination of steps is responsible for recruitment, and another combination is responsible for distortion. More work is required to elaborate how steps common to both processes affect their relative independence. At this time, however, the current version of the RDIE model serves as a model that will stand alongside other proposed LV models [such as the Burkhoff model (3, 9, 42, 50)] in being useful for analyzing certain kinds of LV mechanodynamic data, and further, it will be useful as a starting point for developing an improved LV model.

### Conclusion

The LV possesses an instantaneous elastance property that undergoes transient changes in volumetric distortion when there are changes in LV volume. The dynamics associated with elastance and its transient distortion arise from kinetic processes of myocardial contraction. Elastance element recruitment is driven both by myofilament activation and by LV volume. Elastance element distortion is driven by the time-rate-of-change of volume, and this confers velocity-dependent behavior to the ventricle. There is a cooperative form of activation [symbolized as *k*_{xb}(*t*)] that favors additional recruitment. This cooperative activation strongly influences recruitment dynamics during late systole and is responsible for sustaining systole for longer periods than is possible with Ca^{2+} activation alone. Coupling between distortion and recruitment confers at least three additional nonlinear effects: *1*) exaggerated early-cycle ΔP(*t*) recruitment responses; *2*) diminished late-cycle ΔP(*t*) recruitment responses; and *3*) directional differences (asymmetry) in the ΔP(*t*) response pattern. These nonlinear couplings superimpose on the uncoupled recruitment and distortion processes to integrate the separate processes responsible for overall dynamics. It can be hypothesized that dynamical integration through nonlinear coupling among the recruitment processes and between distortion and recruitment result in an internally ‘tuned’ dynamical system. A change in pattern would not only reflect a change in mechanical properties but also a change in the tuned coordination of the underlying kinetic processes from which whole-organ mechanodynamics emerge. Thus the pattern of LV responses to step-like changes in volume and the model that represents this pattern may serve as a sensitive indicator of coordination or lack thereof in underlying mechanisms responsible for LV function.

## GRANTS

This work was supported in part by a grant from the American Heart Association, Pacific Mountain Affiliate (0555548Z) to K. B. Campbell; by a generous gift from the Washington State Fraternal Order of Eagles, Pasco Erie; and by a PHS grant (RO1 HL-61487/62881) to H. L. Granzier.

## Footnotes

↵1 The true meaning of the word “dynamic” is important to this argument. We use the definition that a system is dynamic if its behavior evolves with time according to the history of that behavior as well as according to system inputs that drive the behavior. Thus the state of the system at previous times is important in determining the state of the system at the current time. Such history-dependent behavior in dynamic systems is traditionally accounted for by differential equation descriptions of the system.

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 © 2008 the American Physiological Society