Implementation of Contraction to Electrophysiological Ventricular Myocyte Models, and Their Quantitative Characterization via Post-Extrasystolic Potentiation

Heart failure (HF) affects over 5 million Americans and is characterized by impairment of cellular cardiac contractile function resulting in reduced ejection fraction in patients. Electrical stimulation such as cardiac resynchronization therapy (CRT) and cardiac contractility modulation (CCM) have shown some success in treating patients with HF. Computer simulations have the potential to help improve such therapy (e.g. suggest optimal lead placement) as well as provide insight into the underlying mechanisms which could be beneficial. However, these myocyte models require a quantitatively accurate excitation-contraction coupling such that the electrical and contraction predictions are correct. While currently there are close to a hundred models describing the detailed electrophysiology of cardiac cells, the majority of cell models do not include the equations to reproduce contractile force or they have been added ad hoc. Here we present a systematic methodology to couple first generation contraction models into electrophysiological models via intracellular calcium and then compare the resulting model predictions to experimental data. This is done by using a post-extrasystolic pacing protocol, which captures essential dynamics of contractile forces. We found that modeling the dynamic intracellular calcium buffers is necessary in order to reproduce the experimental data. Furthermore, we demonstrate that in models the mechanism of the post-extrasystolic potentiation is highly dependent on the calcium released from the Sarcoplasmic Reticulum. Overall this study provides new insights into both specific and general determinants of cellular contractile force and provides a framework for incorporating contraction into electrophysiological models, both of which will be necessary to develop reliable simulations to optimize electrical therapies for HF.


Introduction
Heart failure (HF) affects over 5 million Americans and is characterized by impairment of cellular cardiac contractile function and reduced ejection fraction in patients [1]. Electrical stimulation such as cardiac resynchronization therapy (CRT) and cardiac contractility modulation (CCM) have shown success in treating patients with HF [2,3].
CRT is a current treatment for patients with congestive heart failure (CHF) for a wide QRS complex, where biventricular pacemakers are used to improve electrical synchrony and presumably ejection fraction [2,4]. However, the response to CRT varies greatly among patients. Improvements are quite variable with up to 30% of patients being non-responders to this treatment at all [2,5]. Nevertheless, the efficacy and optimization of CRT continues to being improved [6,7]. For example, it has been shown that extensive electrical remodeling is significantly associated with better survival rates after CRT [6]. Other measurements than electrical dyssynchrony have been investigated as identifications for CRT implantation estimation. It has been suggested that mechanical dyssynchrony is as well essential in prediction of CRT response [7].
An alternate electrical therapy for HF is cardiac contractility modulation (CCM) which delivers electric fields to the heart during its refractory period. Although its mechanism of action was initially thought to result from increased calcium flux from Sarcoplasmic Reticulum (SR) [3,8], recent studies suggest this is not the case, but that increased contractile function is the result of beta adregenergic stimulation [9]. Clinical studies have showed that CCM therapy improved quality of life, exercise capacity, New York Heart Association (NYHA) class, and ejection fraction (EF) during long-term follow up [10]. However, other clinical studies have reported contradictory evidence where CCM did not change effects on hospitalization or mortality [11], and it did not increase myocardial oxygen consumption [12].
It is clear then, that optimization of CRT and CCM therapies require an understanding of the effects of electric fields on myocyte dynamics especially for inotropicity as regulated by intracellular calcium. Computer simulations have the potential to improve these therapies through clarification of the underlying mechanisms [13]. For example, one challenge in CRT is to determine the optimal pacing locations; Miri et al. provided an optimization strategy to find the best pacing sites and timing delays in CRT based on biventricular-paced activation sequences and ECGs obtained from simulations using patient-specific anatomy and pathophysiology models [14]. Computer simulations can also be beneficial to help to predict the response of patients to certain therapies [15,16]. For instance, Niederer et al. have designed an electromechanical heart model based on clinical observations [17]. The model predicted that patients with dyssynchronous electrical activation but effective length-dependent tension regulation at the cellular scale are less likely to respond to CRT treatment than patients with attenuated or no length dependent of tension.
Such simulations of electrical therapy require a robust model of cellular excitation-contraction coupling that occurs in the cardiac myocyte such that the electrical and contraction predictions are both accurate. During the depolarization of an action potential, L-type Ca 2+ channels are activated and the influx of Ca 2+ current into the narrow dyadic space induces a large Ca 2+ release from the Junctional Sarcoplasmic Reticulum (JSR) through the Ryanodine Receptor (RyR). Intracellular Ca 2+ concentration is thus greatly increased immediately following the depolarization of the transmembrane action potential. Some of these Ca 2+ ions bind to affinity sites on Troponin C therefore enable myosin heads, which contain cross-bridges, to attach to actin. Myosin heads 'walk' on actin generating force, transforming chemical energy to mechanical energy resulting in cell contraction. Then during the repolarization of the action potential, Ca 2+ leaves the cytoplasmic compartment into either the Network Sarcoplasmic Reticulum (NSR) through the Ca-ATPase pump or out of the cell membrane though Na-Ca exchangers. The decreased intracellular Ca 2+ concentration causes Ca 2+ to detach from troponin C, inducing uncoupling of Myosin heads and actin, which ends the contraction process [18]. From this Ca 2+ handling process, we can see that Ca 2+ is deeply involved in the two highly non-linear systems, namely the electrical activation and tension generation and Ca 2+ binding with troponin C is the key component linking electrical signals to the activation of tension [18].
While there are close to a hundred cell models describing the detailed electrophysiology (EP) of cardiac cells, the majority of cell models do not include the equations to reproduce contractile force [19,20]. This is a considerable drawback. We suggest the most appropriate and challenging test for Electrical-contractile coupled (ECC) cell models derived from EP models is to reproduce the "postextrasystolic potentiation" (PESP) behavior of the heart [21]. PESP is an example of the changes in stimulation pattern on contractile strength. The effect can be demonstrated with the PESP pacing protocol [22]. First, a train of stimuli is delivered with a fixed basic cycle length called the priming period (PP). Second, an extrasystolic (ES) beat is delivered after the last priming stimulus by an interval called extrasystolic interval (ESI). Following the ES beat, a postextrasystolic (PES) beat is then delivered after another interval called the postextrasystolic interval (PESI). In normal hearts the strength of postrasystles is a function of both ESI and PESI. If ESI is fixed, postexatrasystolic strength increases as PESI is lengthened. If PESI is fixed, postextrasystolic strength increases as ESI is shortened. This effect includes mechanical restitution and pause dependence and it captures the important heart dynamics involving force generation and calcium cycling that occurs over multiple beats. Just like the electrical restitution protocol is a stringent test for electrical rate dependence, we believe that reproducing the PESP protocol is essential for an electromechanical cell model.
In this paper, we first present a systematical methodology to incorporate mechanical models into various existing EP models. We chose fourteen of the most recently developed EP models with multiple cellular compartments, advanced ionic currents and more realistic subcellular dynamics because they include the essential elements required to couple contraction models (e.g., SR dynamics and calcium buffers). We studied these models under the isometric condition by 1) evaluating how their calcium dynamics are affected by the inclusion of the contraction and 2) assessing how well these electromechanical models simulate contraction by comparing their PESP contractile response with experiments from Yue et al. and 3) analyzing the mechanism of PESP by finding the correlation between the calcium release from SR and the contractile strength. Finally we concluded with insights regarding the suitability of various models to reproduce electromechanical activities and the underlying mechanism of PESP, which may provide guidelines for future experiments.

