Computationally efficient model of myocardial electromechanics for multiscale simulations

A model of myocardial electromechanics is suggested. It combines modified and simplified versions of previously published models of cardiac electrophysiology, excitation-contraction coupling, and mechanics. The mechano-calcium and mechano-electrical feedbacks, including the strain-dependence of the propagation velocity of the action potential, are also accounted for. The model reproduces changes in the twitch amplitude and Ca2+-transients upon changes in muscle strain including the slow response. The model also reproduces the Bowditch effect and changes in the twitch amplitude and duration upon changes in the interstimulus interval, including accelerated relaxation at high stimulation frequency. Special efforts were taken to reduce the stiffness of the differential equations of the model. As a result, the equations can be integrated numerically with a relatively high time step making the model suitable for multiscale simulation of the human heart and allowing one to study the impact of myocardial mechanics on arrhythmias.


Introduction
Mathematical modelling of myocardial electromechanics provides a useful tool for understanding processes of excitation, propagation of action potential (AP), Ca 2+ -activation, and contraction in the normal and diseased heart. The modelling approach is particularly important because apart from direct excitation-contraction coupling, there are mechano-electrical and mechano-calcium feedbacks that make the system particularly complex for experimental studying. Detailed electrophysiological models of human ventricular or atrial cardiomyocytes and their clinical applications were reviewed recently [1,2]. Some of these cell-level models [3][4][5][6] reproduce accurately the time-courses of various ionic currents, Ca 2+ transients, and the dependency of the AP duration and amplitude on the stimulation frequency. Such models allowed their authors to simulate various arrhythmias in the human heart and their drug or surgical treatment [1,2,7,8]. One can combine electrophysiological models with a description of myocardial mechanics including both passive elastic and active Ca 2+ -dependent stress. The electromechanical models, for example [9][10][11], allow one to examine numerically the effects of some molecular processes on disorders in muscle contraction or its excitation and can be incorporated into multiscale models of the human cardiovascular system for the simulation of the whole heart or its chambers contraction. The combined models can be used to simulate the pumping function of the heart at different conditions [12] as well as they can be applied to the numerical investigation of electrophysiological activity of the heart [13]. However, the authors of the combined electromechanical models usually did not show any results relating to the influence of electrical stimulation on myocardial mechanics. It is unknown whether the models are capable to reproduce important dependencies of the twitch force on the stimulation frequency, including Bowditch effects and the force response to the variation in the interval between stimuli. The effect of the strain of myocardial tissue on the AP propagation velocity, so-called mechano-electrical feedback, was also not taken into account in most of those models with a few exceptions. The attempts of the modelling of the strong electromechanical coupling including the mechano-electrical feedback were made by a number of research groups [11,[14][15][16][17][18][19], however, the assumptions used in these models were not always supported by experimental observations, or only the mechano-calcium feedback was considered while a direct influence of strain on the speed of activation propagation was not taken into account. Another important feature of the force-frequency relation in cardiac muscle that was not demonstrated by previous simulations is the acceleration of force relaxation at increased stimulation frequency [20]. This acceleration is an important feature of the cardiac muscle as it provides sufficient time for refilling heart chambers even during shortened diastole at high heartbeat frequency.
Electromechanical models with a detailed description of ionic currents are usually expensive computationally as they require a very small time-step for numerical integration of kinetic equations. In contrast, phenomenological models of cardiomyocyte electrophysiology such as the Aliev-Panfilov model [21] are less demanding computationally although allow one to reproduce the AP propagation during normal heartbeats and arrhythmias [22][23][24].
Here we present a model that combines the Aliev-Panfilov model [21] with an improved version of our model of myocardial mechanics [25] extended by the introduction of major calcium currents and the additional balance equations for the Ca 2+ concentrations in the cytosol and sarcoplasmic reticulum (SR) based on the approach suggested in [3]. The mechanical model describes various mechanical phenomena including steady-state force-velocity and stiffness-velocity relations, tension responses to strain or load changes, tension redevelopment, load-dependent relaxation, etc. [25]. Special efforts were made to simplify the equations and make them less stiff. Our model of myocardial mechanics and its regulation by Ca 2+ ions was also modified to reduce computational time and to account for new structural data concerning the interaction of regulatory and contractile proteins. The equations of Ca 2+ balance from the model suggested by ten Tusscher and Panfilov were also simplified using an asymptotic approach to make them less stiff. Although the simplistic electrophysiology model cannot reproduce the time-courses of all ionic currents, the combined electromechanical model reproduces the Ca 2+ -transients in different modes of cardiac muscle contraction. The model considers both mechano-calcium and mechano-electrical feedbacks through the effect of muscle strain on the Ca 2+ -Na + exchange (NCX) and the conductivity and capacitance of cell membranes, respectively.
Our electromechanical model reproduces the effects of the stimulation frequency and interval between stimuli on the twitch force and intracellular free Ca 2+ concentration as well as instantaneous and slow response of twitch force to muscle stretch. Together with its tolerance for relatively high time step for numerical integration, this makes the model suitable for multiscale simulation of the heart including studying the strain-dependence of arrhythmias when the strain and interstimulus interval vary in space and time.

Model description
Electrophysiology and mechano-electrical feedback. The Aliev-Panfilov model [21] was used for the simulation of the membrane potential variation where u is the dimensionless membrane potential related to the potential V: ; v is a phenomenological dimensionless variable that describes the kinetics of ionic currents, and I stim is a dimensionless stimulation current; D ij are the contravariant components of the conductivity tensorD (there is summation over indexes i and j in Eq (1)); r j is a covariant derivative over j-th spatial coordinate in a Lagrangian (material) coordinate system. Values of time scale τ and dimensionless constant parameters k, a, ε, μ 1 , μ 2 are listed in Table 1. In transversally isotropic media, tensorD has the following form whereÊ is the unit tensor,B ¼ l ! � l ! is the anisotropy tensor, and l ! is the unit vector directed along muscle fibres in strained muscle; d 0 and d 1 are isotropic and anisotropic conductivities (d 1 is the longitudinal conductivity along the fibre direction, d 0 is the transversal conductivity).
In order to account for the slowdown of the action potential (AP) propagation upon muscle stretch [26][27][28], the normalized capacitance of the cell membranes C was assumed being strain-dependent [19]. The strain-dependence was shown to be associated with a strain-dependent recruitment of caveolae into the cell membrane causing an increase in its capacitance of and the electrical time constant even when stretch-activated ionic channels were inactivated [27]. Following [19] we assumed that the cell membrane capacitance C normalized for its value at λ = 1 and d 1 depend on the axial strain l ¼ l s l s0 as follows: 8 > < > : Here the values of the constant parameters: Δ C , n C , K C , Δ d , n d , K d are the same as suggested by de Oliveira et al. [19]. The values of the parameters are given in Table 1.
In reality, the caveolae dynamics that controls the changes in the capacity and conductance are not instantaneous processes. Although the time-course of this process was not measured experimentally yet, it seems reasonable to assume that changes in C and d 1 caused by strain are significantly slower than the change in membrane potential during onset and decay of AP in cardiac muscle cells. Therefore, the first equation in Eq (1) can be rewritten as Calcium transport across cell membrane, mechano-calcium feedback. Only two major components of Ca 2+ transport across the membrane of cardiomyocytes were taken into consideration: current through L-type Ca 2+ -channels and calcium-sodium exchange, I CaL , and I NCX , respectively. We neglected any fast and slow gating of the L-type Ca 2+ -channels and simply assumed the flux to be a function of the membrane potential u: where G CaL , u 0 , u � , k u , are constants specified in Table 1. To simplify calculations, the term k CaL ðuÀ uÞ ðexpðk CaL ðuÀ uÞÞÀ 1Þ ¼ p 0 þ p 1 u þ p 2 u 2 in Eq (5) with nominator and denominator approaching zero at u approaching u � was approximated with a second-order polynomial with an accuracy of 3%. Coefficients of the polynomial p i are also listed in Table 1.
The NCX was defined as suggested by ten Tusscher and Panfilov [3] where  In our phenomenological model, the sodium concentration in the cytoplasm was not a variable that obeys a detailed equation for Na + balance as in [3,29]. Instead, it was set as a function of sarcomere length.
where Na i0 = [Na + ] i0 and k NCX are constants specified in Table 1. Such a simple assumption within the framework of a phenomenological model is in line with the results of more detailed modelling of the slow response [30]. The slow response is a graduate increase in twitch amplitude upon a stretch of a cardiac muscle. These and other authors concluded that that the strain dependence of Na + -H + exchange is a major contributor to the slow response [31,32]. We show here that the simplest assumption (6A) is sufficient for explaining the amplitude and the time course of the slow responses of the twitch amplitude and Ca 2+ transients to an increase in muscle preload.
Calcium transport into and from SR. Ca 2+ uptake I up from the cytoplasm to the sarcoplasmic reticulum (SR) was defined similarly to that in [3,33]: where G up , K up are constants specified in Table 1. We have added a variable factor p that is the level of phosphorylation of protein(s) controlling the Ca 2+ uptake to SR by sarcoplasmic reticulum Ca 2+ -ATPase (SERCA) to account for the acceleration of muscle relaxation at high heartbeat rates.
There are several possible candidate mechanisms of posttranslational modification and target proteins responsible for the acceleration of the Ca 2+ uptake to SR upon an increase in the heartbeat rate. Phosphorylation of a small protein phospholamban (PLN) associated with Ca 2+ pump SERCA was the first candidate for frequency-dependent muscle relaxation [34]. Dephosphorylated PLN was shown to inhibit SERCA while its phosphorylation loses the inhibition so that the rate of Ca 2+ uptake to SR increases. Two PLN residues, Ser16 and Ser17, are phosphorylated by different kinases. The second one, Ser17 is phosphorylated by calmodulindependent kinase II (CaMKII) enhanced by high stimulation frequency [35]. Calmodulin is a Ca 2+ binding protein, therefore, the level of CaMKII phosphorylation should increase with an increase in time-average Ca 2+ concentration in cytosol. Later, some CaMKII-dependent acceleration of cardiac muscle relaxation was found in PLN-knockout animals [36] showing that other CaMKII-phosphorylated protein(s) including SERCA itself might be involved in the frequency dependence of cardiac muscle relaxation.
Although the precise mechanism of accelerated relaxation at high stimulation rate is still debated [37,38], the presence of such mechanism, which most probably involves CaMKII, is well established. We do not specify phosphorylation of which particular protein(s) is responsible for the acceleration of the Ca 2+ uptake into SR in our phenomenological model. Instead, we simply specify the dependence of its phosphorylation level p on cytoplasmic Ca 2+ concentrations c as follows: where k p , K p are constants characterizing the phosphorylation rate and the equilibrium constant, respectively; values of both parameters are presented in Table 1.
The major source of Ca 2 in the cytosol is calcium-induced calcium release (CICR) from SR. CICR is a complex process activated by increased Ca 2+ concentration in the narrow subspace between the cell membrane with the L-type Ca 2+ -channels (dihydropyridine receptors; DHPR) and the SR membrane with ryanodine receptors (RyR) as well as by high Ca 2+ concentration in the SR lumen. Several sophisticated models partially reviewed by Hinch et al. [39] as well as a non-Markovian one [40] and a model of interacting clusters of DHPR and RyR [41] were suggested for simulation of different aspects of CICR. Here we neglected all details that are not necessary for describing major feature of CIRC in normally contracting cardiac muscle cell and set the CICR flux in a rather simple form [3]: The kinetics of closing gate R was described by the equation Here G rel , K rel , K 2 , k 4 , K R are constants specified in Table 1. Eqs (9-10) are based on the assumption that CIRC activation is very fast although having relatively low sensitivity for Ca 2+ concentration in the narrow subspace between L-type Ca 2+ channels and RyR, c SS , while CIRC deactivation is much slower although more sensitive to the Ca 2+ concentration in the subspace. In contrast to [3] we assume that neither the CIRC activation nor its inactivation is sensitive to the amount of Ca 2+ in SR as the model was able to describe the frequency-dependence of twitch amplitude without such assumption (see below). Besides, we have modified Eq (10) by taking into account a relatively high sensitivity of RyR to the Ca 2+ concentration in subspace c ss .
Calcium-binding to Tn and kinetics of Tpm-Tn regulation. Muscle contraction is activated by binding of Ca 2+ ions to troponin C (TnC), a subunit of troponin (Tn) complex of three regulatory proteins. The Ca 2+ binding leads to conformation changes causing a shift in the position of another regulatory protein tropomyosin (Tpm), which blocks the binding sites for myosin on actin in the absence of Ca 2+ and opens them upon Ca 2+ binding to TnC. Following Fusi et al. [42] who have shown that Ca 2+ binding to TnC is fast, we assumed that the Ca 2+ binding to TnC in the blocked or activated (closed plus open) states of the regulatory units of the thin filaments is quick and reversible although the affinities of TnC for Ca 2+ are different. Accordingly, the fractions of Ca 2+ -bound TnC molecules in the activated and blocked regulatory units were set as follows: where B and A refer to the blocked and activated states of the regulatory units, respectively; index i = 1, 2 corresponds to the overlap and non-overlap zones of the thin and thick filaments in sarcomeres, respectively. The activated state is a combination of the closed and open states of the Tn-Tpm complex according to the three-state model of McKillop, Geeves [43]. The equilibrium constants K A , K B (given in Table 1) are different as TnC in the activated regulatory units bind Ca 2+ tighter than in the blocked ones [44]. The rate-limiting step of the Ca 2+ activation of a striated muscle is the transition from the blocked to activated state [42]. The transition requires the detachment of the inhibitory domain of the troponin subunit troponin I (TnI) from actin and Tpm and binding of TnI switch segment to the hydrophobic pocket on TnC that opens upon Ca 2+ binding to TnC [45,46].
The kinetics of activation of the regulatory units in the overlap A 1 and non-overlap A 2 zones of the thin filaments was described by equations: @WðlÞ @t Where constant parameters k A , κ, ξ, k l , k n1 (k n2 = 0) and the length of an actin filament in a sarcomere l a are specified in Table 1. W, W max = 0.5(l m −l b ) are the normalized length of the overlap zone of the thin and thick filaments in sarcomere and the maximal length of this zone; l m and l b are the length of the thick filament and its bare zone. These equations are similar although different from those in [25]. Parameters κ, ξ characterize cooperativity between activation of neighbor regulatory units. TnC binds Ca 2+ in both the blocked and activated regulatory units although with different affinities as suggested by [44]. In addition, the Eq (12) is less stiff than that used in [25] and therefore more convenient from the computational point of view. It should be noted that the last "convective" term in Eq (12) can be omitted if the strain rate of sarcomere normalized for the sarcomere length is small compared to the rate constant of Ca 2+ activation k A . Eq (12) then reduces to Calcium-mechanical coupling. We slightly modified the system of the equations for calcium-mechanical coupling based on our model of cardiac cell mechanics [25]. The mechanical model described the kinetics of generation of active tension T a by contractile proteins myosin and actin and Ca 2+ activation of their interaction. We modified the model equations to reduce the stiffness of the system and to describe the force-Ca 2+ relationships more accurately. The expression for active tension was the same as in the previous model.
while the equations for the kinetics of the interaction of the contractile proteins were slightly modified: Here n, δ, A 1 are the fraction of the actin-bound myosin heads, ensemble-averaged distortion of actin-bound heads per sarcomere, and the normalized activation of regulatory tropomyosin-troponin (Tpm-Tn) strand in the overlap zone of sarcomeres, respectively. Functions  F(δ), G(δ), W(λ), θ(δ) and constants E, N M , N xb , h, k 01 , l s0 were the same as suggested previously in [25], and new constant δ � is specified in Table 1. From the mechanical point of view, the threshold δ � defines a transition from viscoelasticity to plasticity at high stretch rate and reduces the stiffness of the differential equations at a high stretch rate. One should note that the binding of myosin heads in Eq (15) is controlled by A 1 squared instead of A 1 itself as it was in [25]. This assumption was motivated by recent structural data showing an interaction between two Tpm-Tn strands via the N-terminal part of troponin T (TnT) that binds the Tpm strand at the region of the overlap junction of two consecutive Tpm dimers arranged on the opposite sides of the thin filament [46]. Such interaction should induce cooperativity of the azimuthal movement of the two Tpm strands so that two regulatory units on both sides of the thin filament undergo the block-to-open transition together when both TnI on the opposite sides of an actin filament are released from actin.
Calcium balance. The calcium concentration in subspace c SS obeys the balance equation from [3] a ss @ @t where α ss �1, B ss , K ss are relative volume of subspace, Ca 2+ buffer concentration in subspace, and its equilibrium constant; G xfer is the rate constant of Ca 2+ diffusion from the subspace to the cytosol. The smallness of α ss makes the system of differential equations stiff, so that it requires a small time step and\or special numerical methods for its integration. On the other hand, the smallness of α ss allows one to substitute Eq (11) with an asymptotic solution at α ss !0 using Tikhonov theorem [47], i.e., to set the right-hand side of Eq (17) to zero provided that the solution of the reduced equation is a stable root of Eq (18).
Eq (18) is a cubic equation for c SS at constant or slowly changing c, c SR , and u: Its solution obtained by the Cardano formulas gives the explicit expression c SS (c SR , c, R, u) that was substituted to Eqs (9, 10) and the equations describing Ca 2+ balance in cytosol and SR (see below). Depending on the model parameters, Eq (13) can have one, two or three real positive solutions. For the case of three solutions, the smallest and highest solutions are stable, while the intermediate one is unstable. More details of the choice of a solution of Eq (19) are given in S1 File.
Both cytosol and SR contain Ca 2+ buffers different from TnC with the capacities (total concentration) B cyt and B sr , respectively, and the equilibrium constants K cyt , K sr , respectively (for values see Table 1). Both buffers were assumed to bind Ca 2+ very quickly and reversibly so that the equation of Ca 2+ balance in cytosol was formulated as follows: where α cyt is a relative volume of cytosol; G leak is the coefficient that characterizes Ca 2+ leak from SR to the cytosol. The total concentration of Ca 2+ bound to TnC is given by expression.
Therefore, the flux of free Ca 2+ from cytosol due to its binding to TnC I Tn is as follows: where C Tn , is the total concentration of Tn. Substituting Eq (22) into Eq (20) one can obtain the equation for Ca 2+ balance in the cytosol. Ca 2+ balance in SR was described by the equation.
Here α sr are relative volumes of SR in cardiac muscle tissue, B SR , K SR are the total concentration of Ca 2+ buffer in SR and its equilibrium constant, respectively.

Results
For simulation of contraction of a single cardiac muscle cell or homogeneously stimulated cylindrical multi-cell samples, the model reduces to a 0D problem specified by the set of ordinary differential equations: Eq (1) without diffusion term as the membrane potential does not depend on coordinates, Eqs (8), (10), (15), (16), (20), (22), (23). The mechanical problem setup and the total tension specification are described in S1 File. The system of the equations was solved numerically together with mechanical equations defining the load or strain of cardiac muscle using forward Euler method with the time step of 0.1, 0.2 and 0.05 ms. The time step of 0.1 ms provided good convergence (see S3 Fig in S1 File). In the simplest case of 0D homogenous isometric contraction at a constant sarcomere length l s , δ = 0 and Eq (17) can be omitted. Isotropic and anisotropic (titin) components of passive elastic stress were the same as in [25].

Steady-state Ca 2+ regulation, effect of sarcomere length
Steady-state dependencies of normalized active tension T A,norm , the normalised fractions of activated Tn-Tpm regulatory units in overlap zone and outside it A 1 , A 2 , and normalized concentration of Ca 2+ -bound Tn in the overlap and non-overlap zones of sarcomere at different constant sarcomere length ([CaTn] 1 , [CaTn] 2 ) are shown in Fig 1. The active tension T a was normalized for its steady-state value at saturating Ca 2+ concentration at full overlap of the thin and thick filaments in sarcomeres. The Ca 2+ -tension relation at different sarcomere length for the model (Fig 1(C), 1(D)) was similar to that observed experimentally in rat [48] and human [49] cardiac muscles. A decrease in the sarcomere length by 0.4 μm lead to a~40% increase in apparent equilibrium constant estimated from normalized Ca 2+ -force relation (Fig 1(D)). The Ca 2+ -force relation obtained from the model is similar to that predicted by the Hill equation although is more steep at low Ca 2+ concentration and less steep at high Ca 2+ concentration as

PLOS ONE
found by Dobesh et al. in [48]. The shift in the Ca 2+ sensitivity caused by myosin cross-bridges can estimated from the difference between A 2 1 and A 2 2 in Fig 1(D). It was similar to that obtained in experiments with substances which inhibit myosin binding to actin [44]. The Ca 2+ -sensitivity of active isometric tension for the model was higher than that observed in experiments with skinned muscle cells and multicellular specimens although was similar to those in intact trabeculae [50]. As the model was designed for describing intact human cardiac muscle, Ca 2+ binding constants for Tn were adjusted to fit available data obtained on intact human cardiac muscles.

Time course of model variables during twitch contractions, effect of sarcomere length
The time course of all model variables and major Ca 2+ flows during an isometric twitch contraction at sarcomere length 2.2 μm at a stimulation rate 1 Hz are shown in Fig 2. The time courses of tension and cytosol Ca 2+ concentration c were similar to those observed experimentally in human cardiac muscle [51]. Sharp instantaneous changes in c SS (Fig 2(A)) were caused by jumps of the Cardano solution of the cubic Eq (19) when system passed a bifurcation point. The difference in the time course of A 1 and A 2 is induced by more tight Ca 2+ binding to Tn in the overlap zone due a cross-bridge Tmp-Tn interaction. Outward Ca 2+ flow thorough NCX during diastole balances its inward flow and flow through the L-type Ca 2+ channels during systole. The total amount of Ca 2+ entering cytosol from SR was about twice higher than that entering cell during systole via NCX and L-type channels (Fig 2(B)).
The effects of sarcomere length on the parameters of twitch contractions are shown in Fig 3. The model describes an increase in the twitch amplitude and duration upon a muscle stretch. Apart from an instantaneous rise in twitch amplitude there was an additional~25% increase that develops for several tens of seconds (Fig 3(A)). A slow response very similar to that shown in Fig 3 was observed in cardiac muscle samples from human [52] and other mammalians [53]. Just after stretch, calculated twitch amplitude increased, while the amplitude of Ca 2+ transient slightly decreased. During the slow response, the peak Ca 2+ gradually increases to its pre-stretch value. The fast decrease in Ca 2+ concentration just after the peak accelerates after the stretch. The slower decrease in the Ca 2+ concentration in cytosol after the peak of the twitch decelerates upon the stretch as observed experimentally [54,55]. The graduate increase in the peak Ca 2+ concentration caused further increase in the twitch amplitude (Fig 3(B), 3 (C)). These changes were caused by more tight Ca 2+ binding to TnC at longer sarcomere length. In our model, the slow response was caused by an increased Ca 2+ influx into cardiac muscle cells via NCX due to a strain-dependent increase in intracellular Na + concentration. As the model does not include detailed description of Na + balance, there is a single parameter k NCX responsible for the slow response in our model.

Effect of stimulation frequency and interval between stimuli on twitch amplitude and duration
Changes in twitch amplitude upon a change in the stimulation rate are shown in Fig 4. The model reproduces so called Bowditch effect, i.e., an initial decrease in the twitch amplitude just after an increase in the stimulation rate followed by its graduate increase (Fig 4, top) and biphasic response to a decrease in the stimulation rate: an increase in the amplitude after first elongated interval and its graduate decrease in the course of subsequent low-frequency stimulation (Fig 4, bottom).
The dependence of the steady-state twitch amplitude on the stimulation frequency and on the duration of the interval after long time stimulation at a frequency of 1 Hz are shown in Fig  5. Two types of experiments were simulated: long-term periodical muscle stimulation at different constant frequency (Fig 5(A)) and the experiment when a stimulus is applied with different intervals after long-time stimulation at 1 Hz (Fig 5(B)).
The force-frequency dependence shown in Fig 5(A) is similar to that obtained for human cardiac muscle samples [20,56,57]. The frequency dependence of the peak Ca 2+ concentration (Fig 5(A)) was also similar to that observed in cardiac muscle samples from non-failing human hearts [51]. The increase in the twitch amplitude with an increase in the stimulation rate in our model was provided by two factors: i) an increase in the time-average Ca 2+ influx via the L-type Ca 2+ channels caused by an increase in the ratio of systole duration to the heartbeat circle duration that led to Ca 2+ accumulation in SR; and ii) an increase in Ca 2+ accumulation in SR due to increased Ca 2+ uptake rate by SERCA caused by a more active phosphorylation of phospholamban or other SERCA-associated proteins at higher average Ca 2+ concentration in cytosol. The model also describes the effect of extra systole (for intervals less than 1 s) and post pause potentiation at interval longer than 1 s (Fig 5(B)).
The time courses of normalized active tension at different stimulation rates are shown in Fig 6. The model describes a decrease in the twitch duration at higher stimulation rate observed experimentally in cardiac muscle specimens from different mammals including humans [20]. In our model, the acceleration of twitch relaxation with frequency is explained by an increase in the rate of Ca 2+ uptake into SR at higher time-averaged Ca 2+ concentration in cytosol due to the enhanced phosphorylation of SERCA-associated protein(s) defined by Eqs (7) and (8).
The results of simulation of isotonic twitch contraction at different load are shown in Fig 7. The lower was the load the higher was shortening amplitude. Importantly, shortening under low load accelerated muscle relaxation that became significantly faster than that for isometric contraction, showing so-called load-dependent relaxation of the left ventricular myocardium [58,59].

