Insights From Computational Modeling Into the Contribution of Mechano-Calcium Feedback on the Cardiac End-Systolic Force-Length Relationship

In experimental studies on cardiac tissue, the end-systolic force-length relation (ESFLR) has been shown to depend on the mode of contraction: isometric or isotonic. The isometric ESFLR is derived from isometric contractions spanning a range of muscle lengths while the isotonic ESFLR is derived from shortening contractions across a range of afterloads. The ESFLR of isotonic contractions consistently lies below its isometric counterpart. Despite the passing of over a hundred years since the first insight by Otto Frank, the mechanism(s) underlying this protocol-dependent difference in the ESFLR remain incompletely explained. Here, we investigate the role of mechano-calcium feedback in accounting for the difference between these two ESFLRs. Previous studies have compared the dynamics of isotonic contractions to those of a single isometric contraction at a length that produces maximum force, without considering isometric contractions at shorter muscle lengths. We used a mathematical model of cardiac excitation-contraction to simulate isometric and force-length work-loop contractions (the latter being the 1D equivalent of the whole-heart pressure-volume loop), and compared Ca2+ transients produced under equivalent force conditions. We found that the duration of the simulated Ca2+ transient increases with decreasing sarcomere length for isometric contractions, and increases with decreasing afterload for work-loop contractions. At any given force, the Ca2+ transient for an isometric contraction is wider than that during a work-loop contraction. By driving simulated work-loops with wider Ca2+ transients generated from isometric contractions, we show that the duration of muscle shortening was prolonged, thereby shifting the work-loop ESFLR toward the isometric ESFLR. These observations are explained by an increase in the rate of binding of Ca2+ to troponin-C with increasing force. However, the leftward shift of the work-loop ESFLR does not superimpose on the isometric ESFLR, leading us to conclude that while mechano-calcium feedback does indeed contribute to the difference between the two ESFLRs, it does not completely account for it.

In experimental studies on cardiac tissue, the end-systolic force-length relation (ESFLR) has been shown to depend on the mode of contraction: isometric or isotonic. The isometric ESFLR is derived from isometric contractions spanning a range of muscle lengths while the isotonic ESFLR is derived from shortening contractions across a range of afterloads. The ESFLR of isotonic contractions consistently lies below its isometric counterpart. Despite the passing of over a hundred years since the first insight by Otto Frank, the mechanism(s) underlying this protocol-dependent difference in the ESFLR remain incompletely explained. Here, we investigate the role of mechano-calcium feedback in accounting for the difference between these two ESFLRs. Previous studies have compared the dynamics of isotonic contractions to those of a single isometric contraction at a length that produces maximum force, without considering isometric contractions at shorter muscle lengths. We used a mathematical model of cardiac excitation-contraction to simulate isometric and force-length work-loop contractions (the latter being the 1D equivalent of the whole-heart pressure-volume loop), and compared Ca 2+ transients produced under equivalent force conditions. We found that the duration of the simulated Ca 2+ transient increases with decreasing sarcomere length for isometric contractions, and increases with decreasing afterload for work-loop contractions. At any given force, the Ca 2+ transient for an isometric contraction is wider than that during a work-loop contraction. By driving simulated work-loops with wider Ca 2+ transients generated from isometric contractions, we show that the duration of muscle shortening was prolonged, thereby shifting the work-loop ESFLR toward the isometric ESFLR. These observations are explained by an increase in the rate of binding