Methods
In this section we introduced the contraction model that we implemented to EP models and the experiment we compared our simulations to, then we described our classification of EP models and corresponding strategy for contraction implementation, the numerical methods used to solve the models as well as the pacing protocols and data analysis. All terminology can be found in Table 1 The Contraction Model and Experiment We chose the Negroni_Lascano_1996 (NL96) [23]contraction model to represent the coupling of cross-bridge dynamics and intracellular Ca 2+ kinetics.
The muscle unit structure and cross-bridge dynamics are shown in Fig 1 in Negroni et al [23]. The muscle unit is composed of an inflexible thick filament (Myosin), a thin filament (actin) and an elastic paralleled element (titin). The cross-bridges are attached to the thick element on one end and can slide on the thin element on the other end. The total force is contributed by two parts: force generated from the elastic element (F p ) and force generated by the cross-bridges (F b ): where K, A and L 0 are constants; L is the sarcomere length; TCa Ã and T Ã are two states associated with troponin C sites on the thin filament (defined below); h is the elongation of the muscle unit. For "isometric" contraction where the sarcomere length is fixed, F p is constant, and the total force is determined only by the force generated from the cross-bridges. Cross-bridge grabbing and sliding on the thin element will generate force F b which is proportional to the elongation (h) and the total number of attached cross-bridges. The number of attached crossbridges is determined by the interaction with Ca 2+ and the relevant buffers, which is represented by a "four-state system" comprised of:1) sites on the thin element with free troponin C (T); 2) sites with Ca 2+ bound to troponin C (TCa); 3) sites with Ca 2+ bound to troponin C and attached cross-bridges (TCa Ã ); and 4) sites with troponin C not bound to Ca 2+ but attached cross-bridges (T Ã ). The system transitions among these four states via the binding and releasing of Ca 2+ and attaching and detaching of cross-bridges (See Fig 2 from Negroni et al [23]). The number of total cross-bridges is the sum of the two states that associated with cross-bridges (  Negroni et al [23] are: ½T ¼ ½Troponin total À ½T Ã À ½TCa À ½TCa Ã ð 5Þ These are reaction equations among the states of free troponin (T), troponin attached with cross-bridges (T Ã ), troponin bound with calcium (TCa), and troponin bound with calcium attached with cross-bridges (TCa Ã ); all with units of mM/ms. [TCa] eff is the effective [TCa] where cross-bridge attachment can occur and it's dependent on the sarcomere length (see Eq 10 in [23]). X = L − h, which is the inextensible length in the muscle unit. Y1 * Y4, Z1 * Z3 and Yd are rate constants. The NL96 model includes feedback from force on Ca 2+ ; therefore it is, in general, a "strongly coupled" contraction model, but under isometric condition it is "weakly coupled" since dX / dt = 0 [24].
There are two reasons why we chose the NL96 model. First, it shares common elements with the EP models; this greatly simplifies the implementation procedure. For example, the four states in the Ca 2+ kinetics are closely related to intracellular Ca 2+ concentration and Troponin concentration, both of which are common elements in advanced EP models. Second, even though NL96 model is a simplified model of the average effect of the interaction between myosin and actin, it can effectively reproduce basic physiological findings, such as time course of isometric force, intracellular Ca 2+ transient and force-length-Ca 2+ relation. Third, our method to incorporate NL96 can be easily applied to other contraction models as long as they use ODEs to describe the Ca 2+ troponin states, like [25][26] [27].
To evaluate the ability of these fourteen "coupled" electromechanical models to reproduce experimental findings, we compared the results of the simulations with the isovolumetric canine ventricular experiments from Yue et al. [22]. We chose Yue et al. because their paper includes not only a comprehensive experimental exploration of the relation between the contractile strength and pacing intervals but also a concise mathematical framework that allows quantitative comparisons of model predictions.
We applied the same PESP pacing protocol employed by Yue et al. [22], which was introduced in the Introduction section. By fixing ESI and varying PESI, Yue et al. [22] measured the rate of pressure change of the PES beat as a function of PESI; then by changing ESI to another value and repeatedly varying PESI, a family of curves can be constructed showing the effect of both ESI and PESI on pressure change (see Fig 2 from Yue et al [22]). In contrast to electrical restitution in which action potential duration is recorded for a single extrasystolic beat, the purpose of the PESP protocol is to study the contractile strength as a function of a single extrasystolic followed by a compensatory pause. The idea is to capture the full cycle Ca 2+ handling, specifically the release of calcium from the SR as a function of timing and the dynamics of the restoration of releasable SR Ca 2+ during (and following) action potential repolarization. The basic theory lies in the fact that Ca 2+ release from the SR is determined by two factors: the recovery of Ryanodine Receptor's (RyR's) which is a function of interbeat interval; and the amount of Ca 2+ in the SR. When ESI>PP, RyRs are more recovered from the last priming beat, resulting in a bigger Ca 2+ release from SR for the ES beat. Consequently, the Ca 2+ content in SR will be smaller compared to the priming period, and may lead to a smaller or larger Ca 2+ release (compared to PP) for the next PES beat depending on the PESI.
In Yue's experiments, they measured the maximum rate of pressure change in fourteen isolated perfused left canine ventricles under isovolumetric conditions. In the NL96 contraction model, force (normalized to muscle cross section area) is the variable to quantify the contractile strength; therefore to correlate with Yue's experiment results, we assumed a linear relationship between force generated from single cell and ventricular pressure [28], and thus we compare the normalized maximum rate of force change in our simulations with the normalized maximum rate of pressure change in Yue's paper.