Aim of the work and its main results
The aim of the work was to build an electromechanical model suitable for multiscale simulation of the human heart. The main requirements for the model were: • its ability to simulate changes in major electrophysiological parameters: duration and propagation velocity of AP as well as instantaneous and slow mechanical responses to strain; • the ability of the model to simulate instantaneous and long-term changes in twitch amplitude and duration caused by variation of the interval between subsequent APs; • the ability of the model to describe all major mechanical properties of cardiac muscle such as force-velocity and stiffness-velocity relations, non-steady responses to changes in muscle length or load including load-dependent relaxation; • the model should be relatively simple computationally for its usage in the multiscale simulations.
A number of models with detailed descriptions of various ion currents, AP generation, and propagation have been suggested [1][2][3][4][5][6]. Mechano-electrical and mechano-calcium feedbacks causing changes in the time course of AP and Ca 2+ transients and the twitch amplitude were also taken into account by some detailed electromechanical models in 0D simulations [10,14,30,60]. A data-based model of the strain-dependence of the AP propagation velocity was suggested [19] although the possible role of this type of mechano-electrical feedback in arrhythmogenesis was not studied.
However, the electromechanical models with detailed descriptions of ion currents require a very small time-step and very fine space mesh for numerical simulation. This makes simulation of heart electromechanics with realistic geometry for several seconds prohibitory expensive computationally making the stimulation of several heartbeat cycles with varying intervals between them or during an arrhythmia with complex pattern of activation and contraction almost impossible. To avoid this problem, we decided to use a popular phenomenological model of cardiac electrophysiology [21] that describes the observed dependence of the AP duration on the stimulation frequency. The model was shown to be able to simulate various cardiac arrhythmias [15,24,61], The mechano-electrical feedback was simulated in our model according to de Oliveira et al. [19] with a slight modification to account for its anisotropy. Cardiac muscle mechanics, calcium-mechanical coupling, and mechano-calcium feedback were described using a modification of our model [25]. The modification accounts for new structural data concerning the molecular mechanics of Ca 2+ regulation of cardiac muscle contraction [46,62] and, importantly, makes the system of equations that describes the Ca 2+ -activation of the thin filaments less stiff. Besides, an equation that describes the kinetics of the actin-myosin interaction was modified to avoid singularity at the high speed of muscle stretch. The model of electro-calcium coupling was adopted from [3] and simplified by neglecting some minor Ca 2+ currents and by using an asymptotic form of the description of Ca 2+ concentration in the narrow subspace between the L-type Ca 2+ channels, the T-tubule, and ryanodine receptors in the SR membrane. A novel aspect of the model is the dependence of the ATPdependent Ca 2+ uptake into SR on the phosphorylation of protein(s) associated with SERCA that is believed to be responsible for the acceleration of muscle relaxation at high stimulation frequency [37].