INTRODUCTION
The end-systolic force-length (or pressure-volume) relationships (ESFLRs) occupy a special place in cardiovascular physiology, often used as a proxy for cardiac contractility. Two distinct relationships are commonly reported, one for isometric (or isovolumic) contractions and a separate one for isotonic (or ejecting) contractions. It is sobering to reflect that the existence of these two distinct ESFLRs has been recognized since the 1899 publication by Otto Frank (1899) but remains incompletely explained. Frank's findings arose from experiments performed on the isolated heart of the frog and clearly reveal the existence of two distinct protocol-dependent end-systolic relations. Despite Frank's insight, and more than 100 years of research since, the mechanisms responsible for the protocol-dependence of the ESFLRs remain incompletely understood. This absence of a definitive explanation stands in stark contrast to the many experimental affirmations of the phenomenon, commencing with Brady (1967) and Gibbs et al. (1967), who independently demonstrated that the isometric ESFLR of isolated papillary muscles lies to the left of the equivalent isotonic (or workloop) curve.
This phenomenon can be readily replicated and a variety of explanations have been proffered. These include the contribution of strain rate-dependent transitions between strong and weak cross-bridge conformations (Landesberg, 1996), velocitydependent shortening deactivation of thin filaments (Iribe et al., 2014), and sarcomeric force-length dynamic interactions (Pironet et al., 2013). Perhaps the most germane of which was that of Brady (1967): "the time-course of the active state, " which focused attention on the duration of the Ca 2+ transient. But it was Housmans et al. (1983) who demonstrated that active shortening retards the decline of the Ca 2+ transient, a phenomenon that their group boldly attributed to the affinity of Ca 2+ for troponin-C, hinting at the role of mechano-calcium feedback. Their insight has been corroborated by many others, in both intact (ter Keurs et al., 1980;Lab et al., 1984;Kentish and Wrzosek, 1998) and skinned (Hibberd and Jewell, 1982) cardiac muscle preparations, and in both relaxed and rigor states, thereby demonstrating the requirement that cross-bridges be attached for the phenomenon to instantiate (Hofmann and Fuchs, 1987;Kurihara et al., 1990). The encompassing phrase "myofilament length-dependent activation" includes the now widely accepted notion of force-dependent binding of Ca 2+ to troponin (Rice and de Tombe, 2004;de Tombe et al., 2010;Solaro and Kobayashi, 2011;de Tombe and ter Keurs, 2016;Li et al., 2016), where force-dependence arises primarily from the consequence of length-dependent overlap of the thick and thin contractile filaments. At the core of this mechanism is the role of Ca 2+ as an activator of contraction and as a target for force feedback-the mechano-calcium mechanism.
A shortcoming of previous experimental studies has been the comparison of isotonic contractions to an isometric contraction at the unique length (L o ) that maximizes force (Housmans et al., 1983;Lab et al., 1984). While such studies have uncovered a number of key properties of actively shortening muscle, the failure to compare isotonic contractions to isometric contractions of equivalent force has hindered the development of a satisfactory explanation to explain their different ESFLRs.
This perceived shortcoming has prompted us to investigate the role of mechano-calcium feedback in the mode-dependence of the ESFLR and to ask whether the Ca 2+ transients of isometricand work-loop-derived ESFLRs are different under equivalent force conditions. If that were to transpire, then a second question would arise: could the activation of work-loop contractions using Ca 2+ transients derived from isometric contractions reconcile these two ESFLRs? To that end, we undertook the assembly of a mathematical model of cardiac excitation-contraction. The model tightly couples three existing components-one for each of electrical excitation (Rogers and McCulloch, 1994), Ca 2+ dynamics (Hinch et al., 2004) and cross-bridge mechanics (Tran et al., 2017). The latter module, which is based on a three-state model of cross-bridge cycling (Rice et al., 2008), is capable of generating both isometric and afterloaded isotonic contractions under either fixed or variable Ca 2+ transients. It thus provides a tool with which to address the two questions raised above.

Model Development
A novel integrated model of cardiac cellular excitation-contraction was developed to investigate the protocoldependence of the isometric and work-loop ESFLRs. The model was designed to encapsulate sufficient biophysical detail for the purposes of the current study, with a focus on computational efficiency to permit integration into large scale tissue modeling studies.
The Tran et al. (2017) mechano-energetics model forms the basis of the coupled model, since it captures the detailed biophysical interdependencies of thin filament activation and cross-bridge cycling. It has been previously parameterized to simulate a wide range of experimentally observed cardiac myofilament behavior, including steady-state, force-length, and force-Ca 2+ relations, as well as dynamic force transients and  Rogers and McCulloch (1994) action potential model. Ca-Trop represents Ca 2+ bound to troponin-C, which is modulated by active force, and therefore mediates the mechano-calcium feedback pathway. metabolite dependence (Rice et al., 2008;Tran et al., 2010). The force-dependent binding of Ca 2+ to troponin-C is also explicitly captured allowing cross-bridge force production to modify the intracellular Ca 2+ transient. The total force generated by the model is composed of an active and a passive component. To enable the required range of simulations, the cross-bridge model was coupled to the Hinch et al. (2004) model of Ca 2+ -induced-Ca 2+ -release (CICR), which simulates Ca 2+ fluxes across the sarcolemmal and sarcoplasmic reticular membranes. The CICR model was driven by a simplified model of the cardiac action potential. Since our study did not require detailed cellular electrophysiology, we parameterized the phenomenological Rogers and McCulloch (1994) model to simulate a rat cardiac action potential and used it to drive the cellular contraction cycle. The mechano-calcium feedback between the cross-bridge and CICR models was implemented using equations 56-64 of Rice et al. (2008). All three components of the model were previously parameterized for rat cardiomyocyte, and were not modified from their original values in this study. The Rogers and McCulloch (1994) action potential model is not temperature dependent, while the Hinch et al. (2004) model was parameterized using data at 22 • C and the Tran et al. (2017) cross-bridge model is appropriate over the temperature range of 23-37 • C. We, therefore, set the temperature of the Tran et al. (2017) cross-bridge model to 23 • C to be consistent with that of Hinch et al. (2004). A schematic of the integrated Tran-Hinch-Rogers (THR) model is presented in Figure 1. The simulated calcium transient therefore arises from interactions between the dynamics of the CICR model and the mechano-calcium feedback of the cross-bridge model.
Each of the Rogers and McCulloch (1994); Hinch et al. (2004), and Tran et al. (2017) models is available in the Physiome Model Repository. 1 We took advantage of the modular features of the CellML standard (Cuellar et al., 2003) and the software tool OpenCOR (Garny and Hunter, 2015) to integrate these disparate models to construct the coupled THR model of cardiac excitation-contraction (see Nickerson and Buist, 2008;Terkildsen et al., 2008 for detailed demonstrations of such model integration). The CellML models, descriptions of simulation experiments, and parameters for reproducing the figures in this study are available from: https://models.physiomeproject.org/ workspace/4ca.