Ventricular models and electrophysiological contraction implementation
We studied fourteen ventricular cell models, which include four species: guinea pig (n = 4) [ [42]. We classify the models into five categories based on the type of their calcium buffers and designed corresponding strategies to implement NL96. The complexity of the implementation of contraction into these five types of models is listed in ascending order. The flowchart in the S1 File illustrates steps of the implementation. We provide the essential equations, initial conditions and choice of parameters below. More details can be found in the S1 File. Table 2 lists general information about all models including classification, species, buffer information, number of variables, stimulus current amplitudes and durations; rate constants for Ca 2+ troponin buffers can be found in Table 3.
Type One: models with contraction. Two of the fourteen models have contraction in their original versions: Iribe_etal_2006 [30] and Matsuoka_etal_2003 [29].
Matsuoka_etal_2003 already has NL96 contraction in the original version. However, to investigate how much influence the mechanics has on this model, we remove the NL96 model We set K on and K off to the values of Y1 and Z1 and B max,troponin to the value in the original model. We call this new model Matuoska_etal_2003 without Contraction.
Iribe_etal_2006 includes the contraction model of Rice et al. (RWH99) in the original version [25]. For RWH99, CaTRPN is expressed by a single-state dynamic equation like Eq 6 but with a dynamic rate constant K off that depends on force, hence it is strongly coupled even for isometric contractions. Also the RWH99 model has six tropomyosin/cross-bridge states with rate constants that are functions of the Ca 2+ troponin buffer. To simulate this model without the contraction part, we eliminated the six tropomyosin/cross-bridge states and set K off for CaTRPN to a constant value by fixing the force in K off expression to half of its maximum value. We then compare simulations of the original model and the version without contraction as well as one with NL96 contraction model (implemented as described in the next "Type Two" section).
Type Two: models with full dynamic buffers. Two of the models incorporate ordinary differential equations (ODEs) for all the buffers: Shannon_etal_2004 [33] and Grandi_e-tal_2010 [37]. The equations for the original Ca 2+ troponin buffers in these models are the same as Eq 6. For models of Type Two we add the four-state NL96 model using Eqs 2-5 above; this is done by splitting the dynamical Ca 2+ troponin buffer in the original model into two states: with (TCa Ã ) or without (TCa) cross-bridges (Eqs 3 and 4): and adding the other two states T and T Ã (Eqs 2 and 5). There are multiple K on and K off in Eqs 2-4: Yd, Y1*Y4, Z1*Z3. They are kept the same as original NL96 paper [23] except for Y 1 and Z 1 , which are chosen to match the values of K on and K off in Eq 6 from the original model (see Table 3). This is to preserve the dynamics of CaTRPN as similar as possible to the original model. Type Three: models with dynamic Ca 2+ troponin buffers. Two of the models use ODEs for CaTRPN but instantaneous forms for all the other buffers: Mahajan_etal_2008 [34] and Iyer_etal_2004 [38].
For this type of models we follow the same step regarding the dynamic Ca 2+ troponin buffer in Type Two above and keep all the other buffers instantaneous. The only difference between Type Three and Type Two is that there are instantaneous buffer factors in Type Three models (see Eqs 9 and 10 below) but not in Type two.
Type Four: models with instantaneous Ca 2+ troponin buffers. Five of the models have all instantaneous Ca 2+ buffers: Hund_etal_2004 [35]; Faber_etal_2000 [31]; Livshitz_etal_2007 [32]; Ohara_etal_2011 [39] and Priebe_etal_1998 [40]. The buffering of Ca 2+ by Troponin is represented using the following equations: I Ca total represents the total intracellular Ca 2+ flux (see Page 18 in the supplement in reference [39] for detailed fomula); β is the instantaneous buffer factor; index j represents each type of intracellular Ca 2+ buffer; B max, j and K d,j are the total concentration and the affinity constant for buffer j. For this type of model we first change the instantaneous intracellular Ca 2+ troponin buffer into a dynamic buffer by eliminating the CaTRPN term from the instantaneous buffer factor β and by adding a dynamic CaTRPN flux into the [Ca 2+ ] i equation: The way to choose K on and K off here is not unique as long as indicates the values from the Shannon_etal_2004 model and s ¼ K d K d;s . We use the values from Shanno-n_etal_2004 as the standard values here because they have been widely used by to simulate dynamic CaTRPN buffers. After changing the instantaneous CaTRPN into a dynamical buffer, we follow the same procedure as for Type Three models.
For this type of models, we first add instantaneous Ca 2+ troponin buffers into the models. If the model has one general Ca 2+ buffer (referred as General) representing the average effect of all intracellular Ca 2+ buffers (e.g. TenTusscher_etal_2006 and Fink_etal_2008), we split the general buffer into two parts: CaTRPN and "Other". We keep K d,TRPN and K d,Other to be the same as K d,General in the original model so that the Ca 2+ affinity of the instantaneous buffer will be retained.
[Troponin] total was set to be 0.07mM, which is a standard value in most models.
[Other] total = [General] total − [Troponin] total so that the concentration of the total intracellular Ca 2+ buffer is the same as in the original model. If the model has other Ca 2+ buffers (such as CMDN) but no CaTRPN (e.g. Fox_etal_2002), we keep the other buffers unchanged and add a CaTRPN buffer. For the new CaTRPN buffer we also set [Troponin] total = 0.07mM and K d, CaTRPN = 0.6μM which are both common values in many models. After adding an instantaneous CaTRPN to the model, we follow the steps for Type Four to implement NL96.

Numerical integration
Except for one model (Iyer_etal_2004), all the code for the original single cell models were downloaded from www.cellml.org in CELLML format. Then they were translated into.mat files (MATLAB files) using a PYCML program. [43] [44] Simulations were run in MATLAB and integrated using the forward Euler integration method with time steps of dt = 0.001ms. The Iyer_etal_2004 model consists of 67 variables and requires a dt integration step as low as 1×10 −5 ms to converge [38] using forward Euler, thus becoming impractically slow to simulate in MATLAB. Therefore as in [45] the Iyer_etal_2004 model was written in FORTRAN using a semi-implicit integration method that allows a much larger (while still convergent) integration time step of dt = 0.005ms.

Pacing protocol and initial conditions
Our pacing protocol is composed of three steps and similar to the pacing protocol from Yue et al [22] except that we chose to use a priming period of 500ms. The stimulus currents and durations for the different models were selected to ensure excitation and the values are provided in Table 2. All simulations were run for isometric contraction of a single cell, where the half sarcomere length is fixed to be 1.05μm.
Step one: Quiescent. In this step, EP models without contractions were run with no stimulation current until quiescent steady states were reached (i.e., no state variable changed by more than 0.01%). Initial values for state variables in each model were directly loaded from CellML code. The time for each model to reach the quiescent state is listed in Table 3, and it ranges from 10min of real time (as in the Mahajan_etal_2008 model) to up to 100min (as in the Hund_etal_2004 model). Note that these are real times and not run times. Steady state quiescent values for voltage, intracellular and JSR (or SR) Ca 2+ concentrations vary considerably among models as shown in Table 4.
Step Two: Priming Cycle. During the priming pacing cycle, all models (with and without contractions) were paced with priming period T = 500ms starting from the quiescent states at the end of Step one, until a new "priming" steady state was achieved. The beat number required to reach this priming steady state for each model is listed in Table 2. For models incorporating NL96 contraction, state variables for the four states of Ca 2+ troponin buffer were calculated by solving Eqs 2-5 for steady state using the quiescent value of [Ca 2+ ] i . The total Ca 2+ troponin buffer ([TCa]+[TCa Ã ]) calculated using this method was verified to be similar to the quiescent value in the original EP model if the original model had a Ca 2+ troponin buffer. We used this method to ensure that the initial conditions for the original models and the corresponding contraction models were as close as possible.
Step Three: Postextrasystolic pacing protocol. For all models, an extrasystolic (ES) and a postextrasystolic (PES) beat were delivered in sequence after the last priming beat as described in the Contraction Model and Experiment section. We held ESI at a constant value and varied PESI to record the transients of transmembrane potential, [Ca 2+ ] i , CaTRPN and generated force; we then changed ESI to a different value and repeated the process using the priming steady state as initial conditions. The shortest PESI is the refractory period of the ES beat so it is just long enough that the PES beat is separated from the ES beat. The longest PESI (ESI) value was chosen as 1200ms, 1500ms or 2000ms, depending on the time it took for the Postextrasystolic MRC pes (MRC es ) to converge to a plateau level, as some models took longer than others to reach it.

Data Analysis
All curves were fitted using MATLAB's built-in function "fit" whose default curve fit algorithm is the Trust-Region method. The value of r 2 was given for each fit and 95% confidence bounds were provided for each fitted coefficient. Our analysis is divided into two parts: Ca 2+ dynamics and Contraction. To analyze Ca 2+ dynamics we compare the Ca 2+ transients for the priming beats and during postextrasystolic potentiation between the original models and models with contraction. To analyze Contraction we generate four characteristic contraction curves and fit two of them into monoexponential functions.
Ca  We quantified the differences of the parameters between original models and models with contractions by the relative difference, defined as Postextrasystolic Potentiation in [Ca 2+ ] i : For each model, we plotted the intracellular Ca 2+ transient for the last 3 priming beats, along with multiple ES/PES beat pairs overlaid on the same graph for each model (without contraction on the left and with contraction of the right).
Contraction. Following the data fitting and analysis from Yue paper [22] with modification, we generate four characteristic curves for each model after the implementation of contrac- where F denotes the force generated from NL96 model; PES denotes that it is the postextrasystolic beat; SS denotes the steady state during the priming period; CR max,pes is the plateau amplitude, termed as Maximum Postextrasystolic Contractile Response; C 0 is the minimum value of dF = dt max ðPESÞ dF = dt max ðSSÞ , in Yue et al. [22] C 0 = 0; t o,pes is the minimum-value axis intercept where MRC pes intercepts the minimum value line dF = dt max ðPESÞ dF = dt max ðSSÞ ¼ C 0 ; and T mrc,pes is the time constant for MRC pes . . Fig 1(A) shows one set of MRC pes generated using equations and parameters in Yue et al. [22].
Extrasystolic (PES) Mechanical Restitution Curves (MRC es ) presents the maximum rates of change of force of extrasystolic beats normalized to the last priming beat (dF / dt max (ES) / dF / dt max (SS)) as a function of ESI. Since it is a subset of MRC pes with ESI equal to the priming cycle length, we do not present them separately.
Postextrasystolic Potentiation Curve (PESPC) is the maximum postextrasystolic contractile response (CR max,pes ) plotted vs ESI. CR max,pes is the plateau amplitude of the above MRC pes for a fixed ESI and this amplitude decreases to a plateau value as ESI increases. Fig 1(B) shows the PESPC and the simultaneously determined MRC es , both generated from the equations in Yue et al. [22]. The equation for PESPC fitting is: where B is the plateau level; A is the amplitude; t o,es is the ESI-axis intercept for the simultaneously determined Extrasystolic Mechanical Restitution Curve; T pespc is the time constant for PESPC.

Results
Our Results section is divided into three parts: Calcium Results, Contraction Results, and the Underlying Mechanism for PESP. In the Calcium Results section we first compared the Priming Ca 2+ transients between original EP models and the corresponding models with the NL96 contraction model to investigate the influence of including contraction on EP models. Then we studied the Postextrasystolic Potentiation behavior in Ca 2+ before and after the implementation of contraction. In the Contraction Results section we compared the four characteristic contractile curves with the Yue et al. experimental results. In the Underlying Mechanism for PESP section we analyzed the mechanism of PESP by finding the correlation between the calcium release from SR and the contractile strength.

Ca 2+ Results
Priming [Ca 2+ ] i properties are undisturbed after the implementation of contraction for models with dynamic buffers. We present the priming [Ca 2+ ] i data in  Table 5  (Table 6) and the relative differences between before and after contraction are shown in Table 7, where double dagger symbols ( ‡) indicate the differences larger than 0.5(50%); single dagger symbols ( †) indicate differences larger than 0.1 (10%) but less than 0.5 (50%).
[Ca 2+ ] i differs substantially among the original EP models (without contraction) as shown in Table 5  The most significant change tended to be an increase in the maximum [Ca 2+ ] i , always accompanied by a decrease in the time of this peak (t peak ). As shown in Table 7, all six models with dynamic Ca 2+ troponin buffers have relative differences of less than 10% in systolic [Ca 2+ ] i and five out them have relative differences less than 10% in t peak . On the other hand, among the eight models with instantaneous or none Ca 2+ troponin buffers, six models show differences larger than 50% and none has a difference less than 10% in systolic [Ca 2+ ] i ; in t peak , one of the eight models (Fink_etal_2008) has a relative difference less than 10% and three of them have differences of more than 50%. In addition, except for one model (Fox_etal_2002), all of them have positive relative differences in systolic [Ca 2+ ] i, indicating a systolic [Ca 2+ ] i increase after the implementation of the NL96 contraction; and all eight models have negative Δt peak , meaning the time of the peak is decreased by the implementation of contraction. The increase in maximum [Ca 2+ ] i and the decrease in t peak in models with instantaneous Ca 2+ troponin buffers after the implementation of NL96 (Fig 2, S2 File) is due to an increased and fast rate change of the intracellular calcium after its dynamics has been modified from Eqs 9 and 10 (before the implementation) to Eqs 11 and 12 (after the implementation) to account for the contraction . Fig 3(A) and 3(B) shows, as an example, the [Ca 2+ ] i and d[Ca 2+ ] i / dt respectively for the Ohara_etal_2011 model before (solid lines) and after (dashed line) the NL96 implementation. It can be seen that after the implementation of NL96: (i) the instantaneous buffer factor β increases due to the elimination of the instantaneous Ca 2+ troponin buffer from the denominator (see Eqs 10 and 12); and (ii) the flux I Ca total is reduced due to the addition of the negative dynamic Ca 2+ troponin flux (see Eqs 9 and 11). The flux I Ca total decreases by less than 50% of the original value (Fig 3(D)) but β increases by more than 200% (Fig 3(C)) resulting in an sharp increase in [Ca 2+ ] i amplitude. And the time for β to reach its peak is clearly shortened, resulting in the decrease of [Ca 2+ ] i peak (t peak ).
Results for the priming [Ca 2+ ] i duration (t 1/2 ) are similar to those of systolic [Ca 2+ ] i (see Table 7). All six models that have dynamic Ca 2+ troponin buffers show little change in t 1/2 (less than 10%). Among the eight models with instantaneous or no Ca 2+ troponin buffers, six have differences more than 50% and none has a difference less than 10%. Seven out of eight have negative The diastolic [Ca 2+ ] i and the Ca 2+ charge (Q) during one beat are the two parameters that do not change much after implementing NL96 contraction in all of the models. Thirteen out of fourteen models have diastolic [Ca 2+ ] i changes of less than 10%, with only one model having a difference larger than 50% (Fox_etal_2002). Q for eleven out of fourteen models changes by less than 10% and only one model (Fox_etal_2002) has a change of more than 50% (see Table 7).This is because, although the systolic [Ca 2+ ] i is greatly increased for models with instantaneous or none Ca 2+ troponin buffers, its increase is localized to a relatively narrow spike with an area that is negligible comparing with that of the rest of the [Ca 2+ ] i transient.   (c) (d) Iribe_etal_2006 without/with contraction (NL96); (e) (f) Shannon_etal_2004 without/with contraction. The first three beats are priming beats; from beat a to beat b are ES beats with different ESIs; from a' to b' are PES beats with a fixed PESI (fully restituted PESI, defined in Step Three in the Pacing Protocol of increase in ES beats nor decrease in PES beats. Note that these types of models, as discussed before, show a change in the shape of the [Ca 2+ ] i when contraction is included. However it doesn't change the postextrasystolic potentiation behavior much, namely, if the original model has a clear PESP trend, after the implementation, the new model will still maintain that behavior and vice versa.

Contraction Results
To investigate how well our EP models with updated contraction can represent postextrasystolic dynamics, we calculated the four characteristic curves, described in the Methods, for all fourteen models with the NL96 contraction and compared them with the experimental data from Yue et al. [22]. Note that we did not show any CaTRPN transient because they are linearly related to the contraction in isometric condition [23].
Postextrasystolic Mechanical Restitution Curve (MRC pes ). Representative examples of MRC pes (dF / dt max (PES) / dF / dt max (SS)) for six models are shown in Fig 6. For all the models, MRC pes are given in S4 File.
For all models MRC pes is well described (r 2 >0.99; averaged over all ESI values) by a monotonically increasing curve with respect to PESI (Eq 15) as shown in Table 8, after the elimination of some data points (see the example of Mahajan_eral_2008 below). However not all models were consistent with the experimental findings of Yue et al. [22] whose C 0 = 0. Among the six models with dynamic Ca 2+ troponin buffers (Type One, Two, Three), four satisfy this condition, while among the eight models with instantaneous or none Ca 2+ troponin buffers (Type Four and Five), only one satisfies C 0 = 0 (see Table 8). Here we provide a brief explanation why C 0 does not go to zero for some models in our simulations. We use Mahajan_e-tal_2008 with NL96 as an example to demonstrate. In Fig 7(A), 7(B) and 7(C) we plot the Action Potential (AP), [Ca 2+ ] i transient and normalized dF / dt transient of Mahajan_e-tal_2008 with NL96 for ESI = 500ms. The first beat is the steady priming beat, the second beat is the ES beat with ESI = 500ms and from a' to d' are PES beats with PESI from 200 to 650ms. We can observe that the dF / dt max increases with PESI from beat c' to d' but there is no welldefined trend before c'. Beat a' and b' are not clearly separated from the ES beat and a' is even almost fused into the ES beat. However, even though beat a' is very close to the ES beat, dF / dt max still does not go to zero. We plot normalized dF / dt max (normalize to steady priming beat) for ESI = 500ms and PESI ranges from 200ms to 1200ms in Fig 7(D). We can see that dF / dt max stops decreasing when PESI is below 300ms and the minimum value is about 0.4. In order to fit this into Eq 15 we exclude data with PESI lower than 300ms and set C 0 = 0.4. That is the reason for Mahajan_etal_2008 and other eight models, that C 0 does not go to zero and also an example of how we choose the shortest ESI and PESI values.
The experimental results of Yue et al. [22] showed that the fully restituted plateau value of MRC pes (CR max,pes ) decreased as ESI increased. From Fig 6 and S4 File we can see that all models with dynamic Ca 2+ troponin buffers show unique plateau values for different ESIs while models with instantaneous or no Ca 2+ troponin buffers converge to the same plateau values for different ESIs; this property will be further discussed in the PESPC subsection.
The time constant for MRC pes (T mrc,pes ) did not vary much as a function of ESI in the experiment of Yue et al. [22]. In addition, a leftward shift of the MRC pes was observed as ESI decreased, presented as an increase of t o,pes as ESI increased. The fourteen models behave the Methods section); a and a' are corresponding ES and PES beats; so are b and b'. Note that each pair of figures (without/with contraction) has the same axis range. (e) (f) TenTusscher_etal_2004 without/with contraction. As in the previous figure, the first three differently regarding these two parameters. We will discuss these in detail in the t o,pes and T mrc,pes subsections.
Postextrasystolic Potentiation Curve (PESPC) . Fig 8 shows representative postextrasystolic potentiation curves (PESPCs) for six models. Figures for all models can be found in S5 File. To provide a better picture of the amplitude and time constant of the PESPC we simultaneously plot the PESPC and the Extrasystolic Mechanical Restitution Curve (MRC es ) as in Yue et al. [22]. Note that for those models with C 0 6 ¼ 0 in the MRC pes , we have shifted the curves upward by C 0 when plotting the PESPC (star Ã ) on the same graph with the MRC es (open circle o) for comparison. We also use 95% confidence bounds as error bars for each data point on PESPC, which is the fitting confidence bound for CR max,pes of each ESI value. Note that confidence bounds are relatively small compared to the scale of the graph.
A list of all parameters for the curve fitting of Eq 16 are listed in Table 8 for all the models along with the experimental data of Yue et al. [22] to allow comparison between our simulation results and those from the experiments. The values without any symbol match the experiment results well. They are either parameters with 95% confidence bounds overlapping the mean plus/minus standard deviation (mean±SD) range from Yue et al. [22] or they have r 2 values that are higher than 0.99. The double and the triple star ( ÃÃ and ÃÃÃ ) means the 95% confidence bound overlaps Yue et al. [22] data's mean±2 Ã SD and the mean±3 Ã SD respectively. The X mark means the confidence bound overlaps the mean±SD but the confidence bound is too large (more than 50% of the center value) so that we do not consider this as a real overlapping. The double X mark (XX) means the confidence bound does not overlap the mean±n Ã SD from Yue experiment with n 3.The single-dagger ( †) means the r 2 is less than 99% but more than 90% and the double-dagger ( ‡) means it's less than 90%.
The experiments of Yue et al.showed that PESPC decreased monoexponentially to a plateau level as ESI increases (Eq 16 and Fig 1(B)), meaning that when PESI is fixed and long enough, the maximum force changing rates of the PES beats will decrease as ESI increases. Models with dynamic Ca 2+ troponin buffers (Type One, Two and Three) have PESPCs well described by Eq 16; all six models have r 2 values larger than 0.99. On the other hand, models with instantaneous or no Ca 2+ troponin buffers (Type Four and Five) are not well fit by Eq 16. Four of the eight models have r 2 values lower than 0.9 with one of them as low as 0.0048 (Priebe_etal_1998). However, the poor quality of the fitting is not evident on Fig 8 or S5 File because the amplitude of the PESPC is much smaller than the scale of the figure.
The plateau levels (B) of CR max,pes for the fourteen contraction models are more similar to the experiment results compared to the other parameters. In Yue et al. [22], the mean value (± SD) for B is 1.05±0.13. In our simulations, two out of the six models with dynamic Ca 2+ troponin buffers (Type One, Two and Three) have 95% confidence bounds overlap the mean±SD range of Yue et al. [22] and four have confidence bounds overlap the mean±2 Ã SD range. On the other hand, only two out of eight models with instantaneous or no Ca 2+ troponin buffers (Type Four and Five) have 95% confidence bounds overlap the mean±2 Ã SD range.
Parameter A has the least matching degree with Yue et al. [22] where the mean value (± SD) is 1.68±0.32. Only one out of the six models with dynamic Ca 2+ troponin buffers (Type One, Two and Three) have 95% confidence bounds overlap the mean±2 Ã SD range from Yue paper and none of the models with instantaneous or none Ca 2+ troponin buffers have confidence bounds overlapping this range (see Table 8). The values of A from our contraction models are significantly smaller than Yue paper. This feature is well depicted in Fig 8 and S5 File in which beats are priming beats; from beat a to beat b are ES beats with different ESIs; from a' to b' are PES beats with a fixed PESI And each pair of figures (without/ with contraction) has the same axis range. the PESPCs of the models with instantaneous buffers have slopes of approximate zero while in the Yue et al. paper there exhibited a pronounced monoexponentially decay (Fig 1(B)). The insufficient change in fully-restituted contractile strength of postextrosystoles with respect to varied ESIs indicates a clear limitation of these models as we will further comment in the discussions section.
t o,es is obtained by calculating where the Extrasystolic Mechanical Restitution Curve (MRC es ) intercepts the line dF = dt max ðESÞ dF = dt max ðSSÞ ¼ C 0 . The mean value (±SD) in the Yue et al. [22] is 284±32ms. Four out of six models with dynamic Ca 2+ troponin buffers (Type One, Two and Three) and five out of eight models with instantaneous buffers (Type Four and Five) are in the mean±2 Ã SD range. For this parameter, it seems that the dynamic buffers do not show much advantage over instantaneous buffers. However, in the Yue et al. paper, this interception is where MRC es intercepts with the line dP = dt max ðESÞ dP = dt max ðSSÞ ¼ 0, so are four out of the six models with dynamic Ca 2+ troponin buffers. But for most models with instantaneous buffers, this parameter is where MRC es intercepts with the line dF = dt max ðESÞ dF = dt max ðSSÞ ¼ C 0 with C 0 6 ¼ 0. Therefore models with dynamic buffers still match this parameter better to the experimental data of Yue et al. [22].
The mean value (± SD) of the time constant for PESPC (T pespc ) is 176±18ms in Yue et al. [22]. Only two out of the six models with dynamic Ca 2+ troponin buffers (Type One, Two, Three) and none of the models with instantaneous or none buffers (Type Four and Five) have 95% confidence bounds that overlap the experimental mean±2 Ã SD range. Ten other models have time constants much larger(50% more) than Yue et al. [22]. This parameter is also the  Table 9.
In the Yue et al. paper [22], t o,pes increased as ESI increased. They explained this trend by suggesting that the refractory period is an increasing function of ESI [22]. Five out of six models with dynamic buffers (Type One, Two and Three) show the same trend as Yue et al.

Underlying Mechanism for PESP
There is a clear resemblance between the intracellular Ca 2+ dynamic and force dynamics. And since the Ca 2+ released from SR (J rel ) is the main source for intracellular Ca 2+ , we propose that J rel dynamics, which depends on RyR, Ca 2þ SR , among others, is the primary determinant of PESP dynamics.
To verify this hypothesis, we analyzed, for each model, the correlation between the dependence of the potentiation in Ca 2+ released from SR (J rel ) and the contractile strength on pacing intervals. In this study we have shown two sets of curves demonstrating the potentiation strength: the Postextrasystolic Mechanical Restitution Curve (MRC pes ) and Postextrasystolic Potentiation Curve (PESPC). The MRC pes demonstrates how postextrasystolic contractile strength changes with respect to various PESI for a fixed ESI. From Fig 6 and S4 File, it can be seen the consistency of the increasing trend of the normalized force changing rate with respect to PESI for all models studied. However, the PESPC, which shows how postextrasystolic contractile strength changes with respect to various ESI values for a fixed and long PESI, differs among models: models with dynamic buffers showed significant variation in contractile strength with respect to ESI while models with instantaneous buffers showed little variation (Fig 8 and S5 File). The amplitude (parameter A) in PESPC captures this difference. Since we are interested in what causes the difference in models, we studied the correlation between the Ca 2+ released from SR (J rel ) and the contractile strength of the postextrasystoles with various ESI and fixed PESI. We designed an analogous parameter A J for the J rel amplitude. Similar to our dF / dt max fit, for a fixed ESI, we first fit the normalized J rel_max (to the prime beats) of the postextrasystolic beats to a mono-exponentially increasing curve with respect to PESI (analogous to MRC pes ). For most models, J rel_max was well described by a mono-exponential increasing curve, similar to the contractile strength. Then we fit the plateau values of these curves to a mono-exponentially decreasing curve with respect to ESI (analogous to PESPC); we define the amplitude of this curve as A J . Models with larger A J suggest that their fully-restituted J rel_max vary with respect to ESI. Therefore the correlation between A and A J can be used as a gauge for the correlation between the dependence of J rel and the contractile strength on pacing intervals. We performed a linear regression for A J and A for thirteen (out of the fourteen) electromechanical models as shown in Fig 11. We found a strong positive correlation (r 2 = 0.848) between the two parameters. In addition, data from models with instantaneous buffers clustered around the origin while data from models with dynamic buffers were distributed along the regression curve. This means for models with instantaneous buffers, the insufficient variation in the fully-restituted postextrasystolic contractile strength responding to various ESI is a direct result of the J rel dynamics. Note that the model excluded from the linear regression was TenTusscher_etal_2006 (upward triangular mark in Fig 11) because the curve fitting for parameter A in this model is not reliable due to its large 95% confidence bound (see Table 8).

Discussion
Simulations of electrical therapies for HF (e.g., CRT and CCM) require models that accurately reproduce the excitation-contraction coupling that occurs in the cardiac myocyte. While action potential duration restitution curves tend to be used to help validate EP models, there has been no systematic method to validate cellular contraction models. We suggest that the PESP protocol can be used in the validation process of mathematical models of cellular electromechanics and we present a methodology to include mechanics into cellular EP models of various types and compare results among the models with and without contraction.
The results from this study suggest four reasons for which electrophysiology models with dynamic Ca 2+ troponin buffers are much better candidates to implement contraction than EP models with instantaneous or none Ca 2+ troponin buffers.
First, the complexity in implementation is greatly reduced since the original models are already equipped with dynamic equations for Ca 2+ troponin buffers. Type Two and Type three models only take one step in the implementation flowchart while Type Four and Type Five require one or two more steps to add in instantaneous Ca 2+ troponin buffers and change the instantaneous ones into a dynamic form (see S1 File).
Second, the inclusion of contraction does not disturb properties in the original models. The [Ca 2+ ] i shape shows minimum change for models with dynamic buffers while it is largely distorted for models with instantaneous or none Ca 2+ troponin buffers (Fig 2 and S2 File). This is closely related to the first benefit as it requires adding extra terms into Type Four and Type Five models to incorporate contraction and therefore Ca 2+ shapes are more likely to be unphysiologically deformed.
Third, models with dynamic Ca 2+ troponin buffers show clear postextrasystolic potentiation in [Ca 2+ ] i with or without contraction implemented (see Figs 4 and 5 and S3 File). Therefore these models show desired variation in contractile strength of extrasystoles and postextrasystoles because of the close relationship between Ca 2+ and contraction.
Fourth, contraction properties of models with dynamic Ca 2+ troponin buffers fit experimental results better [22] (see Figs 6,8,9 and 10). This has been shown in the four characteristic curves (Postextrasystolic Mechanical Restitution Curve (MRC pes ); Postextrasystolic Potentiation Curve (PESPC); Minimum-value Axis Intercept Curve (t o,pes ) and Time Constant Curve (T mrc,pes )) and their various parameters. To unveil the mechanism why model with dynamic buffers reproduce PESP better than models with instantaneous buffers, we studied the correlation between contraction and calcium release from SR. We found a strong positive correlation (r 2 = 0.848). This means for models with instantaneous buffers, the insufficient variation in contractile strength responding to various pacing cycle length is a direct result of the J rel dynamics. The interpretation for this result is twofold. First, Ca 2+ release from SR is the main source for intracellular Ca 2+ therefore J rel can be the direct reason of contraction behaviors. Second, time constants play a crucial role in cell contraction. The time constants of Ca 2+ binding with troponin C in Eq 6 in models with dynamic buffers are in the order of 10ms, which is the same order of the time scales of the upstroke of Ca 2+ , CaTRPN and tension development. The RyR release kinetics is highly timing sensitive and strongly coupled to the CaTRPN binding process. Therefore models with dynamic buffers, which can describe the timing of interaction among different cell compartments more accurately, reproduce correct contraction behaviors. This has biological significance since most electrophysiological models incorporate instantaneous buffering due to their simplicity and lower computational cost.
To assess the robustness of the contraction model we use, we conducted more careful verifications in the supplementary material S8 and S9 Files.
In S8 File, we compared NL96 with two other contraction models: i) NL96 and RWH99 implemented in Iribe_etal_2006; ii) dynamic (original) NL96 and instantaneous NL96 implemented in Ohara_etal_2011. In the first set of comparison, we conclude that under the isometric condition there are no significant differences between NL96 and RWH99 under the scope of this study in the sense that each one matches different parameters with experiments better than the other. In the second set of comparison we find that although including NL96 in models with instantaneous Ca 2+ troponin buffers distort their Ca 2+ shape, the poor contraction behaviors of these models are not directly related to that. Two simulation results support this conclusion. First, their original models do not show clear postextrasystolic potentiation behavior in [Ca 2+ ] i (see Fig 5), which is directly related to the lack of variation in contractile strength and second, Ca 2+ shape is retained after the implementation of the instantaneous NL96 yet the contraction behavior of Ohara_etal_2011 is still not optimal.
Finally in S9 File, we conducted verifications of physiological parameters such as different priming period, low temperature and cooperativity of the contraction model. By comparing simulation results before and after modification of these physiological parameters, we showed that the choice of parameters did not change quantitatively the PESP behavior.