Estimation of parameters and verification of the model
We used standard values of the parameters of the Aliev-Panfilov model [21] in Eq (1) and the same parameters of the mechano-electrical feedback in Eqs (3) and (4) as those suggested by de Oliveira et al. [19]. The parameters of the model for the Ca 2+ flows between cytosol, SR, and subspace as well as the parameters of the Ca 2+ currents via L-type Ca 2+ -channels in Eq (5) and Na + -Ca 2+ exchange in Eq (6) were the same as in [3]. The parameters of the actin-myosin interaction were the same as in [25] with an additional parameter δ � that describes the transition from viscoelasticity to plasticity upon muscle stretch at high velocity. The parameters of the calcium-mechanical coupling describing Ca 2+ binding to TnC in the blocked and activated states and the kinetics of the transition between these states (Eqs 11) and (13) were adjusted to fit the data on the Ca 2+ -tension relations at different sarcomere length [48,49] and the data on the effect of actin-bound myosin heads on the activation of the thin filaments [44]. A correction was made to account for the difference between Ca 2+ binding to Tn in skinned and intact cardiac muscles [50]. The parameter k NCX that describes the strain-dependence of the intracellular Na + concentration was adjusted to fit the slow response amplitude in human cardiac muscle samples [52]. The model was verified by comparison of the calculated frequency dependence of the isometric twitch amplitude (Fig 5) with data obtained on samples of human cardiac muscles [20].