Isometric and Work-Loop Protocols
The THR model was used to simulate isometric and workloop contractions. Isometric simulations consisted of holding the sarcomere length constant throughout the duration of a contraction. Initial sarcomere lengths were varied between 1.95 and 2.3 µm, the length (L o ) at which maximum force is produced. At each sarcomere length, twitches were simulated at 1 Hz frequency until a beat-to-beat steady-state was reached (20 beats). The isometric ESFLR was derived by plotting the peak force as a function of sarcomere length. All forces in the simulations were normalized to the peak L o -isometric force.
The work-loop protocol was designed to mimic the wholeheart pressure-volume loop (Taberner et al., 2011), and was simulated using four contraction phases.

Phase 1: Isometric Contraction
Starting at an initial sarcomere length of L o , stimulation by an action potential initiates CICR and force generation. Muscle length is held constant as force develops. Isometric contraction continues until the total muscle force (the sum of active and passive forces) exceeds a user-specified afterload.

Phase 2: Isotonic Contraction
Following Phase 1, once the muscle force exceeds the afterload, the muscle shortens isotonically against the constant afterload. The shortening velocity is calculated using Equation 38 from Rice et al. (2008). The simulation continues in this phase for as long as the force generated by the muscle matches or exceeds the userspecified afterload.

Phase 3: Isometric Relaxation
When the force generated by the muscle falls below that of the afterload, it is declared to have reached end-systole. The muscle length is held constant while the force decreases until Ca 2+ declines to diastolic levels.

Phase 4: Re-stretch
When Ca 2+ reaches diastolic levels, the now-quiescent muscle is re-stretched, at constant velocity, to its initial length of L o .
Work-loops were simulated using afterloads between 0.2 and 1 (normalized to isometric force at L o ). Work-loops at a normalized afterload of 1 produced contractions that were identical to an L o -isometric contraction, because the afterload exceeded the threshold for muscle shortening. Similar to the isometric protocol, a train of work-loop contractions was generated to attain steady-state at a frequency of 1 Hz. The workloop-derived ESFLR was given by the force and length of the muscle at the end of Phase 2 (end-systole) for each work-loop FIGURE 2 | Comparison of THR model simulations with experimental quick-release data. Experimental quick-release shortening of ferret papillary muscle by Kurihara and Komukai (1995) (left) at 30 • C and equivalent simulations from the THR model (right). Each of the four panels for the experimental (A1-D1) and model (A2-D2) outputs represent different timings of the quick release. Within each panel are traces of outputs from isometric (black lines) and quick-release protocols (gray lines). From top to bottom, the traces represent muscle length, intracellular Ca 2+ concentration, total force development, and difference of Ca 2+ between the two protocols. (A1-D1 were reproduced with permission from John Wiley and Sons).
contraction. The sarcomere lengths for the isometric contractions and afterloads for the work-loop contractions were chosen so that they formed pairs of force-equivalent contractions. This allowed the comparison of force profiles and Ca 2+ dynamics for the two modes of contraction at the same level of force production.