Limitations and future work
Yue's experiments are based on dog ventricles in iso-volume conditions, but our simulations are on isometric single cells. However similar experiment results have been observed in smaller tissue sizes. Wier et al. [28] measured tension in perfused ferret papillary muscles 0.67 ±0.05mm in external diameter and reported similar monoexponential curves as Yue et al.. We do consider the size to be a potential limitation of our work. Qualitatively we get similar characteristic curves in our simulations as in the experiment but quantitatively none of our models show as much variation in contraction strength with respect to changing stimulus rate (see Figs 6 and 8) as in Yue's experiment. One possible reason might be that the coupling and isotropy in tissues reinforce the variation in contractile strength, which suggests tissue simulations should be one direction for future work.
Another limitation on using some of the cell models is the coupling of cross species models and the comparison to experimental data from another species. The argument is that, first, the parameters in the experiment and simulations that we performed quantitatively cross species comparisons on are all normalized to their own priming beats (Figs 6 and 8); the normalization should reduce the species sensitive factors. Second, the quantities that are not normalized are either compared before and after the contraction implementation within each own individual model (Figs 2, 4 and 5) or compared qualitatively with the experiment (trends with respect to ESI in Figs 9 and 10). Nevertheless we admit this is a limitation in our work due to the lack of experiment data for all species.
In our study, we chose a hybrid approach to implement contraction where we only changed CaTRPN into a dynamic buffer while kept all the other buffers instantaneous. The reason for that is because we wanted to modify the original model as little as possible. Since only CaTRPN is presented in the contraction model, we changed only that into a dynamic buffer. We however did conduct simulations on O'Hara et al model to investigate the effect of dynamic changes in other buffers, specifically we changed both of the instantaneous intracellular calcium buffers in the model (troponin and calmodulin) into dynamic buffers using the same strategy described in the Method section. We found that doing these changes not only affected the calcium as before but in addition the voltage AP, both by large amount (not shown).