Reduction of previous models
To make the model as simple as possible computationally, while capable to simulate all major phenomena of the excitation-contraction coupling in human cardiac muscle, we neglected all Ca 2+ transmembrane currents except those through the L-type Ca 2+ channels, I CaL , and the Na + -Ca 2+ exchange (I NCX ). Compared to [3] I CaL was set as a function of membrane potential according to Eq (5) assuming that the fast opening of L-type channels is instantaneous while their slow closure can be neglected. The Eq (6) for I NCX was the same as in the model by ten Tusscher and Panfilov [3]. As subspace volume is very small (α ss �1 in Eq (17)) in [3], it requires a very small time step for numerical integration. To avoid this problem, we employ Tikhonov theorem [47] and substituted differential Eq (17) with cubic algebraic Eq (19).

Mechano-electrical and mechano-calcium feedbacks
Some modern models describe so called mechano-electical feedback, i.e. the changes in the AP durations upon changes in muscle length. Such changes could be caused by changes in Ca 2+ concentration in the cytosol accompanying Ca 2+ binding or release from Tn [14]. Some models include a description for the stretch activated channels [10,16]. Some other models take into account the strain-dependence of myocardial electrical conductivity due to a change in geometry at constant specific conductivity [17,60,61]. A more detailed data-based model by de Oliveira et al. [19] describes a decrease in the AP propagation velocity upon muscle stretch caused by the strain-dependence of membrane capacitance and cell-to-cell conductivity caused by a recruitment of caveolae to the sarcolemma as found by Pfeiffer et al. [27]. Although this process cannot be instantaneous, we assume, as a first approximation, the capacitance and the conductivity to be functions of current sarcomere length according to Eq (3). This equation was the only mechano-electrical feedback in our model, stretch-activated channels were neglected. The mechano-calcium feedback was described by two terms: the strain-dependence of the transition from the blocked to the activated state of the Tn-Tpm complex in Eqs (12) and (13) and the strain-dependence of Na + concentration in cytosol, [Na + ] i , in Eq (6).
Although an increase in [Na + ] i is possibly caused by strain-dependent changes in Na + exchange, and other mechanism(s) can also be involved to slow response to stretch [30], we used the simplest approximation that can describe this phenomenon.