Model Validation
We used the coupled THR model to examine the role of Ca 2+ in the ESFLR of cardiac muscle contraction. The Ca 2+ and force dynamics of the model were validated by comparing simulated output to quick-release shortening experiments reported by Kurihara and Komukai (1995). In those experiments, isolated cardiac muscles were initially held at a length (L o ) that produced maximum force, and then was rapidly shortened to a new length (0.92 L o ) at four different predetermined times following the stimulus. The resulting force and Ca 2+ transients were measured and compared to those of an isometric contraction where the length was held constant at L o throughout the duration of the twitch (Figure 2). The sudden length change due to the quickrelease caused an immediate drop in developed force and a near-simultaneous increase of intracellular Ca 2+ , evidenced by the extra surge of Ca 2+ during the declining phase of the Ca 2+ transient. The magnitude of the extra Ca 2+ elicited by the quickrelease twitch is shown in the lower-most trace in each panel of Figure 2.
We simulated the quick-release protocol in our THR model and semi-quantitatively reproduced the experimental data from Kurihara and Komukai (1995). Note that none of the parameters in the THR model were optimized to reproduce the experimental data. The quick release from L o to 0.92 L o caused a decrease in force and an increase in intracellular Ca 2+ . The rates of force redevelopment in the model simulations were greater than those of the experiments, as is evident in panels C1 and D1 (experiment) versus C2 and D2 (model) of Figure 2, where the length change was applied after the peak of the Ca 2+ transient. This may be due to the inherently faster cross-bridge cycling rates of the rat (model) versus those of other larger mammalian species, such as the ferret (Hasenfuss et al., 1991;Palmer and Kentish, 1998). The Ca 2+ surges in panels C2 and D2 were smaller than those recorded in the experiments (panels C1 and D1), because the recovery force was much larger. Nevertheless, the force recovery and Ca 2+ surge were semi-quantitatively similar to the experimental data in panels A1 and B1, thereby validating the ability of the THR model to mimic force-dependent binding of Ca 2+ from troponin-C and the consequent effect on the simulated Ca 2+ transient.

