A Mathematical Model of the Mouse Ventricular Myocyte Contraction

Mathematical models of cardiac function at the cellular level include three major components, such as electrical activity, Ca2+ dynamics, and cellular shortening. We developed a model for mouse ventricular myocyte contraction which is based on our previously published comprehensive models of action potential and Ca2+ handling mechanisms. The model was verified with extensive experimental data on mouse myocyte contraction at room temperature. In the model, we implemented variable sarcomere length and indirect modulation of the tropomyosin transition rates by Ca2+ and troponin. The resulting model described well steady-state force-calcium relationships, dependence of the contraction force on the sarcomere length, time course of the contraction force and myocyte shortening, frequency dependence of the contraction force and cellular contraction, and experimentally measured derivatives of the myocyte length variation. We emphasized the importance of the inclusion of variable sarcomere length into a model for ventricular myocyte contraction. Differences in contraction force and cell shortening for epicardial and endocardial ventricular myocytes were investigated. Model applicability for the experimental studies and model limitations were discussed.


Introduction
Cardiac cell functions include the interaction of several major subsystems, including those responsible for the generation of electrical activity, Ca 2+ dynamics, and cardiac contraction. Experimental data from diseased hearts or obtained at fast pacing rates show that the changes in one of the subsystems can lead to abnormal behavior in others. For example, dysfunction of the Ltype Ca 2+ channel, as in Timothy syndrome when the channel's inactivation is significantly reduced, affects Ca 2+ handling in cardiac cells [1,2] resulting in cardiac arrhythmias. Heterogeneities in cellular electrical activities in the heart, dysfunction of K + channels, or acidosis can also produce pro-arrhythmic behavior, such as action potential propagation block, re-entry, Ca 2+ alternans, and irregular contractions [3,4]. In particular, instability of Ca 2+ dynamics (alternans) can lead to the action potential alternans [5] and alternans in mechanical contraction [6]. Therefore, understanding interactions of the major cardiac cell subsystems and mechanisms of their pro-arrhythmic activity is of great importance.
Mathematical modeling of electrical activity, Ca 2+ dynamics, and cardiac contraction is a supplementary tool for experimentalists in order to understand mechanisms of pro-arrhythmic activity in the heart. There are several models for cardiac myocyte contraction that have been developed for different species. Such models were developed for guinea pig [7,8], rabbit [9], canine [10], and mouse [11] ventricular myocytes. The models include experimentally-verified sets of ionic currents, Ca 2+ dynamics, and contractile parameters for cardiac cells of the particular species.
Myocyte contraction is a complex process which involves activation of ionic currents, including L-type Ca 2+ current, through which Ca 2+ enters the cell and causes Ca 2+ release from the intracellular Ca 2+ store, the sarcoplasmic reticulum [12]. High intracellular Ca 2+ concentration leads to an increase in Ca 2+ bound by intracellular proteins (troponin, calmodulin) and changes the myofilament configuration, resulting in force development. Force generation involves conformational changes in thick (myosin) and thin (actin, tropomyosin, and troponin) filaments (Fig. 1A) resulting in an increase in their overlap. Myosin represents a polypeptide chain with globular heads, which constitute crossbridges that interact with thin filaments. Thin filaments are composed of long tropomyosin polypeptide chains, on which globular actin molecules aggregate in double-stranded helix with crossbridge binding sites. In a non-active configuration, troponin blocks crossbridge binding sites. Upon Ca 2+ binding to troponin, troponin-tropomyosin complex exposes crossbridge binding sites which interact with myosin globular heads, thereby creating weak bonds. ATP molecules bound to actin release a phosphate group and transform weak bonds into strong bonds. This transformation results in a change of crossbridge conformation to a bent position and forces thick filaments to slide relative to thin filaments.
Because of the complexity of the contraction mechanism, most mathematical models use a significantly simplified description of this process [13]. They explore the Huxley two-state crossbridge model [14], extend it to a larger number of crossbridge states, and include direct and indirect interaction with troponin and variable sarcomere lengths [13,15]. Such simplified description, for example, does not involve energy metabolism and interaction with mitochondria. The crossbridge models are further incorporated into cellular models, which include electrical activity, comprehensive Ca 2+ dynamics [7,9,11], and energy metabolism [16,17].
A mathematical model which describes mouse ventricular cell contraction for body temperature (310uK, or +37uC) was developed and published recently by Land et al. [11]. The model adjusted the cellular contraction model of Rice et al. [9] for mice and is primarily devoted to simulation of the whole heart contraction. While some experimental data is recapitulated by the model, such as time course of the myocyte contraction, the deviation of the simulated from the experimental data and model limitations are also discussed by the authors [11]. For example, absolute values of the simulated cellular contraction force are considerably larger than the measured contraction force for both physiological and room temperatures [18]. In addition, Land et al. [11] did not study contraction force-frequency relationships at the cellular level.
In this paper, we developed a new electromechanical model for mouse ventricular myocyte contraction at room temperature (298uK, or +25uC). We employed our previously published models for action potential and Ca 2+ dynamics in mouse ventricular myocytes [19,20,21,22], which were also developed for room temperature (298uK, or +25uC), and incorporated a myocyte contraction model from Rice et al. [15]. These models were successfully employed for simulations of proarrhythmic activities in mouse cardiac cells and tissues [21,22]. In addition, in the Rice et al. [15] model, we implemented sarcomere length variation during twitch. We also explored the effects of heterogeneity of the electrical activity and Ca 2+ dynamics in epicardial and endocardial cells on the contraction force generation and cell shortening. The resulting model was adjusted to fit experimental data on mouse ventricular cell contraction. Our model successfully reproduced steady-state force-calcium relationships for different sarcomere lengths; time courses of the Ca 2+ transients, developed force, and cellular shortening; peak force-frequency and cell shorteningfrequency relationships; and time-to-peak force and time-to-50% force relaxation. We also investigated and emphasized the importance of using variable sarcomere lengths in models of myocyte contraction. In the simulations, we compared both the  [19]. Transmembrane currents are the fast Na + current (I Na ), the L-type Ca 2+ current (I CaL ), the sarcolemmal Ca 2+ pump (I p(Ca) ), the Na + /Ca 2+ exchanger (I NaCa ), the rapidly recovering transient outward K + current (I Kto,f ), the slowly recovering transient outward K + current (I Kto,s ), the rapid delayed rectifier K + current (I Kr ), the ultrarapidly activating delayed rectifier K + current (I Kur ), the noninactivating steadystate voltage activated K + current (I Kss ), the time-independent K + current (I K1 ), the slow delayed rectifier K + current (I Ks ), the Na + /K + pump (I NaK ), the Ca 2+ -activated chloride current (I Cl,Ca ), the Ca 2+ and Na + background currents (I Cab and I Nab ). I stim is the external stimulation current. The Ca 2+ fluxes within the cell are uptake of Ca 2+ from the cytosol to the network sarcoplasmic reticulum (SR) (J up ), Ca 2+ release from the junctional SR (J rel ), Ca 2+ flux from the network SR (NSR) to junctional SR (JSR) (J tr ), Ca 2+ leak from the SR to the cytosol (J leak ), Ca 2+ flux from the subspace volume to the bulk myoplasm (J xfer ), Ca 2+ flux to troponin (J trpn  absolute value of the contraction force and cellular shortening, and their normalized dependences to fit existing experimental data.

Methods
A mathematical model for mouse ventricular myocyte contraction is a natural extension of our previously published models for action potential and Ca 2+ dynamics in mouse ventricular myocytes [19], with model improvements from [20,21,22] (Fig. 1A), developed for room temperature (298uK, or +25uC). In this paper, we explore mouse ventricular myocyte models from the epicardial and endocardial regions of the heart [22]. Endocardial cells have more prolonged action potentials and larger intracellular [Ca 2+ ] i transients compared to epicardial cells [22]. We incorporated the Rice et al. [15] contraction model 4 in our model of electrical activity and Ca 2+ handling [19,20,21,22] (See Appendix S1) and adjusted model parameters to fit experimental data on myocyte contraction obtained for room temperatures.
The Rice et al. [15] model links Ca 2+ dynamics and myocyte contraction (Fig. 1B). The model contains two nonpermissive tropomyosin states (N0 and N1) and four permissive tropomyosin states (P0, P1, P2, and P3). N0, N1, P0, P1, P2, and P3 are functions of time that describe probabilities of finding the model in that particular state. N0 is the rest state of the model, with no strongly bound crossbridges. When Ca 2+ binds to the tropomyosin, it changes its conformation to a permissive state without strongly bound crossbridges (P0), which allows for strong binding of one (P1), two (P2), or three (P3) crossbridges. The model also includes one nonpermissive state with one strongly bound crossbridge even without a bound Ca 2+ ion. All transition rates in the model are Ca 2+ -independent, except for k NP , which depends on the concentration of troponin with Ca 2+ bound to a low-affinity binding site. Detailed analysis of several contraction models and the plausibility of different cooperative mechanisms was performed in [15]. The model which we adopted for the mouse ventricular myocyte contraction (Model 4 from [15]) gave the best fit to the existing experimental data for mice. The contraction model parameters for epicardial and endocardial cells are presented in the Supporting Information (Appendix S1).
Contraction force F contr (in mN/mm 2 ) was calculated using the equation [15]: where P2 max~f P3 max~f 01 f 12 f 23 S , ð6Þ In equation (1), F contrn is the normalized contraction force, and the coefficient 273.26 is obtained from fitting absolute values of the steady-state and dynamic experimental forces. For simulation steady-state force-calcium relationships (F-Ca), we used fixed values of the sarcomere lengths (SL), so that d(SL)/dt = 0, and changed intracellular Ca 2+ concentration. We simulated F-Ca relationships for sarcomere lengths 1.9, 2.1, and 2.3 mm. In this case, F contrn has time-independent magnitude.
For simulation twitch contraction, where F contrn is timedependent, we used Hooke's law, the linear relationship between contraction/relaxation force and cell shortening/extension: where SL 0 is the initial value of SL. In this case, sarcomere length becomes a function of time. To compare simulation data with experiments for cellular contraction, we estimated the variable cell length by where initial cell length L0 = 100 mm. For all simulations, we used extracellular Ca2+ concentration [Ca2+]o = 2 mM. The electromechanical cardiac cell models were stimulated with different frequencies using a stimulus current (I stim = 80 pA/pF, t stim = 0.5 ms) for at least 200,000 ms to reach a quasi-steady state. Simulated data of intracellular [Ca 2+ ] i transients, myocyte contraction force F contr , and sarcomere length SL on the interval from 192,000 to 200,000 ms were compared to extensive experimental data.
The model consists of 51 ordinary differential equations (see Appendix S1). Differential equations are solved by fourth-order Runge-Kutta method with time step 0.0001 ms. The model was implemented as an original Intel FORTRAN 90 code, which was run under SUSE Linux on a Dell Precision Workstation T3500 (Intel Xeon Processor W3670, 3.20 GHz, 8 GB RAM).

Steady-state Force-calcium Relationships
We first simulated steady-state force-calcium relationships. Both epicardial and endocardial cell models demonstrated the same simulation data for the steady-state force, as the contraction model Rice et al. [15] depends on intracellular [Ca 2+ ] I concentration. Figure 2A shows the Ca 2+ -dependence of the absolute value of contraction force obtained by Prabhakar et al. [23] for two different sarcomere lengths, 1.9 and 2.3 mm, from skinned mouse ventricular myocytes. For both cases, the force represents an increasing sigmoid function of calcium concentration. There is a relatively small increase in the saturation force from 48.8 to 57.2 mN/mm 2 when sarcomere length increases by about 20%, from 1.9 to 2.3 mm. Figure 2B shows simulation of the steady-state force-calcium relationships for three sarcomere lengths, 1.9, 2.1, and 2.3 mm. Our model is able to closely reproduce the saturating value of the force for corresponding sarcomere lengths. However, there are some differences between simulated and experimental data in sensitivity to external Ca 2+ , as simulated force saturates at smaller values of Ca 2+ concentrations. Such differences are due to a decrease in Ca 2+ sensitivity of skinned compared to intact cardiac cells [24].
Our model is also able to reproduce a shift in Ca 2+ sensitivity for steady-state force-calcium relationships shown for three sarcomere lengths (Fig. 2D). Such a shift can be clearly seen for normalized steady-state force-calcium relationships. Simulations show that an increase in sarcomere length leads to smaller half-saturation values of Ca 2+ concentrations, demonstrating an increase in Ca 2+ sensitivity (Fig. 2D). A similar shift in Ca 2+ sensitivity is also observed experimentally for mouse cardiac cells (Fig. 2C) [23,25].
In addition to the skinned mouse ventricular myocytes, our simulation data is also compared to the available experimental data on steady-state force-calcium relationships from intact cells, shown in Fig. 2C and 2D with unfilled circles [26]. Figure 2D shows that our simulations are in good agreement with the experimental data. IC 50 and Hill coefficient h obtained by fitting steady-state force-calcium relationships from McCloskey et al. [26]

Dynamic Behavior of Contraction Force
To test the ability of our model to reproduce the time behavior of the contraction force developed by mouse ventricular myocytes, we first stimulated the model cells with a constant frequency of 0.5 Hz. The time course of force in epicardial and endocardial cell simulations is plotted in Fig. 3A by red solid and dashed lines, respectively. As endocardial cells show larger [Ca 2+ ] i transients than epicardial cells, we obtained that the former develops stronger contraction force and larger shortening than the latter. The time behavior of the contraction forces obtained experimentally is shown by black solid lines with symbols [18,27,28,29].
There are significant differences in the experimental data obtained from different experimental groups on the time behavior of force, both in peak values and residual forces ( Table 1). Comparison of the time behavior of normalized simulated and experimental forces, both for epicardial and endocardial cells, shows a clear similarity in the time-to-peak values and relaxation of the simulated forces (Fig. 3B) [18,27,28,29,30].
Our model includes changes in sarcomere length during myocyte contraction. The time behavior of normalized sarcomere shortening for simulated cells is shown in Fig. 3C by red solid and dashed lines for epicardial and endocardial cells, respectively. The models do not show large differences in time-to-peak shortening and relaxation times. They closely reproduce myocyte shortening obtained in different experiments with mice (solid lines with  [44], Kogler et al. [28], and McCloskey et al. [29] (lines with symbols). For comparison, the initial sarcomere length in the model simulation is set to 2.1 mm, extracellular [Ca 2+ ] i concentration is 2 mM, and the frequency is 0.5 Hz, the frequency used most in the experimental data (see Table 1 Fig. 3C) [31,32]. For comparison of the time scales of contraction force and Ca 2+ dynamics, we also plot the time courses of the simulated and experimental intracellular Ca 2+ transients by red lines and black solid lines with symbols in Fig. 3D, respectively. In each case, there is a delay in force development following the peak of the Ca 2+ transient (compare times to peaks in Fig. 3B and 3D).

Force-frequency Relationships
In order to investigate force-frequency relationships, we also stimulated model cells with different frequencies ranging from 0.25 to 2.0 Hz. Frequency dependences of intracellular Ca 2+ transients, contraction force, and cell shortening are shown in Fig. 4. Our simulated peak [Ca 2+ ] i -frequency relationship (red solid and dashed lines in Fig. 4a) is within the variability of experimental data (solid lines with symbols in Fig. 4A) [18,29,33,34]. Note, that the simulated amplitudes of [Ca 2+ ] i transients for epicardial and endocardial cells are verified by the experimental data obtained by Dilly et al. [35] (Fig. 4D). The models are able to reproduce peak contraction force-frequency relationships for mouse ventricular myocytes in the frequency range from 0.5 to 2.0 Hz (Fig. 4B). The experimental data shows biphasic behavior of the peak force, with a decrease from 0.25 to 0.5-1.0 Hz, followed by an increase from 1.0 to 2.0 Hz [29], with a clear minimum in force-frequency relationships (however, see data of Ito et al. [34] were the minimum is less apparent). Our model reproduced such biphasic behavior of the force-frequency relationships for epicardial cells. Peak contraction force for endocardial cells increases with stimulation frequency.
Finally, we are able to simulate peak lengthening-frequency relationships (red lines in Fig. 4C). While some experimental data shows consistent decrease in cellular shortening with frequency [32], other data follows biphasic behavior [33,34] (solid lines with symbols in Fig. 4C). Our modeling data demonstrates biphasic behavior in cell shortening for epicardial cells, which is consistent with the biphasic behavior of the contraction force and [Ca 2+ ] i transients (red solid lines in Fig. 4A, 4B, and 4C). Model endocardial cells show only an increase in cell shortening as well as in [Ca 2+ ] i (red dashed lines in Fig. 4A and 4C, respectively).
Simulated time courses for contraction forces, sarcomere lengths, and sarcomere shortenings for three different resting sarcomere lengths (1.9, 2.1, and 2.3 mm) for epicardial and endocardial cells are shown in Fig. 5. As seen from the figure, an increase in the resting sarcomere length increases twitch force and relative sarcomere shortening. Similar behavior is also observed experimentally and from the simulations of others [9,11]. At comparable sarcomere lengths, the endocardial cells develop larger contraction force and sarcomere shortening than the epicardial cells (Fig. 5).

Constant versus Variable Sarcomere Length
While steady-state simulations show that peak force is dependent on the initial sarcomere length, there is also a dynamic relationship between force and sarcomere length. Our models use a variable SL when calculating the transition rate from nonpermissive to permissive states, as well as in the detachment rates in permissive states. To see the effect of using a variable SL in the transition rate equations, we ran simulations in which a constant SL replaced the variable SL in the calculation of the normalized sarcomere length SL norm~S L{1:3mm 2:3mm{1:3mm ð11Þ which is used in the detachment rates and transition rates in Markov model (Fig. 1B) g 10SL~gxbSL , g 21SL~2 g xbSL , g 32SL~3 g xbSL , ð12Þ g xbSL~gminxb 2{ SL norm ð Þ 1:6 , ð14Þ Ntm~5z3SL norm , where K Ca = k 2 ltrpn /k + ltrpn , and the constants can be found in the Supporting Information (Appendix S1).  Figure 6A shows force development in epicardial cells at a stimulation rate of 1 Hz in the simulation using a constant SL (dashed line) versus the simulation using a variable SL (solid line) (see also Fig. 6B). Data for endocardial cells displays similar behavior and is shown in Fig. 6C and 6D. The peak force when using a constant SL is clearly higher, while the residual force appears to be about the same. However, simulations run at various frequencies show that the peak and residual force when using a constant SL (Fig. 7E) is always higher than corresponding forces when a variable SL is used (Fig. 7C). Even though there is a difference in the magnitude of force, the frequency dependence of peak force when using a constant SL (black dashed line in Fig. 7F) is similar to the frequency dependence when a variable SL is used (black solid line in Fig. 7F). For comparison, Figs. 7B and 7D show simulated data on cell shortening and contraction force at different stimulation frequencies for endocardial cells, using variable sarcomere length (data on constant SL is not shown). As seen from the figures, both peak contraction force and cell shortening are larger for the endocardial cells than the epicardial cells.
In both cases, constant and variable SL, we observed a decrease in time-to-peak and time to 50% relaxation rate for the contraction force with an increase of stimulation frequency starting from 0.5 Hz. A similar increase in the residual contraction force at the larger stimulation frequencies is also observed experimentally [36].

Frequency Dependence of dL/dt and dF/dt
The frequency dependencies of the peak force and cell shortening are shown in Fig. 4. As might be expected, dL/dt and dF/dt also show frequency dependence. Simulated time courses for dL/dt (Fig. 8A and 8B) and dF/dt (Fig. 8C and 8D) are shown for various frequencies from 0.25 Hz to 4.0 Hz, both for epicardial ( Fig. 8A and 8C) and endocardial ( Fig. 8B and 8D) cells. A negative dL/dt value indicates cell shortening during a  contraction, while a positive dL/dt corresponds to relaxation. A positive dF/dt indicates the increase in force during a contraction, while a negative dF/dt corresponds to relaxation. The epicardial cell demonstrates a monotonic increase in the magnitudes of peak values for dL/dt and dF/dt in the frequency range from 0.25 to 4 Hz (Fig. 8A and 8C). In contrast, the endocardial cell shows a biphasic behavior in the peak magnitudes of the derivatives: an increase when the stimulation frequency changes from 0.25 to 2 Hz, and a decrease in the frequency range from 2 to 4.0 Hz (Fig. 8B and 8D).
The frequency relationship for +dL/dt max (solid lines) and 2dL/dt max (dashed lines) is shown in Fig. 9A. Both values show biphasic behavior. For the epicardial cell, +dL/dt max and 2dL/ dt max decrease at stimulation frequencies from 0.25 to 0.5 Hz, and then increase for stimulations frequencies up to 4 Hz. For the endocardial cell, +dL/dt max and 2dL/dt max increase at stimulation frequencies from 0.25 to 2.0 Hz, and then decrease for stimulations frequencies up to 4 Hz. When compared to the experimental data, our model tended to show, on average, peak contraction rates approximately equal to experimental data (open symbols) [32,37,38]. However, the model showed somewhat slower relaxation, thus lower values of +dL/dt max , than experimental data (solid symbols) [32,37,38]. Figure 9B shows the frequency relationship for +dF/dt max (solid lines) and 2dF/dt max (dashed lines). As with corresponding values for +dL/dt max and 2dL/dt max , the +dF/dt max and 2dF/dt max show biphasic behavior for both epicardial and endocardial cells.
To compare experimental and simulated data quantitatively, we plotted experimental and simulated results on time-to-peak and time-to-50% relaxation of the contraction force and intracellular [Ca 2+ ] i transients in Fig. 10. Simulated data are shown for both epicardial and endocardial cells (black and red, respectively, in Fig. 10B and 10D). Simulated data for time-to-peak force shows good agreement with the experimental data (compare Fig. 10B  and 10A), while time-to-50% relaxation are somewhat longer in the simulated data than those obtained in the experiments (compare Fig. 10D and 10C). Experimental data for time-to-peak and time-to-50% relaxation of [Ca 2+ ] i transients are somewhat longer than those from simulations, but the simulated time-to-50% relaxations approach the experimental values at larger frequencies. Epicardial and endocardial cells show similar simulated values for time-to-peak and time-to-50% relaxation of [Ca 2+ ] i transients, and for time-to-50% relaxation of contraction force. However, there are moderate differences between the cells for time-to-peak of the contraction force (Fig. 10B). The simulation data with constant SL is shown only for epicardial cells, as the data for endocardial cells is similar. The frequency dependence of force for an epicardial cell when a variable SL parameter is used is not as pronounced as the frequency dependence of force when a constant SL parameter is used (C and E). The initial SL for each simulation is 2.1 mm, but the residual force for higher frequencies leads to significant shortening (A and B). Frequency dependence of peak force for epicardial and endocardial cells with variable SL and for epicardial cells with constant SL is shown in (F). doi:10.1371/journal.pone.0063141.g007

Discussion
In this paper, we developed a new model for mouse ventricular myocyte contraction. This model is based on our previously published models for epicardial and endocardial cells [19,20,21,22], which includes a comprehensive description of action potential, ionic currents, and Ca 2+ dynamics. For a description of myocyte contraction, we adopted Model 4 developed by Rice et al. [15] by fitting experimental data on contraction for mice.
Mice demonstrate considerably faster heartbeats than many other species. Their contraction rate is about 10 beats per second [39], which is, for example, faster than the rabbit (4 Hz, [40]) and human (1 Hz, [41]) heart contraction rates. In addition, the action potential duration in mouse ventricular myocytes is also considerably shorter (APD 50 , 4.5 ms in mice [19] versus ,200 ms in rabbits [9] and ,300-400 ms in humans [42]). These differences suggest different time characteristics for contractions in mouse, compared to human or rabbit, ventricular myocytes.
In a mouse cardiac cell, at moderate stimulation rates, an increase in action potential is followed by an increase in [Ca 2+ ] i and a delayed increase in force. The peak value of Ca 2+ transient occurs after almost complete repolarization of action potential. In addition, peak contraction force appears after a significant decline of [Ca 2+ ] i . Our model replicates this relationship. Figure 11 shows normalized values for epicardial action potential (solid line), [Ca 2+ ] i (dashed line), and force (dotted line) over a 0.5 second interval for a simulation at 1 Hz. In larger species, such as rabbit, time scaling of the action potential, [Ca 2+ ] i and contraction force transients is different (Fig. 9 in [9]). For rabbits, [Ca 2+ ] i transient, in significant part, overlaps with the action potential and contraction force transient, while the peak sequence is the same as in mice.
Mouse ventricular myocytes, unlike other species, demonstrate biphasic frequency dependence of intracellular [Ca 2+ ] i transient and peak force [18,29] (however, see data of Ito et al. [34] were biphasic behavior is less apparent). Stuyvers et al. [18] suggested a qualitative mechanism which explains this biphasic behavior based on frequency-dependent Ca 2+ dynamics. The minimum occurs at the crossroad of the descending frequency trend of the Ca 2+ load into the sarcoplasmic reticulum during diastole and ascending trend in Ca 2+ entry into the cell through L-type Ca 2+ channels. They used a simplified description of Ca 2+ dynamics for mouse ventricular myocytes. Our model for an epicardial cell, which includes a comprehensive description of the electrical activity and Ca 2+ dynamics in mouse myocytes during cell twitch, was also able to reproduce this physiological phenomenon. In our model, myocyte contraction force is related to Ca 2+ dynamics through the Markov model for crossbridge kinetics. While both peak [Ca 2+ ] i transients and peak contraction force show minimum values as functions of stimulation frequency, these minimum frequency values are slightly different (Fig. 4). This trend is also confirmed by the experimental data of McCloskey et al. [29]. However, our model for the endocardial cell does not show biphasic behavior in the frequency-dependence of both peak [Ca 2+ ] i transients and peak contraction force. There are also some experimental data in which non-monotonic increase in peak [Ca 2+ ] i transients and myocyte shortening in mice is less apparent: even saturation and decrease in myocyte shortening amplitude at relatively large stimulation frequencies occur [34]. Our model for the endocardial cell, at least qualitatively, reproduced saturation and even decrease in sarcomere shortening and contraction force amplitude at 4-Hz stimulation ( Fig. 7B and 7D). This effect can be explained by the larger peak and diastolic values of [Ca 2+ ] i transients in endocardial cells compared to epicardial cells, which shift the operation interval of intracellular Ca 2+ towards a smaller slope in force-calcium relationships (Fig. 2D). Figure 9. Frequency dependence of dL/dt max and dF/dt max . (A) The simulated frequency dependences of (+dL/dt) max (solid lines) and (2dL/ dt) max (dashed lines). Experimental data from Chu et al. [37], Flagg et al. [38], and Huang et al. [32] are shown by symbols. We consider (2dL/dt) to correspond to cell shortening. (B) The simulated frequency dependence of (+dF/dt) max (solid lines) and (2dF/dt) max (dashed lines). We consider (+dF/ dt) to correspond to contraction. The initial SL for the simulations in (A) and (B) is 2.1 mm. Data for epicardial and endocardial cells are shown in black and red, respectively. doi:10.1371/journal.pone.0063141.g009 While there are no specific experimental studies of contraction force and cell shortening in mouse epicardial and endocardial ventricular myocytes, there are a few studies of the differences in action potentials and Ca 2+ handling in these cells [35,43]. The studies show that the endocardial cells demonstrate significantly larger [Ca 2+ ] i transients, and our modeling predicts larger contraction force and shortening in these ventricular myocytes.
Our electromechanical model for mouse ventricular myocyte contraction includes a variable sarcomere length during cell contraction, the effect that occurs in most experiments. Simulations with variable sarcomere length produce significantly smaller contraction force than the simulations with constant sarcomere length despite the same time course and amplitude of [Ca 2+ ] i transient during twitch. This suggests the importance of the inclusion of cell shortening in the model for cardiac myocyte contraction. Note that a similar result was obtained with a more complex model of Rice et al. [9], developed for rabbit ventricular myocytes, who also studied the effects of variable and fixed sarcomere length on the force development.
Several models for cardiac myocyte contraction have been developed to date [7,9,11,15,17] (see also review [13]). Earlier models did not include sarcomere shortening during twitch [7,15,17]. They are primarily focused on simplification of the description of crossbridge kinetics, their dependence on Ca 2+ dynamics, and careful reproduction of the existing experimental data on steady-state and dynamic force-calcium relationships. Most of these models have limitations due to this and other simplifications.
Rice et al. [15] investigated five Markov models describing contraction mechanisms in cardiac myocytes. Two of the models consisted of four tropomyosin states and transitions between them (N0, N1, P0, and P1, see Fig. 1B). These models differ by the mechanisms of modulation of the transition rates (in Fig. 1B they are defined as k NP and k PN ). In Model 1, rates k NP and k PN are independent of the developed force, while in Model 2 the rates depend on the developed force. In both models, Ca 2+ binding to troponin directly affects tropomyosin shifting, i.e., rates k NP and k PN . Model 3 includes an indirect connection of the Ca 2+ binding to troponin and tropomyosin shifting, as shown by dashed arrows in Fig. 1B (see also [15]), and only four states (N0, N1, P0, and P1). Models 4 and 5 were extended to up to three crossbridge bindings, which resulted in four permissive tropomyosin states, P0, P1, P2, and P3, (Fig. 1B and [15]). The only difference between Models 4 and 5 is the modulation of the k 2 ltrpn rate by generated force. Because Model 4 and Model 5 yielded an approximately equal description of myocyte contraction, we implemented Model 4 in our electrophysiological model, as Model 5 led to unstable solutions.
Our model of mouse ventricular myocyte contraction also has some limitations due to the simplification of the biophysical mechanism of contraction. In particular, the model uses a simplified description of the relationships between contraction force and cellular shortening in the form of Hook's law, while the real dependence is more complicated [9]. It does not describe the effects of cellular shortening on Ca 2+ transients, as does the model of Rice et al. [9]; however, this effect is relatively small. Also, our model, as most other models, did not take into account intracellular spatial inhomogeneities of Ca 2+ concentration and crossbridge binding sites.
Nevertheless, despite the limitations, our electromechanical model of mouse ventricular myocyte contraction was extensively verified by experimental data obtained for mice. It reproduced reasonably well a significant amount of the existing experimental data. The model can be used for cells from two different regions of the heart (epicardium and endocardium). As with most other models, it uses a simplified description of the contraction force generation. We employed a six-state Markov model for tropomyosin dynamics and separate Ca 2+ binding to troponin (Fig. 1B) to describe force development. More comprehensive models will be necessary to develop a better simulation of more extended experimental data sets.

Supporting Information
Appendix S1 Model Summary.