Novelty
The coupling of different models of cardiac mechanics and electrophysiology is a common approach for developing electromechanical models. We would like to highlight the following aspects of our model.
1. Our previous model of myocardium mechanics [25], which modification was used in the study, reproduced a large set of mechanical features observed in experiments performed on cardiac muscle samples or single cells, while being computationally simple. The model is specified by a system of ordinary differential equations for the kinetics of the interactions of the cardiomyocyte contractile and regulatory proteins and is similar, in this aspect, to some others model of myocardium mechanics [63][64][65]. However, the performance of our mechanical model is slightly better, and the set of experimental data reproduced by the model is larger than ones shown in the other model studies (see [25] for details). New modifications of the model described in Materials and Methods section were shown to fit some experimental data (force-calcium steady-state relationship) even better than the previous version of the model [25]. In addition, the modifications allowed us to improve computational efficiency of the model.

2.
In order to reproduce the dependence of the twitch tension on the stimulations frequency, one should couple the mechanical block with the model electrophysiology. This includes the description for the variation of the AP and the specifications of the electromechanical AP-Ca 2+ coupling. We describe the AP with the simple two-variable phenomenological model of Aliev-Panfilov [21] that reproduce the general form of the AP time-course and the AP dependence on the stimulating frequency well and was used and approbated in numerous publications. Our description of AP-Ca 2+ coupling is mostly based on the one introduced in the model by ten Tusscher and Panfilov [3]. However, we simplified the equations from [3] significantly making sure that the simplified equations were still relevant to the physical processes that the equations describe. An important feature of our model, which was absent in [3], is the introduction of calcium-dependent Ca 2+ uptake into sarcoplasmic reticulum via Ca-dependent phosphorylation of protein(s) involved into the uptake. While there are several models of heart electromechanics in the literature, some of which use the models [21] or [3] for the electrophysiology description and in some cases include detailed descriptions of myocardium mechanics [9][10][11], we were not able to find the results of reproducing the dependence of the twitch force and its relaxation rate on stimulation frequency.
3. Another novel feature of our model is the mechano-electrical feedback provided by the strain-dependence of the cell membrane capacitance and cell conductivity taken from [19].
In numerical experiments reported here, the model variables did not depend on spatial coordinate. However, preliminary results of the simulation of 2D muscle contraction demonstrate significant changes in the excitation-contraction waves caused by the strain-dependent electrophysiology. These results will be published elsewhere.
To summarize, we have not used some brand new approaches while developing our model, with exception for a couple of modifications of our previous model equations and some equations from the model [3]. However, we made an effort to combine and modify the models developed earlier to obtain new results that reproduce the important features of myocardium contraction with the least computational cost possible.