Simulating Isometric and Work-Loop ESFLRs Using Dynamic Ca 2+
The ESFLRs for isometric and work-loop contractions were simulated using the THR model following the protocols outlined in Section "Materials and Methods" (Figure 3A). The isometric ESFLR was derived from isometric contractions at different Frontiers in Physiology | www.frontiersin.org sarcomere lengths while the work-loop ESFLR was derived from afterloaded shortening contractions, all of which commenced at the same initial length (L o ). These simulation results compare well with the experimental data from rat trabeculae (Tran et al., 2016) in Figure 3B and are consistent with other previous studies on rat cardiac myocytes (Iribe et al., 2014), rat trabeculae (Han et al., 2014a,b), rabbit papillary (Sørhus et al., 2000), ferret papillary (Hisano and CooperIV, 1987), and rat papillary (Gülch, 1986), the simulated work-loop ESFLR lies to the right of its isometric counterpart. In Figure 3A, the afterloads for the simulated work-loop contractions (black circles) correspond to isometric contractions of equivalent peak force (gray circles). For a given afterload, a work-loop contraction did not shorten to the sarcomere length of its force-equivalent isometric contraction.
Previous studies have attributed this difference to mechanocalcium feedback, a phenomenon whereby active shortening of a muscle during a contraction reduces the affinity of Ca 2+ for troponin-C, thus offloading Ca 2+ into the cytosol and reducing the level of activation of the cross-bridges (Housmans et al., 1983;Lab et al., 1984;Baan et al., 1992). The principal mechanism for this phenomenon is attributed to the force-dependent binding of Ca 2+ to the troponin-C complex, whereby the affinity of Ca 2+ for troponin-C increases with increasing force (Hibberd and Jewell, 1982;Hofmann and Fuchs, 1987). For our simulated workloop contractions, shortening during the isotonic phase (Phase 2) is expected to offload Ca 2+ from troponin-C and enhance the descending portion of the intracellular Ca 2+ transient. The level of enhancement is also predicted to increase with decreasing afterload.
Model simulations of the Ca 2+ transients associated with work-loop contractions are presented in Figure 4. The workloop contractions (black lines) did indeed prolong the Ca 2+ transient compared to the isometric contraction at L o . The Ca 2+ transients widened with decreasing afterload, but the extent of widening was small relative to the effect of decreasing sarcomere length on the isometric contractions (Figure 4, gray lines). The narrowest of the isometric Ca 2+ transients corresponds to an isometric contraction at L o (gray line in Figure 4 inset). Note that because all of the work-loop contractions originated at L o , the Ca 2+ transients for the work-loop contractions matched those of the L o -isometric contraction up to the onset of Phase 2 (i.e., up to t = 45 ms). Decreasing the muscle length of the isometric contractions resulted in a more significant widening of the Ca 2+ profile relative to the effect of decreasing afterload in the work-loop contractions. For clarity, the Ca 2+ transients for the L o -isometric contraction and the work-loop contraction at the lowest afterload are shown in the inset of Figure 4. Additional simulations of work-loop ESFLRs at initial sarcomere lengths below L o were consistent with recent experimental data on rat cardiac trabeculae showing the preload dependence of work-loop contractions , and simulation of associated Ca 2+ transients were consistent with the trends observed in Figure 4 (Supplementary Material).
The variation of the Ca 2+ transients within and between the two modes of contraction (Figure 4) is a result of the interaction between cross-bridge force production and the force-dependent affinity of Ca 2+ for troponin-C. In Figure 5, the intracellular Ca 2+ concentration, force, and the flux of Ca 2+ binding to troponin (Ca-troponin) as functions of time are presented for the L o -isometric contraction, as well as for a selected pair of isometric and work-loop contractions at equivalent peak forces. The forceequivalent isometric contraction was simulated at SL = 0.91 L o , while the force-equivalent work-loop was simulated at a normalized afterload of approximately 0.45. The narrowest of the three Ca 2+ transients was generated by the L o -isometric contraction, while the widest was generated by the 0.91 L oisometric contraction. The variations in the Ca 2+ profiles were directly due to the differential responses of the Ca-troponin fluxes to the contractions.
The Ca-troponin fluxes began diverging at t = 20 ms when the active force in each of the three contractions started to rise ( Figure 5; t 1 ). The faster rate of force development in the L oisometric and work-loop contractions transiently increased the Ca-troponin flux above that of the 0.91 L o -isometric contraction, which had a lower rate of force development over that initial period. As a result, less Ca 2+ was bound to troponin for the 0.91 L o -isometric contraction, which was reflected in the widening of its Ca 2+ transient. Intracellular Ca 2+ concentration directly affects the level of cross-bridge activation and the ability of the muscle to maintain force. For given equivalent peak forces, the model simulations reveal that isometric contractions were characterized by wider Ca 2+ transients, relative to those of workloop contractions (Figure 4). We hypothesized that driving workloop contractions with wider Ca 2+ transients would result in greater isotonic shortening and a leftward shift of the work-loop derived ESFLR toward the isometric ESFLR.
Simulating Work-Loop ESFLRs Using Fixed Ca 2+ Transients Figure 6 shows simulation results whereby Ca 2+ transients associated with isometric contractions (Figure 4) were used to drive force-equivalent work-loop contractions. The activation of work-loop contractions using the wider Ca 2+ transients derived from isometric contractions caused a leftward-shift of the work-loop derived ESFLRs (Figure 6). Specifically, the greater level of activation increased the proportion of cross-bridges in forceproducing states, leading to greater force production at lower sarcomere lengths.
The extent of the shift increased with decreasing afterloads. At high afterloads, the wider Ca 2+ transients were not sufficient to shift the work-loop end-systolic point to coincide with its isometric counterpart. However, at afterloads below approximately 0.4 normalized total force, the work-loop endsystolic point was shifted beyond its isometric counterpart. The two ESFLRs therefore crossed over at around the midpoint. These simulations demonstrate that the work-loop derived ESFLR can be shifted toward the isometric ESFLR by increasing the duration of the Ca 2+ transient to prolong the shortening phase. But the contraction-mode difference in the Ca 2+ transients cannot completely explain the discrepancy between the isometric and work-loop ESFLRs.
FIGURE 5 | Simulations highlighting the mechano-calcium feedback mechanism. Simulated intracellular Ca 2+ transient (A), total force production (B), and Ca-troponin flux (C) for an L o -isometric contraction, and a pair of isometric (SL = 0.91 L o ) and work-loop (normalized afterload at 0.45) contractions at the same peak force. Inset: close-up of the Ca-troponin flux at the start of the twitch. The dotted vertical lines indicate the times at the start of force production (t 1 ) and at the start of the isotonic phase in the work-loop contraction (t 2 ).

DISCUSSION
Using computational modeling, we have demonstrated that mechano-calcium feedback in cardiac muscle produces Ca 2+ transients, which are dependent on the mode of contraction. We compared the force and Ca 2+ dynamics of work-loop contractions with their force-equivalent isometric counterparts ( Figure 5) and revealed that the Ca 2+ transients generated by isometric contractions were of longer duration (Figure 4). Our computational modeling approach enabled the application of a novel "cross-over" technique wherein work-loops are driven with Ca 2+ transients generated by force-equivalent isometric contractions (Figure 6), a feat which has yet to be realized experimentally. However, the difference in the duration of Ca 2+ transients is insufficient to completely explain the protocoldependent difference between the isometric-and work-loopderived ESFLRs. At higher afterloads, the increase in the Ca 2+ transient duration was insufficient to overcome the shortening deactivation that is present in the cross-bridge model (Rice et al., 2008). Specifically, muscle shortening reduces cross-bridge distortion, which promotes cross-bridge detachment and leads to a loss of activation and a reduction in tension. At lower afterloads, the Ca 2+ transients generated by force-equivalent isometric contractions were longer, overcompensating for the shortening deactivation, resulting in a shift of the work-loop endsystolic points past those generated from isometric contractions.
The mechanism of force-dependent affinity of Ca 2+ to troponin-C is responsible for the prolongation of the Ca 2+ transient duration for both modes of contraction (Figure 5). Model simulations revealed that reducing the initial muscle length of an isometric contraction produced a much more significant widening of the Ca 2+ transient than decreasing the afterload for a work-loop contraction. Isometric contractions therefore operate at a greater level of Ca 2+ activation, explaining the leftward-shift of the isometric ESFLR. Our conclusion is in accord with the study of De Tombe and Little (1994) on rat cardiac trabeculae at 27 • C, which showed that the ESFLR of shortening contractions is left-shifted with respect to isosarcometric contractions when the duration of activation is prolonged by substitution of Ca 2+ with Sr 2+ . Our simulations reveal that the difference in the ESFLR of these two modes of contraction cannot be solely explained by the differing levels of activation elicited through the mechanism of force-dependent affinity of Ca 2+ binding to troponin-C.

Force-Equivalent Contractions
Previous experimental studies investigating the protocoldependence of the ESFLR compared isotonic contractions to a single isometric contraction at L o (Housmans et al., 1983;Lab et al., 1984). Housmans et al. (1983) initially hypothesized that: "The rightward location of work-loop end-systolic points is due to sarcomere shortening abbreviating the Ca 2+ transient responsible for excitation-contraction coupling" but rejected this hypothesis when their experimental data showed that the isotonic contraction actually exhibited a wider Ca 2+ transient than the isometric case. Instead, they proffered the hypothesis that the "extra" Ca 2+ released during the relaxation phase of the Ca 2+ transient appears too late in the twitch cycle to have an effect on prolongation of force development.
Our simulation results are consistent with these two studies provided that we replicate their protocol by comparing the Ca 2+ transient from an afterloaded work-loop contraction to that of a single isometric contraction at L o (Figure 4 inset). The Ca 2+ transient of the work-loop contraction was indeed wider than that of the L o -isometric contraction. Compared to the experimental results of Housmans et al. (1983) and Lab et al. (1984), the magnitude of the increase of Ca 2+ in our simulations was less pronounced. This is likely to be because the lowest afterload of our work-loop simulations is not as low as those achieved in the experiments owing to the lower passive force component in experimental preparations relative to that of our model. But this reveals only part of the story, since we showed that the greater prolongation of the Ca 2+ transient for a work-loop contraction holds only when compared to the L o -isometric contraction. For isometric contractions at reduced muscle lengths, the predicted Ca 2+ transient is of a longer duration than that of the work-loop contractions (Figure 4). In other words, a work-loop does have an abbreviated Ca 2+ transient when compared to an isometric contraction that is of equivalent peak force (Figure 5).

Time-Course of Isometric and Work-Loop Twitches
The end-systolic point of an isometric contraction coincides with the peak force of the twitch, whereas the end-systolic point of a work-loop contraction is defined as the point where the muscle can no longer shorten. As a result, the end-systolic point of a work-loop contraction always occurs later in the time-course of the Ca 2+ transient when its intracellular Ca 2+ concentration is lower (Figure 5). Therefore, the width of the Ca 2+ transient, governed by the rate at which it declines, is a key factor that determines the positioning of the work-loop derived ESFLR. By widening the Ca 2+ transient, more time is available for shortening to occur before Ca 2+ levels fall too low (Figure 6).
The different time-courses of isometric and work-loop contractions also explain their varying effects on the width of the Ca 2+ transient. Decreasing the afterload of a workloop contraction has a relatively small impact on the Ca 2+ transient compared to decreasing the initial length of an isometric contraction (Figure 4). Every work-loop starts Phase 1 at an initial length of L o and behaves like an L o -isometric contraction during Phase 1. The mechanical feedback of isotonic shortening onto the affinity of Ca 2+ for troponin-C does not occur until Phase 2, which is approximately 45 ms after the start of the Ca 2+ transient (Figure 5; t 2 ). This time point is well into the declining phase of the Ca 2+ transient; hence, the enhancement of the transient from this point on is relatively small. In contrast, for an isometric contraction, the muscle length is changed prior to stimulation of the muscle. The commencement of the twitch takes place 25 ms after the start of the Ca 2+ transient (Figure 5; t 1 ). Hence, for isometric contractions at reduced initial muscle lengths, the mechanical feedback begins to affect the Ca 2+ transient from an earlier time point, leading to a greater increase in the width of the Ca 2+ transient relative to decreasing the afterload in work-loop contractions. Mechanical feedback from isometric contractions arises from the positive relationship between muscle length and the initial rate of tension development. At decreasing muscle lengths, the rate of rise of tension development also decreases (Han et al., 2014a,b;Tran et al., 2016), which lowers the Ca-troponin flux, allowing enhancement of the entire declining phase of the Ca 2+ transient (Figure 5).