Limitations
The main limitation of the model is the absence of a detailed description of the ion currents and their contribution to the AP. For this reason, the model cannot describe some important details of heart electrophysiology, and we could not validate patient-or species-specific parameters of the electrical part of the model. On the other hand, the Aliev-Panfilov model [21] used here reproduces many major phenomena of AP propagation in normal and diseased human heart including reentrant arrhythmias [22][23][24]. Also, only a few Ca 2+ flows and currents were included in the model while changes in concentrations of Na + , K + and other ions were neglected. Therefore, any factors affecting these neglected processes cannot be reproduced by the model. Complicated cooperative processes involved in CIRC [39,66] were not taken into account. Instead, a simple approach of ten Tusscher and Panfilov [3] was further simplified. We cannot exclude that our model is unable to reproduce changes in CICR at some heart diseases. The processes involved in the mechano-electrical and mechano-calcium feedbacks probably are not instantaneous. There is a delay between changes in sarcomere length and changes in the electrical and biochemical processes. The delay was not taken into account here. Nevertheless, the model reproduced the slow response to stretch of cardiac muscle that takes tens of seconds (Fig 3).

Conclusions
The model suggested here combines simplified and modified version of previously described models [3,21,25]. The simplifications allowed us to make model more effective computationally due to the reduction of the number of differential equations and their stiffness compared to predecessor models. The modifications included the description of the mechano-electrical and mechano-calcium feedbacks and the dependence of Ca 2+ uptake into SR on Ca 2+ -dependent phosphorylation of protein(s) involved into the uptake.
The model describes various mechanical experiments with cardiac muscle including loaddependent relaxation (Fig 7), steady-state contractions, nonstationary isometric twitches (Fig  2), both instantaneous and slow mechanical responses to muscle stretch (Figs 1 and 3), and the dependence of the twitch amplitude and duration on the interstimulus interval (Figs 4-6).
Because of the above-mentioned modification and simplifications of the model, the time step of the numerical integration by the explicit Euler method can be increased to 0.1-0.2 ms. As the Aliev-Panfilov [21] model and the mechanical model do not require small mesh, the model appears to be suitable for 3D simulation of heart electro-mechanics.
Supporting information S1 File. Additional details of the model. The rule of the choice of the solution of Eq (19) for Ca 2+ concentration in subspace, the 0D mechanical problem setup for muscle sample contraction, and the results of the convergence test for the numerical method are specified. (PDF)