Time-Course of Ca 2+ Transients
The duration of our simulated Ca 2+ transients increased with decreasing sarcomere length and decreasing afterload for isometric and work-loop contractions, respectively (Figure 4). The prolongation of the Ca 2+ transient associated with isometric contractions at shorter muscle lengths is consistent with experimental studies on isolated rat cardiac tissues (Allen and Kurihara, 1982;Kentish and Wrzosek, 1998). At shorter muscle lengths, the affinity of Ca 2+ for troponin-C is lower, hence less Ca 2+ is bound to troponin-C, resulting in a wider Ca 2+ transient. The effect of an afterloaded contraction on the Ca 2+ transient has typically been compared only to an L o -isometric contraction, as has been discussed above.
Despite a prolongation of the Ca 2+ transient due to either a reduction of initial muscle length in isometric contractions or a reduction of afterload in work-loop contractions, the peaks of the simulated Ca 2+ transients remained relatively constant for both modes of contraction. Studies on rat cardiac trabeculae have reported no change in the peak Ca 2+ transient in isometric contractions at different muscle lengths (Allen and Kurihara, 1982;Backx and Ter Keurs, 1993;Kentish and Wrzosek, 1998). Studies on the effect of a shortening contraction on the peak of the Ca 2+ transient are not so clear. Peak Ca 2+ has been shown to be higher for a shortening contraction compared to an isometric contraction in rat cardiac myocytes (Yasuda et al., 2003) and trabeculae (Janssen and de Tombe, 1997). In contrast, no differences were reported for either guinea pig myocytes (White et al., 1995), or for papillary muscles of either ferret (Lab et al., 1984) or cat (Housmans et al., 1983). In rat cardiomyocytes undergoing auxotonic contractions, Hongo et al. (1996) showed no change in the Ca 2+ transient peak immediately after a stretch, but a gradual increase in peak force and peak Ca 2+ transient 3 min later.
In our simulations, the peak of the Ca 2+ transient remained relatively constant, because only a fraction of peak force development had occurred at the time of peak Ca 2+ (Figure 5). The completion of force development occurs after the peak Ca 2+ has passed and is reflected in the prolongation of the declining phase of the Ca 2+ transient (Allen and Kurihara, 1982).

Comparison to Previous Computational Modeling Studies
A number of studies have developed computational approaches to investigate the ESFLR of cardiomyocytes (Iribe et al., 2014) and the end-systolic pressure-volume relation (ESPVR) of the whole heart (Landesberg and Sideman, 1994;Landesberg, 1996;Shimizu et al., 2002). Iribe et al. (2014) were able to reconcile their protocol-dependent ESFLR data collected from rat cardiomyocytes by modifying an excitation-contraction model to include a term that represents the shortening velocity-dependent inactivation of the thin filament. However, despite the absence of this mechanism in our THR model, we were still able to produce the contraction-mode dependence of the ESFLR, suggesting that velocity-dependent inactivation of thin filaments is not the sole mechanism. In our THR model, the protocol-dependence of the ESFLR was attributed to force-dependent affinity of Ca 2+ to troponin-C, a mechanism that is in accord with that of many other modeling studies that describe thin filament activation and cross-bridge force generation (Katsnelson et al., 1990;Izakov et al., 1991;Markhasin et al., 2003;Niederer et al., 2006;Syomin and Tsaturyan, 2017).
Our model simulations of the mode-dependency of the ESFLR (Figure 3) and the effect of changing initial sarcomere length on the Ca 2+ transient from isometric contractions (Figure 4) are consistent with the model simulations of Landesberg and Sideman (1994), which incorporated force-dependent affinity of troponin-C for Ca 2+ as well as velocity-dependent transitions between weak and strong cross-bridge states in their fourstate model of cross-bridge contraction. In subsequent work, Landesberg (1996) showed that the ESPVR of work-loop contractions could be shifted by varying the scale of the velocitydependence parameter. However, in that model, the effect of changing initial sarcomere length on the Ca 2+ transient from isometric contractions was reversed, i.e., longer initial sarcomere lengths resulted in wider, rather than narrower, Ca 2+ transients. We are unable to determine why the behavior of the Ca 2+ transient is inconsistent between those two publications by Landesberg except to point out that in Landesberg (1996), a simple thin-wall model of the LV was added.
In experiments on canine hearts, Shimizu et al. (2002) reported a single ESPVR for isovolumic and ejection contractions and invariant shape of the Ca 2+ transient for both modes of contraction. As such, they used a single Ca 2+ transient to drive their model of cross-bridge force development. They concluded that length-dependence of the Ca 2+ binding affinity for troponin-C was not an important factor in the canine heart, which is at odds with our findings using the THR model of rat myocardium. This disparity is likely due to the peculiarity of canine hearts, which have been reported to exhibit an ESPVR that is independent of contraction protocol (Taylor et al., 1969;Suga and Sagawa, 1974;Weber et al., 1976;Suga et al., 1981). Contraction-mode independence of the ESPVR is likely an artifact arising from performance of work-loop contractions at low diastolic left ventricular volumes . These findings stand in contrast to those from the hearts of other smaller mammalian species (ferret, rabbit, and rat) (Hisano and CooperIV, 1987;Sørhus et al., 2000;Han et al., 2014a,b;Tran et al., 2016), all of which show contraction-mode dependence.

Model Limitations
The THR model of excitation-contraction utilizes a simple phenomenological model of rat cardiac electrophysiology to capture the action potential waveform, and couples it to detailed biophysical models of Ca 2+ handling and cross-bridge kinetics. The phenomenological nature of the action potential waveform means that electrophysiological length-dependent pathways, such as stretch-activated channels (SACs) which underlie the slow force response (Calaghan and White, 2004;Niederer and Smith, 2007;Dowrick et al., 2019), are not explicitly captured. But despite their absence, the THR model is able to reproduce the widely observed difference in isometric and work-loop ESFLRs, suggesting that length-dependent electrophysiological pathways play a minor role, if any, in eliciting this phenomenon.
The contraction-mode dependence of the cardiac ESFLR has been observed to be independent of experimental temperature (Gülch, 1986;Hisano and CooperIV, 1987;Sørhus et al., 2000;Han et al., 2014a,b;Tran et al., 2016), providing confidence in the predictions of our model, which was parameterized using data at room temperature. Experimental reports from cat papillary (Downing and Sonnenblick, 1964) and frog myocytes (Parikh et al., 1993) at room temperature, which purported to show contraction-mode independence of the ESFLR, was recently refuted in a re-examination of the data using a novel end-systolic zone framework .

Implications for Cardiac Physiology
The mechanical loading experienced by individual cardiomyocytes in the heart wall during a contraction cycle is heterogeneous. The tension-time-courses generated by cardiomyocytes vary transmurally from endocardial to the epicardial surfaces (Khokhlova et al., 2020). Mechanical feedback mechanisms allow communication of, and adaptation to, changing loading conditions across individual cells and over regions of myocardium. In the case of mechano-calcium feedback, facilitated by binding of calcium to troponin C, this mechanism allows modulation of twitch dynamics in response to changes in local loading conditions.
The impact of pathological changes on mechano-calcium feedback on cardiac function, in particular, on the ESFLR is not well known. Our recent work demonstrates that in healthy cardiac tissues, the end-systolic zone (defined as the region enclosed by isometric and work-loop ESFLRs; see Han et al., 2019) increases in the presence of isoproterenol-induced enhancement of contractility (Tran et al., 2020). We would expect the end-systolic zone to also be impacted in a disease setting. We note that mutations in the gene that encodes for cardiac troponin C (TNNC1) have been associated with hypertrophic cardiomyopathy (Liang et al., 2008;Pinto et al., 2009;Parvatiyar et al., 2012), restrictive cardiomyopathy (Parvatiyar et al., 2010), and ventricular fibrillation (Parvatiyar et al., 2012). While these pathologies are associated with the role of troponin C as an activator of contraction, changes in Ca 2+ binding