Transmural cellular heterogeneity in myocardial electromechanics

Abstract

Myocardial heterogeneity is an attribute of the normal heart. We have developed integrative models of cardiomyocytes from the subendocardial (ENDO) and subepicardial (EPI) ventricular regions that take into account experimental data on specific regional features of intracellular electromechanical coupling in the guinea pig heart. The models adequately simulate experimental data on the differences in the action potential and contraction between the ENDO and EPI cells. The modeling results predict that heterogeneity in the parameters of calcium handling and myofilament mechanics in isolated ENDO and EPI cardiomyocytes are essential to produce the differences in Ca²⁺ transients and contraction profiles via cooperative mechanisms of mechano-calcium-electric feedback and may further slightly modulate transmural differences in the electrical properties between the cells. Simulation results predict that ENDO cells have greater sensitivity to changes in the mechanical load than EPI cells. These data are important for understanding the behavior of cardiomyocytes in the intact heart.

Introduction

The term “heterogeneity” quite often signifies various features in cardiac physiology. In this paper we address electromechanical heterogeneity of cardiomyocytes from different regions of the ventricular wall. This heterogeneity includes both electrical and mechanical properties of the cells, specified as follows.

Distinctions in the intracellular expression of various ion channel isoforms ensure a pattern of cellular electrical heterogeneity. In particular, cardiomyocytes in the ventricular wall may have distinct expressions of the voltage-dependent Na⁺ channels [1], K⁺ channels [2], and L-type calcium channels [3] that result in significant differences in the morphology and duration of the action potential (AP) [4].

Distinctions in the intracellular properties of the cardiomyocytes responsible for their ability to shorten and generate the force, together assemble the cellular mechanical heterogeneity. The distinctions include different stiffnesses of cytoskeleton proteins [5], and various intracellular ratios of the fast (V1) and slow (V3) myosin isoforms [6, 7], revealing themselves in different kinetics of crossbridge attachment/detachment, etc. Differences in the Ca²⁺ sensitivity of myofilaments [8] also contribute to the mechanical heterogeneity, since Ca²⁺ activation of the crossbridges directly affects cardiomyocyte contractility.

All such specific regional features of the cellular electromechanical coupling in cardiomyocytes together contribute to distinctions in integrative functional characteristics of cardiomyocytes, such as AP morphology and duration [4], ‘force–velocity’ and ‘force–length’ relationships [5, 8, 9], and so on. Importantly, not only does electrical intracellular heterogeneity directly contribute to the AP distinctions between the cells, but mechanical heterogeneity also may contribute to the latter via cellular mechano-electric feedback mechanisms [10]. Analysis of mechanisms of cellular transmural electromechanical heterogeneity within the ventricular wall is a subject of the present work.

Regional heterogeneity in the electrophysiological and mechanical properties of myocardium in the normal heart has been observed in isolated cardiomyocytes [11, 12], wedge preparations [4, 13] and the whole heart [14, 15]. For example, AP duration is longer in cardiomyocytes from the subendocardial (ENDO) region than from the subepicardial (EPI) region of human ventricles [16], and different animal species such as dogs [12], guinea pigs [11], rats [17], mice [18], etc.

Several studies have reported a population of special M-cells found between the midmyocardial (MID) and EPI regions of the lateral wall, between the ENDO and MID regions of the anterior wall and in the deep regions of the papillary muscles, trabeculae and septa of human [13], and swine and canine ventricles [19, 20]. A distinctive feature of these M-cells is their significantly longer AP duration (APD) than is seen in EPI and ENDO cells at a low-frequency pacing rate [20]. Notably, the morphology of APs in cells from transmural regions also differs in various animal species. For instance, APs in EPI cardiomyocytes have a prominent initial repolarization phase with a deep notch (spike-and-dome configuration) in man, dog, rabbit and pig, but not in guinea pig, rat or mouse [11, 17, 18]. This configuration is less prominent in MID cells and is absent in ENDO cells.

Numerous experimental data demonstrate that significant differences in ion channel expression are responsible for the heterogeneity observed in AP morphology and duration [1, 21]. The transmural gradient in APD also correlates and may be associated with inherent transmural differences in calcium handling [3, 12]. Despite the large amount of experimental evidence for the presence of transmural heterogeneity in the electrophysiological properties of cardiomyocytes, much less is known about the heterogeneity in the mechanical function of these cells and the underlying mechanisms of excitation–contraction coupling. A few studies have focused on transmural variations in the passive and active mechanical properties of isolated myocytes, and regional differences in these properties have been reported [5, 22, 23].

Mathematical models may serve as a tool for identifying mechanisms underlying the regional heterogeneity of the myocardium in normal and pathological conditions. Several mathematical models take into account the transmural gradient in the electrical properties of the myocardium at different levels of myocardial tissue organization [24,25,26,27,28]. However, the majority of these models are based and focused on the electrophysiology only, and either do not describe the mechanical function at all or take into account regional differences in the electrical function and calcium handling without considering the transmural gradient in the mechanical activity of the cardiomyocytes. In contrast to the electrophysiology, rather few modeling studies have addressed the transmural gradients in the cellular mechanical activity [29] or have analyzed the effects of the intra- and intercellular mechano-electric feedback on the myocardium function [30,31,32,33].

As we know, only one team of authors developing regional cellular models paid attention to the transmural distinction in intracellular mechano-electric feedback mechanisms, namely that in the crossbridge kinetics [29]. However, even their model did not account for some existing experimental data, e.g., a smaller Ca²⁺ release from the sarcoplasmic reticulum [3], a greater myofilament Ca²⁺ sensitivity [8], and a higher stiffness [5, 8, 9] in the ENDO cells relative to that in the EPI cells. So, in this work we were trying to assess the contribution of these mechanisms to cellular heterogeneity. Moreover, in contrast to other models simulating either unloaded cell shortening [29], or isometric twitches at a constant sarcomere length [28], we have studied effects of intracellular heterogeneity on the specific features of cellular contraction under more physiological conditions of imposed mechanical load, i.e., simulating afterloaded cellular contractions. One more specific aim of this study was to assemble in the models as much as possible available experimental data on the cellular heterogeneity specific to the guinea pig heart.

We previously published our first steps toward developing integrative mathematical models that describe the transmural features of the electromechanical coupling in cardiomyocytes from the ENDO and EPI regions of the guinea pig left ventricle (LV) [34, 35]. We have taken into account the regional features of potassium currents [34] and Na⁺/Ca²⁺ exchange currents [35] in ENDO and EPI cells. In our previous study [34] we fitted parameters to the current–voltage characteristics of rapid and slow components of delayed K⁺ current (i _Kr and i _Ks) registered by Bryant and co-workers in guinea pig cardiomyocytes [11] and used the currents in the ENDO/EPI models. This allowed us to reproduce shorter APD in the EPI cells than in the ENDO cells [34]. We have also tested the hypothesis that the transmural gradient in i _NaK may underlie the heterogeneity in i _NaCa in isolated ENDO and EPI myocytes [36]. Following the experimental findings, we suggested a higher i _NaK amplitude in the ENDO model and showed a resulting smaller density of i _NaCa in the reverse mode of the exchanger in the ENDO model, which leads to an increase in the differences in APD between the ENDO and EPI models [35]. These two ionic mechanisms, which may underlie the differences in electrical properties of the ENDO and EPI cells, were utilized in the cellular models presented here.

However, when only the above mechanisms were used in our models, the models failed to reproduce experimental data on significant transmural differences in both electrical and mechanical function of ENDO and EPI cardiomyocytes registered experimentally [23, 37]. We hypothesized that this mismatch might be resolved by accounting in the models for differences in calcium handling [3, 13, 38] and intracellular mechanical properties inherent to the cells from different myocardial regions [5, 9]. In this paper, we propose improved ENDO and EPI cellular models with an extended set of different parameters, which reflect experimental findings on the transmural mechanical heterogeneity between the cells (Table 1).

Table 1 Experimental data on the transmural differences in intracellular electrical and mechanical properties between the ENDO and EPI cells taken into account in the cell subtype models

Our study is aimed at the development of the integrative electromechanical models of the ENDO and EPI cardiomyocytes taking into account mechano-calcium-electric feedbacks, and based on the species specific experimental data available at the moment. The models were used to assess contribution of individual cellular mechanisms to the overall transmural gradient in the cellular electrical and mechanical performance.

Materials and methods

EO model of cellular electromechanics as a predecessor of ENDO and EPI models

Our ENDO and EPI models of the electrical and mechanical activity in cardiomyocytes from ENDO and EPI regions are based on an Ekaterinburg–Oxford mathematical model (EO model), which integrates the ionic and contractile mechanisms of excitation–contraction coupling in the guinea pig ventricular cardiomyocyte [39].

The main components of the original EO model are described below. Figure 1 shows schematic diagrams of the EO model for electrophysiological, calcium handling and mechanical blocks of the model.

Fig. 1

Scheme of a single cell in the EO model. a Scheme of the ionic dynamics for a single-cell model. Simulated ionic currents (i _x ) involved in the AP generation are shown: i _Na fast Na⁺ current, i _pNa persistent Na⁺ current, i _bNa background Na⁺ current, i _NaCa Na⁺/Ca²⁺ exchange (NCX) current, i _CaL Ca²⁺ current via L-type Ca²⁺ channels, i _bCa background Ca²⁺ current, i _to transient outward K⁺ current, i _Kr , i _Ks rapid and slow components of delayed K⁺ current, i _K1 rectifier K⁺ current, i _NaK Na⁺/K⁺ pump current, DS dyadic space, TC terminal cisterns, LR longitudinal reticulum, i _rel Ca²⁺ release flow from the sarcoplasmic reticulum (SR), i _up SERCA pump flow from cytosol to the SR, i _tr Ca²⁺ diffusion between the LR and the TC, CaS Ca²⁺ complexes with calsequestrin, CaTnC Ca²⁺ complexes with TnC, B ₁ , B ₂ Ca²⁺ complexes with fast and slow buffer ligands. b Rheological scheme for a single cell model, where a contractile element (CE) is connected to an in-series and parallel passive elastic element (SE, XSE, PE), and a viscous element (VS) is inserted in parallel to the PE. Variables l, l ₁ and l _ex define deformations of PE, CE and XSE, respectively, relative to their slack lengths (see text for details)

An electrophysiological block of the EO model is based on the Noble′98 ionic model [40]. The transmembrane potential is calculated with the following equation:

$$\begin{aligned} {\frac{dV}{dt} = \frac{ - 1}{{C_{\text{m}} }}\left( {i_{\text{stim}} + i_{\text{Na}} + i_{\text{pNa}} + i_{\text{bNa}} + i_{\text{NaCa}} + i_{\text{CaL}} + i_{\text{bCa}} + i_{\text{to}} + i_{\text{Kr}} + i_{\text{Ks}} + i_{\text{K1}} + i_{\text{NaK}} } \right.),} \hfill \\ i_{\text{stim}} = - 3.0{\text{ nA }}\left( {{\text{every cardiac cycle on at 6}}0{\text{ ms}};{\text{ off at 62}}. 5 {\text{ ms}}} \right). \hfill \\ \end{aligned}$$

(1)

where C _m is the membrane capacitance, i _stim is the stimulation current, i _Na is the fast Na⁺ current, i _pNa is the persistent Na⁺ current, i _bNa is the background Na⁺ current, i _NaCa Na⁺/Ca²⁺ is the exchange current, i _CaL Ca²⁺ is the current via L-type Ca²⁺ channels, i _bCa is the background Ca²⁺ current, i _to is the transient outward K⁺ current, i _Kr, i _Ks are the rapid and slow components of the delayed K⁺ current, i _K1 is the rectifier K⁺ current, and i _NaK Na⁺/K⁺ is the pump current (Fig. 1a).

Figure 2 shows a scheme of links between mechanisms of the excitation contraction coupling and mechano-calcium-electric feedback taken into account in the EO model. AP generation governs cytosolic Ca²⁺ transient (direct links, solid lines) and in turn depends on Ca²⁺ kinetics via transmembrane Ca²⁺ fluxes (feedback links, dashed lines).

Fig. 2

Schematic links between mechanisms of electromechanical coupling and mechano-calcium-electric feedback in the EO model. Solid lines show direct links between the mechanisms of excitation–contraction coupling and dashed lines show feedback links. Cooperative mechanisms of calcium activation of myofilament RUs (Xb–CaTnC, CaTnC–CaTnC and RU end-to-end cooperativity) are described in the text

The equation for the concentration of free cytosolic Ca²⁺ ([Ca²⁺] _i ) describes Ca²⁺ fluxes between the cytosol and the dyadic space (DS), longitudinal reticulum (LR) and terminal cisterns (TC) of the sarcoplasmic reticulum (SR) and takes into account Ca²⁺ binding to myofilament regulatory protein troponin C (TnC) and other Ca²⁺ binding ligands (Fig. 1a) [39]:

$$\begin{aligned} \frac{{d[{\text{Ca}}^{2 + } ]_{i} }}{dt} & = \frac{ - 1}{{2 \cdot V_{\text{Cell}} \cdot V_{{i\_{\text{ratio}}}} \cdot F}} \cdot \left( {i_{{{\text{CaLCa}}_{\text{cyt}} }} + i_{\text{bCa}} - 2 \cdot i_{{{\text{NaCa}}_{\text{cyt}} }} } \right) + [{\text{Ca}}^{2 + } ]_{\text{ds}} \cdot 10 \cdot V_{{{\text{ds}}\_{\text{ratio}}}} \\ & \quad + \frac{{V_{{{\text{rel}}\_{\text{ratio}}}} }}{{V_{{i\_{\text{ratio}}}} }} \cdot i_{\text{rel}} - \frac{{d[{\text{CaTnC}}]}}{dt} - \frac{{d[B_{1} ]}}{dt} - \frac{{d[B_{2} ]}}{dt} - i_{\text{up}} , \\ \end{aligned}$$

(2)

where V _Cell is the cell volume, V _{i_ratio} is the volume fraction of cytosolic space, V _{ds_ratio} is the volume ratio of DS, V _{rel_ratio} is the volume ratio of the SR release site, F is Faraday’s constant, i _CaLCacyt is the Ca²⁺ current via L-type Ca²⁺ channels, i _bCa is the background Ca²⁺ current, i _NaCacyt is the Na⁺/Ca²⁺ exchange current, [Ca²⁺]_ds is the Ca²⁺ concentration within the DS, i _rel is the Ca²⁺ release flux from the SR, [CaTnC] is the concentration of Ca²⁺ complexes with TnC, [B ₁], [B ₂] are the concentrations of Ca²⁺ complexes with fast and slow binding ligands, and i _up is the sarcoplasmic reticulum (SR) Ca²⁺ ATPase (SERCA) pump flux from the cytosol to the SR (Fig. 1a).

Note that [Ca²⁺]_ds ˑ10ˑV _{ds_ratio} in Eq. (2) is a term used in the model Noble′98 [40] to describe Ca²⁺ diffusion from DS to cytosol. Strictly speaking, the formula D ([Ca²⁺]_ds − [Ca²⁺] _i ) would be a more correct description of the diffusion between DS and the cytosol (with diffusion coefficient D) in this equation, allowing for a reverse flow if Ca²⁺ concentration is higher in the cytosol than in the DS. Nevertheless, (as we checked many times) this situation actually does not occur during the excitation–contraction cycle in the model Noble′98, probably owing to the passing of the input calcium flow into the cytosol only through the DS. Therefore, “unidirectional diffusion” of Ca²⁺ from the DS to the cytosol in Eq. (2) seems to be a permissible simplification inherent from the model Noble′98, so we left it unchanged when we incorporated it into the EO model.

The cooperative mechanisms of calcium activation of myofilament regulatory units (RUs) accounted for in the model and respective formulations are identified and justified in our previous papers [39, 41]. Here, we describe them briefly (see Fig. 2 for the scheme):

Xb–CaTnC cooperativity The off-rate of CaTnC dissociation decreases with an increase in the fraction of force-generating cross-bridges (Xbs) per single CaTnC complex.

CaTnC–CaTnC cooperativity The off-rate of CaTnC complex dissociation decreases with increasing numbers of CaTnC complexes in proximity.

RU end-to-end cooperativity Ca²⁺ binding by TnC, located within one regulatory unit (RU) of actin–tropomyosin–troponin on a thin filament affects the neighboring RUs through the tropomyosin end-to-end conformational interaction, thus contributing to the opening of the active actin sites for attachment of myosin heads.

The equation for the CaTnC kinetics governs Xb kinetics [see Eq. (4) below for the fraction N of force-generating Xbs] and cooperatively depends on this kinetics (note nonlinear dependence of the off-rate on N via Xb–CaTnC cooperative mechanism), making Ca²⁺ kinetics mechano-sensitive as well [39]:

$$\frac{{d[{\text{CaTnC}}]}}{dt} = k_{\text{on}} \left( {[{\text{TnC}}]_{\text{tot}} - [{\text{CaTnC}}]} \right) \cdot [{\text{Ca}}^{2 + } ]_{i} - k_{\text{off}} ([{\text{CaTnC}}],N) \cdot [{\text{CaTnC}}],$$

(3)

where $k_{\text{on}}$ and $k_{\text{off}} ([{\text{CaTnC}}],N)$ are on- and off-rate “constants” of CaTnC complex formation. Note nonlinear dependence of the off-rate on [CaTnC], accounting for CaTnC–CaTnC cooperativity in the RUs.

Equations describing calcium handling in cardiomyocytes (Eqs. 2, 3) play a crucial role in simulating both the electromechanical and mechano-electrical coupling in the EO model (see Fig. 2). As the direct link, the time course of CaTnC complexes via cooperative mechanisms of myofilament RU controls kinetics of force-generating Xbs [see Eq. (4) below] and thus, the mechanical behavior of the cardiomyocyte (Fig. 2, solid lines). As the feedback link, the mechanical activity of the sarcomere via the Xb–CaTnC cooperative mechanism affects the kinetics of the CaTnC complexes, which together with the CaTnC–CaTnC cooperative mechanism (Eq. 3) affects kinetics of free intracellular calcium, which can provide the feedback on calcium-dependent ionic currents and AP generation (Fig. 2, dashed lines). These Ca-RU cooperative mechanisms underlie the mechano-calcium-electric feedback in the model.

In the rheological scheme of the model (Fig. 1b), the active contractile element CE is associated with the sarcomere ensemble of the cardiomyocyte. Sarcomeres generate cellular force and shortening during contractions due to active cycling of force-producing Xbs between actin and myosin, which are allowed by the conformations of the myofilament RUs triggered by the binding of activating Ca²⁺ ions by troponin C (Fig. 2, see also [39, 41]). Active CE interacts with the elastic and viscous elements (SE, PE and VS), which mainly determine the mechanical properties of the passive myocardium, but may also modulate the active myocardial mechanics [39].

The equation for the fraction N of force-generating Xbs directly governs the mechanical behavior of the contractile element, but it also depends itself on the sarcomere mechanics [sarcomere length (SL) l ₁ and shortening velocity v] [39]:

$$\frac{dN}{dt} = k_{\text{on}} ([{\text{CaTnC}}],l_{1} ,v) \cdot (1 - N) - k_{\text{off}} (v) \cdot N,$$

(4)

where $k_{\text{on}} ([{\text{CaTnC}}],l_{1} ,v)$ and $k_{\text{off}} (v) = k_{{{\text{m}}_{\text{v}} }}$ are on- and off-rate “constants” of force-generating Xb cycling that depend on the SL and the shortening velocity. Note the nonlinear dependence of the on-rate for Xb attachment on [CaTnC], accounting for direct effects of [CaTnC] and for the RU end-to-end cooperativity (Fig. 2):

$$M([{\text{Ca}}_{\text{TnC}} ]) = \frac{{\left( {\frac{{[{\text{CaTnC}}]}}{{[{\text{TnC}}]_{\text{tot}} }}} \right)^{\mu } \cdot \left( {1 + 0. 6^{\mu } } \right)}}{{\left( {\frac{{[{\text{CaTnC}}]}}{{[{\text{TnC}}]_{\text{tot}} }}} \right)^{\mu } + \,0. 6^{\mu } }},$$

(5)

where [CaTnC] is the concentration of CaTnC complexes, [CaTnC]_tot is the total TnC concentration, and µ is the parameter of cooperativity.

Therefore, nonlinear dependence of the on-rate for Xb attachment on [CaTnC] (via RU end-to-end cooperativity) and the off-rate for CaTnC dissociation on the Xb (via the Xb–CaTnC cooperative mechanism) and CaTnC (via the CaTnC–CaTnC cooperative mechanism) formalize key cooperative mechanisms of myofilament Ca²⁺ activation in the model. These mechanisms underlie a number of essential properties of cardiac muscle, such as the dependence of contraction and relaxation on myofilament length and load.

Suppose that l ₁ is the relative change in the SL of the cell against its slack length (deformations of the CE, normalized by the slack length) and l is a relative change in the cell length per sarcomere (deformations of PE normalized by the sarcomere slack length). The myofilament force generated by the CE (F _CE) is defined by the fraction of force-generating Xbs N and depends on the velocity of the CE shortening/stretching ($v = \frac{{dl_{1} }}{dt}$) (Fig. 2) [39]:

$$F_{\text{CE}} = \lambda \cdot N \cdot p(v),$$

(6)

where $\lambda$ is the model parameter. Function $p(v)$ is the dependence of the average Xb force on the sarcomere shortening/lengthening velocity (for more details, see the “Appendix”).

Passive force of the cell in the model as the force of the PE element is calculated with the following nonlinear stress–strain dependence [39]:

$$F_{\text{PE}} = \beta_{2} \cdot \left( {e^{{\alpha_{2} \cdot l}} - 1} \right),$$

(7)

where β ₂, α ₂ are model parameters (for more details, see the “Appendix”).

The elastic element SE and the damper VS also modulate the passive and active mechanical properties [39]. The element SE generates a force proportional to the difference between l and l ₁. The damper VS, being parallel to PE, generates a force proportional to the velocity of the cell shortening (see the “Appendix”).

The CellML model representation of the original EO model [39] can be found at http://models.cellml.org/e/b9/ and run using the Cellular Open Resource at http://cor.physiol.ox.ac.uk/ [42].

In this study, we have fitted parameters of the EO model to simulate specific features of ENDO and EPI cells registered experimentally (Table 1). In the next Section, “ENDO and EPI cellular model development”, we describe selection of model parameters, which differ in the ENDO and EPI models. The complete set of equations and parameters for each cell type model are provided in the “Appendix”.

ENDO and EPI cellular model development

We developed our electromechanical ENDO and EPI models to account for experimental data on transmural differences observed in intracellular electrical and mechanical properties between cells from different regions of the LV (Table 1). Justification of the model parameters that differ in the ENDO and EPI models is given in detail, as follows.

Heterogeneity in the electrophysiological properties

The late (persistent) sodium current, i _pNa

Osadchii and co-workers investigated the protein expression distribution of voltage-dependent Na⁺ channels (Na_v) across the ventricular wall in the guinea pig heart. They demonstrated higher Na_v protein expression in the subendocardium than in the subepicardium in both heart ventricles and concluded that this may contribute to greater excitability and higher susceptibility to repolarization alternans in the endocardium [1]. Zygmunt and co-workers showed a somewhat larger i _pNa density in the subendocardium than in the subepicardium in canine hearts [43]. Based on these data, the maximum membrane i _pNa conductance was assigned a higher value in the ENDO model than in the EPI model to simulate experimental data on the transmural differences in the amplitudes of i _pNa (≈13%) at 0 mV [43] (see the Appendix “Persistent Na⁺ current, i _pNa ”, Appendix Table 6).

Fast and slow components of delayed potassium current, i _Kr and i _Ks

Bryant and co-workers demonstrated that the differences in the APD in the ENDO and EPI cells of the guinea pig can be determined by differences in the density of i _Kr and i _Ks. They found that the density of i _Kr and i _Ks is higher in the EPI than in the ENDO cardiomyocytes that is consistent with shorter AP of the EPI cells [11]. In our recent study, we took into account the heterogeneity in potassium currents in the guinea pig, reproducing the current–voltage characteristics of i _Kr and i _Ks as registered experimentally in the EPI and ENDO cells [34]. Here, we used these findings in the cellular models (see Appendix “Rapid delayed rectifier potassium current, i _Kr ”, “Slow delayed rectifier potassium current, i _Ks ”, Appendix Table 6).

L-type calcium current, i _CaL

No differences were found in the current–voltage characteristics of i _CaL between ENDO and EPI cells in the mouse [18] and dog [45]. Nonetheless, the expression of the L-type calcium channel protein DHP_α1 was greater in the subendocardium than in the subepicardium in the guinea pig LV [3]. According to these data, the membrane permeability for Ca²⁺ through L-type channels was set as larger in the ENDO model than in the EPI model (see Appendix “L-type Ca²⁺ current, i _CaL ”, Appendix Table 6). The slightly higher maximum membrane i _bCa conductance chosen for the ENDO cell maintains a greater Ca²⁺ current into the cell (see Appendix “Background Ca²⁺ current, i _bCa ”, Appendix Table 6).

Na⁺/Ca²⁺ exchange and Na⁺/K⁺ pump currents, i _NaCa and i _NaK

Experimental data on the regional differences in i _NaCa and i _NaK in LV are still controversial. Laurita and co-workers found no significant differences in the expression of the dominant isoform of NCX₁ in the dog [38], while other authors registered i _NaCa as smallest in ENDO cells in canine and mouse hearts [18, 44]. Quinn and co-workers have demonstrated a higher i _NaCa density and NCX protein expression in ENDO cells in rabbit [21]. Gao and co-workers found no differences in NCX protein expression within the canine LV wall, but revealed a transmural gradient in the corresponding Na⁺ concentration [Na⁺] _i [36] and suggested that this transmural gradient in i _NaCa may be generated by a transmural gradient in the expression of Na⁺/K⁺ ATPase.

In our recent study, we tested the hypothesis that the transmural gradient in i _NaK underlies the heterogeneity in i _NaCa in isolated ENDO and EPI myocytes. We assigned a higher i _NaK amplitude in the ENDO model and showed a resulting smaller density of i _NaCa in the reverse mode of NCX activity in the ENDO model, which leads to an increase in the differences in APD between the subtype models [35] (see Appendix “Na⁺–K⁺ pump, i _NaK ”, Appendix Table 6). This mechanism is utilized in the cellular models presented here.

Heterogeneity in Ca²⁺ handling

Based on the experimental data showing lower expression of SERCA2 in ENDO cells than in EPI cells of the LV [3, 13, 38], the rate constant of the SERCA Ca²⁺ uptake was given a smaller value in the ENDO model than in the EPI model (see Appendix “Ca²⁺ uptake by the SR-pump, i _up ”, Appendix Table 7). Wan and co-workers suggested that less Ca²⁺ was released from the SR in the ENDO region of the guinea pig [3]. We achieved this by suggesting a lower sensitivity of the binding sites of the ryanodine receptors to cytosolic Ca²⁺ in the ENDO model (see Appendix “Ca²⁺ release from the SR, i _rel ”, Appendix Table 7).

Heterogeneity in the mechanical properties

The ENDO cells in guinea pig and canine hearts demonstrate slower time to peak shortening and a larger relaxation time constant, when compared to EPI cells [12, 23, 37]. Transmural differences in the active mechanical properties of cardiomyocytes may be associated with regional differences in Xb cycling kinetics. Ait Mou and co-workers demonstrated in skinned cardiomyocytes of guinea pig LV at saturated Ca²⁺ concentration that maximal rate constant of force redevelopment K _tr after release is higher, and restretch to the initial cell length is more rapid in EPI than in ENDO cells (Fig. 3), suggesting a higher rate constant of Xb cycling in EPI cells [8]. Similar data were also obtained in skinned multicellular preparations of myocardium from different regions of the LV in the pig [7].

Fig. 3

Ratio of maximal rate constant of force redevelopment K _tr in EPI to ENDO cells registered in experiments [8] and produced by the models at SL₀ = 2.3 μm. K _tr was determined by monoexponential fitting to the maximal tension for each pCa (= −log[Ca²⁺]). Ratio of mean values of K _tr from experimental data [8] is used. The latter data we calculated as a ratio of the K _tr values corresponding to equal values of pCa [see left and right plots in Fig. 3b (for SL₀ = 2.3 μm) in the cited article]

To take into account these experimental findings and to achieve higher velocity of contraction-relaxation in the EPI cells than in ENDO cells, we assigned a greater maximal velocity of unloaded sarcomere shortening (see Appendix “Kinetics of force generating Xbs, N”, “Mechanical variables, F _CE, F _VS, F _SE, F _PE, F _XSE, l, l ₁, l _ex ”, Appendix Table 8) and a higher rate constant of Xb cycling in the EPI model (see Appendix “Mechanical variables, F _CE, F _VS, F _SE, F _PE, F _XSE, l, l ₁, l _ex ”, Appendix Table 8). These modulations in model parameters allowed the models to reproduce qualitatively experimental data on higher K _tr in EPI cardiomyocytes (Fig. 3). Note that quantitative mismatch between experimental data and simulations may be attributed to the differences in parameters of skinned cardiomyocytes used in experiments and that of simulated intact cardiomyocytes.

Ait Mou and co-workers showed in skinned cells that Ca²⁺ sensitivity of force at SL of about 2.3 μm was higher in ENDO cells than in EPI cells, suggesting a greater molecular cooperativity of myofilament Ca²⁺ activation in the ENDO cells [8]. Based on these data, we set a different degree of nonlinearity μ for the RU end-to-end cooperativity in the cells (Eq. 5, see also Appendix “Kinetics of force generating Xbs, N ”, Appendix Table 8). As we have shown recently, a lower parameter µ in Eq. 5 for the ENDO model makes the effects of RU end-to-end interaction on the on-rate of cross-bridge attachment higher in the ENDO model than in the EPI model [46]. This allowed us to reproduce a steeper Ca²⁺-force relationship in the ENDO cells at the physiological range of cytosolic Ca²⁺ concentrations as shown experimentally. The modeling results suggest that such a difference in myofilament cooperativity may also contribute to the regional gradient in the active tension development between the cells [46].

Transmural differences in the passive mechanical properties of cardiomyocytes have been observed in several animal species [5, 9, 22]. The ‘SL-passive tension’ curve is much steeper in ENDO cells than in EPI cells, specifying their greater stiffness [5, 8, 9]. These authors suggested that differences in the passive characteristics of cardiomyocytes from the ENDO and EPI regions might be associated with transmural differences in the distribution of titin isoforms [6, 9]. Titin is a large myofilament protein that extends from the Z- to M-line of the sarcomere, and is thought to be the major determinant of passive mechanical properties in cardiomyocytes [22]. In our models, we have reproduced experimental data on the ‘SL-passive tension’ relation registered in isolated ENDO and EPI cardiomyocytes from the guinea pig LV [9] (Fig. 4 and see Appendix “Mechanical variables, F _CE, F _VS, F _SE, F _PE, F _XSE, l, l ₁, l _ex ”, Appendix Table 8).

Fig. 4

‘SL-passive tension’ curves. The ‘SL-passive tension’ curve in ENDO and EPI cells derived from experimental data (dashed-dotted linear regression for each group of cells from [9]) against the model fitting (solid lines)

In silico experimental protocols

External ionic concentrations in the ENDO and EPI models are assumed to be constant. In single cell experiments, the cells were prestretched to have a slack sarcomere length SL_o of 2.0 µm. For simulations we used a Euler integration time step of 0.01 ms. All the results are shown in their steady state, achieved by allowing the models to run for 100 s at a 1-Hz pacing rate. Isometric twitches and afterloaded isotonic contractions were examined.

Results

Transmural differences in action potential characteristics

At first, we estimated individual contributions of regional differences in the model parameters of ionic currents to the differences in APD, and the time course of calcium transient and contraction between the ENDO and EPI cells (Table 2). In these simulations, parameters of the ENDO model were all the same as the parameters of the EPI model except the only ionic mechanism under examination which had its parameters set the same as in the complete ENDO model (see the Appendix Table 6, ENDO parameters).

Table 2 Contribution of ENDO/EPI differences in individual ionic currents to differences in action potential, Ca²⁺ transient, and contraction between the ENDO/EPI cellular models

Figure 5 shows the time course of AP in the complete ENDO and EPI models and superposition of the ionic currents that primarily determine the repolarization differences between ENDO and EPI cells. Table 2 and Fig. 5 show that higher density and amplitude of depolarizing i _CaL in the ENDO model most significantly contributes to the prolonged plateau phase and to the delayed initial rapid repolarization of the AP. Although the difference in i _NaK between the models is viewed as negligible (Fig. 5), this difference led to the regional heterogeneity in i _NaCa, as mentioned above (see the Na⁺/K⁺ pump and Na⁺/Ca²⁺ exchanger section). The modeling results suggested that the smaller outward i _NaCa in the ENDO model contributes essentially to the delayed initial rapid repolarization and more prolonged plateau phase in the AP (Table 2; Fig. 5). Our simulation data are consistent with the experimental hypothesis that also the higher outward i _Kr and i _Ks and the smaller inward i _pNa density in the EPI model contribute to the discrepancy in APD observed between the cells [11, 43], providing for the faster repolarization phase in the EPI AP (Table 2; Fig. 5). The simulation results showed that the gradients in the ionic currents between the ENDO and EPI cellular models had no effect on the gradient in the resting membrane potential (RP) (Table 2), in accordance with experimental recordings [11, 37].

Fig. 5

Simulation of regional differences in the electrical properties between ENDO and EPI cells. Action potential (V) in the ENDO and EPI models during an isometric twitch (top panel) and the ionic currents (i _ion) mainly underlying the discrepancy in APs between the models (bottom panel). The stimulus was applied at 60 ms

To reveal the contribution of the heterogeneity in calcium handling and mechanical properties to the transmural gradient in electromechanical activity of the cells, we compared the EPI model with three “ENDO” models accounting for different cellular mechanisms of cellular heterogeneity: (a) ENDO (E) model that differed from the EPI model in the parameters of ionic currents only (see Appendix Table 6), (b) ENDO (E + Ca) model that accounted for the heterogeneity in the electrical and calcium handling properties between the cells (see Appendix Tables 6, 7), and (c) the complete ENDO model accounting for the heterogeneity in the electrical, calcium handling and mechanical properties (see Appendix Tables 6, 7, 8). First, it is seen that the specific features of ionic currents in the ENDO (E) cell produced essentially longer AP in the ENDO cell, which had no effect on Ca²⁺ transient but ensured slower contraction (Table 3). However, the gradient in characteristics of contraction between the cells was not as pronounced as registered experimentally (Table 4). Heterogeneity in calcium properties [ENDO (E + Ca) model, Table 3] substantially prolonged the 90% decay phase of Ca²⁺ transient (CaT₉₀) due to slower calcium uptake from the cytosol, and essentially delayed T _max and TR₇₅ in the ENDO cell (Table 3). Finally, the change in the mechanical parameters for the complete ENDO model produced the most prolonged APD₉₀, CaT₉₀, and T _max (Table 3). As compared to the ENDO (E + Ca) model, delayed contraction in the complete ENDO model essentially prolonged the calcium transient and, via mechano-calcium-electric feedback, slightly prolonged further the AP (Table 3). In additional simulations, we tested the ENDO model without Xb–CaTnC and CaTnC–CaTnC cooperativity, excluding mechano-calcium-electric feedback in our models [ENDO (–MEF), Table 3]. In this case, delayed T _max and TR₇₅ had no effect on CaT₉₀ and AP as compared to the ENDO (E + Ca) model. Moreover, APD in this model proved to be even smaller than APD produced by either ENDO (E + Ca) or ENDO models. Thus, we showed quantitatively that the gradient in calcium handling parameters is most essential for the distinctions in Ca²⁺ transient and contraction between the cells, while different mechanical parameters further enhance the transmural gradient in mechanical function and slightly influence the differences in electrical function between ENDO and EPI cells.

Table 3 Contribution of ENDO/EPI differences in the electrical, calcium handling and mechanical properties to the characteristics of action potential, Ca²⁺ transient, and contraction at isometric twitches

Table 4 Characteristics of action potential and low-loaded shortening in the ENDO and EPI models, compared to experimental data [23]

Note that all simulations presented below were produced by the complete ENDO model and compared to the EPI model.

Comparison of AP duration in the developed ENDO and EPI models to experimental data registered in the ENDO and EPI myocytes isolated from the guinea pig LV [23, 37] showed good consistency (Table 4).

Ca²⁺ handling in ENDO and EPI single cell models

When compared to the EPI model, the ENDO model exhibited a slightly longer time to peak and prolonged 90% decay phase of the Ca²⁺ transient (Table 3) that qualitatively corresponded to experimental recordings in guinea pig ENDO and EPI cardiomyocytes (~10% difference in Ca²⁺ transient duration at 90% decay (CaT₉₀) [3]) and canine transmural wedges (~20% in CaT₉₀ [3, 13, 38]).

In accordance with the experimental data [3], the diastolic [Ca²⁺] _i was similar in the ENDO and EPI models (9.18 and 9.48 nM under isometric conditions, respectively). The lower maximal velocity of the SR pump led to a lower SR Ca²⁺ content in the ENDO model than in the EPI model (1.37 and 1.49 mM under isometric conditions, respectively). These model predictions are consistent with a higher SR content reported in isolated EPI cardiomyocytes in canines, as estimated by caffeine application for a period of 1 s under voltage-clamp conditions [12].

Our ENDO and EPI models produced almost the same peaks of [Ca²⁺] _i (1.78 vs 1.68 μM, respectively). Experimental data on the differences in peak [Ca²⁺] _i are controversial. Wan and co-workers demonstrated a lower amplitude of the Ca²⁺ transient in the ENDO cells of guinea pig LV than in the EPI cells, which is consistent with lower expression of SERCA2 in the subendocardium and a resulting lower Ca²⁺ content in the SR [3]. More data were observed in cardiomyocytes of rabbit, dog and man, where no significant differences were found in the peak [Ca²⁺] _i between the cell subtypes, despite lower expression of SERCA2 in the subendocardium [13, 47, 48].

We determined the mechanisms underlying the differences in Ca²⁺ transients between ENDO and EPI cells by analyzing the amount of Ca²⁺ (∆[Ca²⁺] _i ) that enters into and exits from the cytosol over the duration of the calcium transient (≈500 ms) during isometric twitches (Fig. 6).

Fig. 6

Amount of Ca²⁺ (∆[Ca²⁺] _i , µM) that enters into and extrudes from the cytosol via ionic channels and exchangers in ENDO and EPI cells over the duration of the calcium transient (≈500 ms) during an isometric twitch. Numbers show percentage of the values in the ENDO model with respect to the EPI model. Integrals of the following simulated ionic currents and fluxes (i _x ) carrying Ca²⁺ are shown: i _CaL Ca²⁺ current via L-type Ca²⁺ channels, i _bCa background Ca²⁺ current, i _NaCa(f) NCX current during the forward mode, i _NaCa(r) NCX current during the reverse mode, i _rel Ca²⁺ release flow from the SR, i _up SERCA pump flow from cytosol to the SR

A greater amount of Ca²⁺ (+34.4%, 9.4 µM in the ENDO model vs 7.0 µM in the EPI model) entered the ENDO cell via i _CaL, while a lesser amount of Ca²⁺ entered the cytosol via i _NaCa(r) (−34.9%, 1.7 vs 2.6 µM) and a lesser amount of Ca²⁺ extruded from the cytosol via i _NaCa(f) (+2.2%, 17.6 vs 18.0 µM) in the ENDO cell. Thus, in total, about the same amount of Ca²⁺ entered the ENDO and EPI cells via transmembrane currents over the duration of calcium transient. Despite a smaller amplitude of i _rel in the ENDO cell, the Ca²⁺ release from the SR was longer lasting, resulting in the release of a greater amount of Ca²⁺ from the SR (+8.6%, 84.7 vs 78.0 µM). Therefore, the latter mechanism provided for almost identical Ca²⁺ amplitude, despite the lower SR Ca²⁺ content, and contributed to the slower calcium transient in the ENDO cell compared to the EPI cell (Table 3).

Transmural differences in the cellular mechanical function

Figure 7 shows the concentration of CaTnC, the fraction of force-generating Xbs, and the isometric tension generated by the ENDO and EPI models. In the ENDO model, the longer AP and slower Ca²⁺ transient contributed to a delayed time to peak contraction (T _max) and a slower relaxation of contraction (longer period of time from peak to 75% of relaxation, TR₇₅) when compared to the EPI model. Model simulations predicted that either of the ionic mechanisms individually (see Table 2) or their combination [see Table 3, ENDO (E)] might contribute essentially to the delayed contraction and relaxation in the ENDO cells. The most prominent effect was predicted from enhanced i _CaL in the ENDO model, which alone provided for about 10% discrepancy between T _max and about 20% between TR₇₅ in the cells. Transmural differences in the inherent parameters of intracellular Ca²⁺ handling essentially increased discrepancy between T _max and TR₇₅ during contraction [compare ENDO (E) to ENDO (E + Ca) in Table 3]. Finally, the mechanical properties also directly contributed to the heterogeneity observed in the Xb kinetics between the cells [see ENDO (E + Ca + M) in Table 3; Fig. 7], which, via mechano-dependence of Ca²⁺ kinetics, modulate the Ca²⁺ transient and Ca²⁺ currents and, consequently, both AP and contraction profile in the cells. In the EPI model, the higher rate constant of Xb cycling contributed to the faster contraction of the cell (Fig. 7). In the ENDO model, the higher RU end-to-end cooperativity additionally contributed to slower attachment-detachment of the Xbs, which further slowed down contraction in the ENDO model versus the EPI model (Fig. 7).

Fig. 7

Cellular mechanics. Time dependent concentration of CaTnC ([CaTnC]) (top panel), the fraction of force-generating Xbs (N, dotted lines), and the isometric tension (F, solid lines, bottom panel) in the ENDO and EPI models. The stimulus was applied at 60 ms

Our EPI and ENDO models produced similar contraction amplitudes at an initial SL_o = 2.0 μm in both high-loaded (isometric) and low-loaded isotonic contractions at different afterloads (Fig. 8).

Fig. 8

Contractions of ENDO (left) and EPI (right) cells under a high (isometric) and low loads (at afterload of 20, 40, 60, and 80% of maximal isometric force F _max) at isotonic mode of contractions. Tension (F) developed by cells (top panel) and cell deformations (bottom panel, ΔL in  % of the initial length L ₀) during the cardiac cycle in ENDO and EPI models. The stimulus was applied at 60 ms

In both cell types, a decrease in the afterload resulted in a prolongation of T _max (ENDO: 20% increase in the T _max during the afterloaded contraction at a low load of 20% F_max compared to T _max during isometric contraction, EPI: 14%) and a decrease in the time of relaxation (ENDO: 46% decrease, EPI: 40%, see Fig. 8), with an increase in APD₉₀ (ENDO: 10% increase, EPI: 6%). Thus, our models predict that ENDO cells demonstrate greater sensitivity to changes in afterload. We could not find any experimental data on transmural differences in the characteristics of contractions at different afterloads, but load dependence of contraction and relaxation was shown in isolated cells and in multicellular myocardial preparations [49,50,51]. The mechanisms underlying the load dependence of both T _max and TR₇₅ were analyzed in detail in our previous papers [41, 52]. Briefly, reduction in the afterload increased the amplitude and the velocity of sarcomere shortening, making CaTnC complexes dissociate faster due to the cooperative mechano-dependence of their kinetics, which resulted in a prolongation of Ca²⁺ transients [53]. Consistent with this, a decrease in afterload prolonged the APD, as seen experimentally [49]. Thus, our models account for specific transmural features of the cellular mechanics that modulates the electrical properties of the cells via mechano-calcium-electric feedback.

Discussion

Experimental data show that cardiac muscle is structurally and functionally heterogeneous, and that this heterogeneity manifests itself at all levels of functional integration, from the molecular to the whole-organ level [4, 12, 15, 16]. Studies on the transmural gradient in cellular mechanics and the possible effects of mechano-calcium-electric feedback on regional differences in electrophysiological function are still in their infancy.

Computational models of electromechanical coupling in ventricular myocytes are useful for investigating the mechanisms of functional differences across the LV wall. We developed mathematical models of the electromechanical coupling in cardiomyocytes from the ENDO and EPI regions of the guinea pig LV wall, which reproduce the physiological experimental data obtained in isolated myocardial cells. The ENDO and EPI models we developed are based on a widely validated Ekaterinburg–Oxford mathematical model (EO model), which simulates the electrical and mechanical activity of cardiomyocytes during isometric and afterloaded twitches under different cardiomyocyte lengths and mechanical loads. This model reproduces a wide range of the effects of cardiac mechano-calcium-electric feedback in healthy heart and under pathological conditions [10, 39, 54, 55]. When developing the next model generation, the ENDO and EPI model parameters were fitted quantitatively to available experimental data specific to guinea pig (Table 1). Unlike other models that have described transmural heterogeneity in the electrical properties of the myocardium [25,26,27], our models suggest the mechanisms which may explain the morphology of the AP in different myocardial regions by accounting for specific features of the mechanisms of excitation–contraction coupling that affect Ca²⁺ transients and the time course of contraction. The use of different parameters of Ca²⁺ handling and myofilament mechanics (Xb cycling) for ENDO and EPI models allowed us to reproduce a longer duration of the Ca²⁺ transient and a lower velocity of contraction and relaxation in ENDO against EPI cells, as observed experimentally (Table 4). Moreover, we showed that mechanical heterogeneity contributes essentially to the distinctions in the time course of the Ca²⁺ transient and, via mechano-calcium-electric feedback, slightly modulates AP properties of the cells.

Electrical heterogeneity

It was shown previously in canine and human hearts that heterogeneity in AP morphology reflects significant differences in ion current expression, particularly potassium currents [2, 4]. Significantly larger transient outward current (i _to) was observed in EPI and MID cells. This current is responsible for the prominent spike-and-dome morphology in these cells. ENDO cells have a much lower i _to that eliminates the notch and produces a longer APD at physiologic heart rates. Modeling results also considered i _to as a primary mechanism underlying transmural differences in cardiomyocyte function [29]. In the guinea pig ventricle, i _to is absent [56]. Our modeling results suggested that the higher i _CaL in guinea pig ENDO cardiomyocytes mainly underlies the longer APD and the longer contraction in ENDO cells (Table 2).

Ca²⁺ handling heterogeneity

A wide range of studies have described various aspects of spatial heterogeneity in the electrophysiological properties of cardiomyocytes and underlying ionic mechanisms. By contrast, less is known about the heterogeneity in the intracellular Ca²⁺ dynamics, which is a key mechanism of excitation–contraction coupling in cardiomyocytes. Experimental and simulation studies show that the regional differences in AP contribute to, but do not solely govern, the inherent differences in Ca²⁺ transients between ENDO and EPI myocytes [3, 18, 26, 38]. Some authors have suggested that the heterogeneity in Ca²⁺ handling is associated with heterogeneous SERCA2 expression and expression of ryanodine release channels [3, 13, 38]. In this study we have, for the first time to our knowledge, included in the mathematical model such transmural distinctions in calcium handling as smaller Ca²⁺ release from the SR in the ENDO region of the guinea pig [3].

We showed that changes only in the electrophysiological parameters did not allow the ENDO model to reproduce significantly slower Ca²⁺ transient against the EPI model (Table 3) as observed experimentally in isolated cardiomyocytes [3, 13, 38].

Our models predict that despite the lower SR Ca²⁺ content in the ENDO cells, a greater amount of Ca²⁺ released from the SR (Fig. 6) in the cells may provide for similar Ca²⁺ amplitudes in both cell types, which together with slower uptake of Ca²⁺ to the SR may result in a slower calcium transient in the ENDO compared to the EPI cell.

Mechanical heterogeneity

Passive and active mechanical properties of ventricular myocardium show significant spatiotemporal heterogeneity during the heartbeat, which is a prerequisite for normal cardiac activity. Recently, Ashikaga and co-workers have reported on significant transmural gradients in the amount and the time course of regional strains in canine hearts in vivo [15, 57]. However, much uncertainty remains, especially regarding the heterogeneity in the mechanical properties of the single cells. We used experimental data showing the significant differences in the passive mechanical properties obtained in isolated and skinned cardiomyocytes from guinea pigs [8, 9]. Cazorla and co-workers reported distinctions in the ‘SL-active tension’ relationship underlying the Frank–Starling law between the ENDO and EPI cells isolated from rat, ferret and guinea pig hearts [5, 9]. They have shown that the slope of the ‘SL-active tension’ relationship is steeper in ENDO myocytes, specifying their greater contractility. Similar studies have been recently carried out on the ENDO and EPI cells from the right and left ventricles of the guinea pig by P. Kohl’s group [58], who found no significant differences in the slopes of the “cell length-tension” curves between ENDO and EPI cardiomyocytes under auxotonic contractions. These controversial experimental data should be analyzed further within the models and experiments.

To our knowledge, there is the only one modeling study reproducing the transmural gradient in the mechanical activity of isolated cells from canines [29]. The authors suggested that Xb cycling kinetics primarily underlies transmural variation of myocyte contractile function [29]. This finding may be partially explained by the different regional combinations of myosin isoforms in the ventricular wall. Several studies suggest a transmural gradient in the expression ratio of the faster V1 isoform to the slower V3 isoform of myosin, with V1 expression being highest at the subepicardium compared to the subendocardium [6, 7]. Our simulation results are consistent with a previous modeling study [29]. We demonstrated that heterogeneity in Xb kinetics not only provides differences in contractile function between the ENDO and EPI cells but, via mechano-calcium-electric feedback, essentially enhances the transmural gradient in the duration of Ca²⁺ transient between the cells, which further may modulate electrical diversity.

To clarify the role of different intracellular mechanisms in the heterogeneity of cellular function, and to justify our conclusions, we assessed individual contributions of regional differences in ionic currents to the differences in the AP and mechanical characteristics between the ENDO and EPI cells (Table 2). We showed that transmural differences in ionic currents are important but not sufficient to produce the differences in the mechanical activity as experimentally observed (Table 4). Changes in the parameters of any individual ionic current from EPI to ENDO values produced about 10–15% difference in the temporal characteristics of contraction and relaxation between the models (Table 2), with highest response to the changes in the parameters of i_CaL, producing a difference in T _max and TR₇₅ between the cells of 13 and 22%, respectively. Note that the whole set of ENDO electrical parameters produced less effect on contraction-relaxation (6% for T _max and 13% for TR₇₅) than that of some individual currents, e.g., i _CaL and i _Ks. However, neither of the changes in the electrical parameters produced the roughly 25–30% difference in contraction characteristics as registered experimentally (Table 4).

Then, we compared the relative contribution of the differences in ionic currents with that combined with differences in calcium dynamics, and then additionally with the distinctions in the mechanical properties as well (Table 3). We revealed that distinctions in calcium handling and mechanical parameters slightly enhance the transmural gradient in AP characteristics between the ENDO and EPI cells (Table 3) and produce essential gradients in Ca²⁺ transient and contraction between the cells. Then, by excluding mechano-calcium-electric feedback in the ENDO model [Table 3, ENDO (−MEF)], we demonstrated explicitly that these feedback mechanisms underlie a significant mechanical effect on the duration of Ca²⁺ transient in the ENDO cell, which slightly affects the gradient in APD between the cells. These results emphasize that mechanical activity may be important not only as such for cardiac output, but also for the electrical and calcium signaling, which is altered via mechano-calcium-electric feedback. If this cross-link had been incorporated in the electrophysiological models, which account for the transmural differences in AP generation [25,26,27] referred to in our work, the incorporation would have led to a higher gradient in Ca²⁺ transient between the ENDO and EPI cells in all these models. Thus, at least parametrical refitting of these models would be needed after the hypothetical incorporation of the mechano-calcium-electric feedback, to keep their relevance even within the field of the electrophysiology and calcium handling. In our models, the whole set of distinctive electrical and mechanical parameters assigned to the ENDO and EPI cells allowed us to reproduce the combination of experimental data on the gradients in the profiles of AP, Ca²⁺ transient and contraction observed in isolated ENDO and EPI cells, and to explain underlying mechanisms.

There is the lack of experimental data on transmural differences in mechanical function of ENDO and EPI cardiomyocytes at different loading conditions. We could not find any experimental data in the literature about ENDO and EPI cardiomyocyte contraction at different afterloads. Our models predicted that ENDO and EPI cells responded differently to a decrease in the afterload, with greater load sensitivity of ENDO cells. This makes the modeling results a predictive source for further experimental verification. These data are important for analyzing load dependence of the whole heart where the dynamic mechanical interaction between cardiomyocytes leading to load redistribution should affect their activity in a way that depends crucially on the temporal activation sequence and inherent electromechanical properties of cells. Under normal physiologic conditions, heterogeneous myocardial elements confer system homogeneity to optimize cardiac function [53, 59]. Conversely, pathological conditions increase or decrease heterogeneity that may reduce cardiac efficiency [60]. These hypotheses may be analyzed further within the cardiac tissue models using the cellular models presented here.

Limitations

The models presented here have several limitations. One of these is that no account is made for temperature dependence. This problem arose due to the lack of a whole package of experimental data obtained at the same temperature. Most cellular model parameters were fitted and verified against data from experiments carried out at the physiological temperature of 34–37 °C (Table 1, Table 4). However, some characteristics, such as ‘SL-passive tension’ (Fig. 4), were specified at room temperature. Temperature can possibly modulate transmural differences in the electrical and mechanical functions of cells [61]; therefore, the models should be further verified against temperature effects.

In our study, we have focused on the mechano-calcium-electric feedback, such as cooperative Ca²⁺ activation of contractile proteins, and have not accounted for mechano-sensitive ionic channels and mechano-dependent SR Ca²⁺ release in the models, which may also contribute to the regulation of cellular contraction. We will include these pathways of mechano-electric feedback in future studies to analyze their roles in generation of a transmural gradient in myocardial properties.

Note that most of the numerical parameters of the models we developed were directly derived from experimental data obtained in guinea pig hearts, where such data were available (Table 1). Other species (dog, mouse, rat, ferret, etc.) are also compared in this table to show that the respective values are quite similar for various animals (for instance, passive stiffness, SERCA2). However, for fitting of i _pNa conductance and maximum i _NaK current (to reproduce experimental data on transmural differences in i _NaCa) we have used data from the canine heart, as such data is not available in the guinea pig. This is a rather commonly used approach in developing mathematical models, when a lack of experimental data for a particular species makes the authors use data from other species. Of course, such borrowing means that reasonable caution is necessary to quantitatively asses the modeling results. In our case, we believe that certain electrophysiological reasons could allow us to use data obtained in dogs for i _pNa and i _NaCa for a guinea pig model. In particular, it is known that the Na⁺ channel isoform Na_{v 1.5} is a major contributor to i _pNa [1, 62], and NCX₁ is a major contributor to i _NaCa [36, 63] in both guinea pig and canine ventricular cardiomyocytes, suggesting that the role of i _pNa and i _NaCa in AP generation is comparable in these species. The other source of model uncertainty is the use of qualitative data, which could not be applied directly to fit model parameters. An example of such parameters in our study is the sensitivity of the binding sites of ryanodine receptors to cytosolic Ca²⁺, or degree of RU end-to-end cooperativity (Table 1). In these cases, the parameter fittings were based on the overall model signals that could be compared to the experimental data (the time course of AP, contraction, etc.).

Conclusion

In this study we have developed detailed computational models of cardiomyocytes from different transmural regions based on experimental data. The models were validated against experimental data on the electrical and mechanical activity in isolated ENDO and EPI cells. The modeling results predict that heterogeneity in the parameters of calcium handling and myofilament mechanics in isolated ENDO and EPI cardiomyocytes via cooperative mechanisms of mechano-calcium-electric feedback are essential to produce the differences in Ca²⁺ dynamics and contraction profiles, and may further modulate transmural differences in the electrical properties between the cells. Our models predict that ENDO and EPI cells respond differently to changes in the afterload, with greater load sensitivity of ENDO cells. These data are important for understanding the behavior of cardiomyocytes in the whole heart.

Appendix

The single ENDO and EPI cardiomyocyte model equations and formulations based on the Ekaterinburg–Oxford mathematical model (EO model) [39] are presented here.

Electrical block of the model

Membrane currents

Fast Na⁺ current, i _Na

$$i_{\text{Na}} = g_{\text{Na}} \cdot m^{3} \cdot h \cdot (V - E_{\text{mh}} ),$$

$$\frac{{{\text{d}}m}}{{{\text{d}}t}} = \alpha_{m} \cdot (1 - m) - \beta_{m} \cdot m,$$

$$\frac{{{\text{d}}h}}{{{\text{d}}t}} = \alpha_{h} \cdot (1 - h) - \beta_{h} \cdot h,$$

$$\alpha_{m} = \left\{ {\begin{array}{*{20}l} {2000} \hfill & {{\text{if }}|V + 41| < \delta_{m} ,} \hfill \\ {\frac{200 \cdot (V + 41)}{{1 - e^{ - 0.1 \cdot (V + 41)} }}} \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$

$$\beta_{m} = 8000 \cdot e^{ - 0.056 \cdot (V + 66)} ,$$

$$\alpha_{h} = 20 \cdot e^{ - 0.125 \cdot (V + 75)} ,$$

$$\beta_{h} = \frac{2000}{{1 + 320 \cdot e^{ - 0.1 \cdot (V + 75)} }}.$$

Persistent Na⁺ current, i _pNa

$$i_{\text{pNa}} = \frac{{g_{\text{pNa}} }}{{1 + e^{{\frac{ - (V + 56.51)}{6.62}}} }} \cdot (V - E_{\text{Na}} ).$$

L-type Ca²⁺ current, i _CaL

$$i_{\text{CaL}} = d \cdot f \cdot f_{2,ds} \cdot (i_{\text{CaLCa}} + i_{\text{CaLNa}} + i_{\text{CaLK}} ),$$

$$i_{\text{CaLCa}} =\, P_{\text{Ca}} \cdot \frac{4 \cdot F}{R \cdot T}\frac{1}{{1 - e^{{\frac{ - 2 \cdot (V - 50) \cdot F}{R \cdot T}}} }}\left( {[{\text{Ca}}^{2 + } ]_{i} \cdot e^{{\frac{100 \cdot F}{R \cdot T}}} - [{\text{Ca}}^{2 + } ]_{o} \cdot e^{{\frac{ - 2 \cdot (V - 50) \cdot F}{R \cdot T}}} } \right) \cdot (V - 50),$$

$$i_{\text{CaLNa}} = 0.01 \cdot P_{\text{Ca}} \cdot \frac{F}{R \cdot T}\frac{1}{{1 - e^{{\frac{ - (V - 50) \cdot F}{R \cdot T}}} }}\left( {[{\text{Na}}^{ + } ]_{i} \cdot e^{{\frac{50 \cdot F}{R \cdot T}}} - [{\text{Na}}^{ + } ]_{\text{o}} \cdot e^{{\frac{ - (V - 50) \cdot F}{R \cdot T}}} } \right) \cdot (V - 50),$$

$$i_{\text{CaLK}} = 0.002 \cdot P_{\text{Ca}} \cdot \frac{F}{R \cdot T}\frac{1}{{1 - e^{{\frac{ - (V - 50) \cdot F}{R \cdot T}}} }}\left( {[K^{ + } ]_{i} \cdot e^{{\frac{50 \cdot F}{R \cdot T}}} - [K^{ + } ]_{\text{o}} \cdot e^{{\frac{ - (V - 50) \cdot F}{R \cdot T}}} } \right) \cdot (V - 50),$$

$$\frac{{{\text{d}}d}}{{{\text{d}}t}} = \alpha_{d} \cdot (1 - d) - \beta_{d} \cdot d,$$

$$\frac{{{\text{d}}f}}{{{\text{d}}t}} = \alpha_{f} \cdot (1 - f) - \beta_{f} \cdot f,$$

$$E_{{0{\text{d}}}} = V + 19,\;E_{{0{\text{f}}}} = V + 34,$$

$$\alpha_{\text{d}} = \left\{ {\begin{array}{*{20}l} {360} \hfill & {{\text{if }}|E_{{0{\text{d}}}} | < 0.0001,} \hfill \\ {\frac{{90 \cdot E_{{0{\text{d}}}} }}{{1 - e^{{ - \frac{{E_{{0{\text{d}}}} }}{4}}} }}} \hfill & {\text{otherwise}} \hfill \\ \end{array} ,} \right.$$

$$\beta_{\text{d}} = \left\{ {\begin{array}{*{20}l} {360} \hfill & {{\text{if }}|E_{{0{\text{d}}}} | < 0.0001,} \hfill \\ { - \frac{{36 \cdot E_{{0{\text{d}}}} }}{{1 - e^{{\frac{{E_{{ 0 {\text{d}}}} }}{10}}} }}} \hfill & {\text{otherwise}} \hfill \\ \end{array} ,} \right.$$

$$\alpha_{\text{f}} = \left\{ {\begin{array}{*{20}l} {7.5} \hfill & {{\text{if }}|E_{{0{\text{f}}}} | < 0.0001,} \hfill \\ { - \frac{{1.875 \cdot E_{{0{\text{f}}}} }}{{1 - e^{{\frac{{E_{{0{\text{f}}}} }}{4}}} }}} \hfill & {\text{otherwise}} \hfill \\ \end{array} ,} \right.$$

$$\beta_{\text{f}} = \frac{3.6}{{1 + e^{{\frac{{ - E_{{0{\text{f}}}} }}{4}}} }},$$

$$\frac{{{\text{d}}f_{{2,{\text{ds}}}} }}{{{\text{d}}t}} = 20 \cdot \left( {1 - \left( {\frac{{[{\text{Ca}}^{2 + } ]_{\text{ds}} }}{{0.001 + [{\text{Ca}}^{2 + } ]_{\text{ds}} }} + f_{{2,{\text{ds}}}} } \right)} \right).$$

Na⁺–Ca²⁺ exchange current, i _NaCa

$${i_{{{\text{NaCa}}_{\text{cmpt}} }} = {\text{Frac}}_{\text{cmpt}} \cdot i_{{{\text{NaCa}}\_{ \hbox{max} }}} \cdot \frac{{e^{{\frac{0.37 \cdot V \cdot F}{R \cdot T}}} \cdot [{\text{Na}}^{ + } ]_{i}^{3} \cdot [{\text{Ca}}^{2 + } ]_{\text{o}} - e^{{\frac{ - 0.63 \cdot V \cdot F}{R \cdot T}}} \cdot [{\text{Na}}^{ + } ]_{\text{o}}^{3} \cdot [{\text{Ca}}^{2 + } ]_{\text{cmpt}} }}{{1 + \frac{{[{\text{Ca}}^{2 + } ]_{\text{cmpt}} }}{0.0069}}}},$$

where cmpt is used to denote either cytosole (cyt) or dyadic space compartments (ds) of the cell;

$${\text{Frac}}_{\text{cyt}} = 0.999,\;\;{\text{Frac}}_{\text{ds}} = 0.001,$$

$$i_{\text{NaCa}} = i_{{{\text{NaCa}}_{\text{cyt}} }} + i_{{{\text{NaCa}}_{\text{ds}} }} .$$

Na⁺–K⁺ pump, i _NaK

$$i_{\text{NaK}} = i_{{{\text{NaK}}\_\hbox{max} }} \cdot \frac{{[{\text{K}}^{ + } ]_{\text{o}} }}{{1 + [{\text{K}}^{ + } ]_{\text{o}} }} \cdot \frac{{[{\text{Na}}^{ + } ]_{i} }}{{30 + [{\text{Na}}^{ + } ]}}_{i} .$$

Rapid delayed rectifier potassium current, i _Kr

$$i_{\text{Kr}} = x_{r} \cdot g_{\text{Kr}} \cdot \sqrt {\frac{{K_{0} }}{5.4}} \cdot \frac{{V - E_{K} }}{{1 + e^{{\frac{{V - V_{\text{oKr}} }}{{k_{\text{Kr}} }}}} }},$$

$$\alpha_{\text{xr}} = \frac{10}{{1 + e^{{\frac{{ - \left( {V - 10} \right)}}{370}}} }},$$

$$\beta_{\text{xr}} = A_{\text{Kr}} \cdot e^{{\frac{{ - \left( {V - 20} \right)}}{{k_{{ 1 {\text{Kr}}}} }}}} ,$$

$$\frac{{{\text{d}}x_{\text{r}} }}{{{\text{d}}t}} = \alpha_{\text{xr}} \cdot \left( {1 - x_{\text{r}} } \right) - \beta_{\text{xr}} \cdot x_{\text{r}} .$$

Slow delayed rectifier potassium current, i _Ks

$$i_{\text{Ks}} = \frac{{g_{\text{Ks}} }}{{\sqrt {[K]_{0} } }} \cdot x_{\text{s}}^{2} \cdot \left( {V - E_{\text{Ks}} } \right),$$

$$\alpha_{\text{xs}} = \frac{14}{{1 + A_{\text{Ks}} \cdot e^{{\frac{{ - \left( {V - 68} \right)}}{{B_{\text{Ks}} }}}} }},$$

$$\beta_{\text{xs}} = e^{{\frac{ - V}{88}}} ,$$

$$\frac{{{\text{d}}x_{\text{s}} }}{{{\text{d}}t}} = \alpha_{\text{xs}} \left( {1 - x_{\text{s}} } \right) - \beta_{\text{xs}} x_{\text{s}} .$$

Inward rectifier K⁺ current, i _K1

$$i_{\text{K1}} = g_{\text{K1}} \cdot \frac{{[{\text{K}}^{ + } ]_{\text{o}} }}{{[{\text{K}}^{ + } ]_{\text{o}} + 10}} \cdot \frac{{(V - E_{\text{K}} )}}{{1 + e^{{\frac{{1.25 \cdot \left( {V - E_{\text{K}} - 10} \right) \cdot F}}{R \cdot T}}} }}.$$

Transient outward K⁺ current, i _to

$$i_{\text{to}} = g_{\text{to}} \cdot s \cdot r \cdot (V - E_{\text{K}} ),$$

$$\frac{{{\text{d}}r}}{{{\text{d}}t}} = 333 \cdot \left( {\frac{1}{{1 + e^{{\frac{{ - \left( {V + 4} \right)}}{5}}} }} - r} \right),$$

$$\alpha_{\text{s}} = 0.033 \cdot e^{{\frac{ - V}{17}}} ,$$

$$\beta_{\text{s}} = \frac{33}{{1 + e^{{ - 0.125 \cdot \left( {V + 10} \right)}} }},$$

$$\frac{{{\text{d}}s}}{{{\text{d}}t}} = \alpha_{\text{s}} \cdot \left( {1 - s} \right) - \beta_{\text{s}} \cdot s.$$

Background Na⁺ current, i _bNa

$$i_{\text{bNa}} = g_{\text{bNa}} \cdot (V - E_{\text{Na}} ).$$

Background Ca²⁺ current, i _bCa

$$i_{\text{bCa}} = g_{\text{bCa}} \cdot (V - E_{\text{Ca}} ).$$

Reversal potentials, E _Na, E _K, E _Ks, E _Ca, E _mh

$$E_{\text{Na}} = \frac{R \cdot T}{F} \cdot \ln \frac{{[{\text{Na}}^{ + } ]_{\text{o}} }}{{[{\text{Na}}^{ + } ]_{i} }},$$

$$E_{\text{K}} = \frac{R \cdot T}{F} \cdot \ln \frac{{[{\text{K}}^{ + } ]_{\text{o}} }}{{[{\text{K}}^{ + } ]_{i} }},$$

$$E_{\text{Ks}} = \frac{R \cdot T}{F} \cdot \ln \frac{{[{\text{K}}^{ + } ]_{\text{o}} + P_{\text{KNa}} \cdot [{\text{Na}}^{ + } ]_{\text{o}} }}{{[{\text{K}}^{ + } ]_{\text{i}} + P_{\text{KNa}} \cdot [{\text{Na}}^{ + } ]_{i} }},$$

$$E_{\text{Ca}} = \frac{R \cdot T}{2 \cdot F} \cdot \ln \frac{{[{\text{Ca}}^{2 + } ]_{\text{o}} }}{{[{\text{Ca}}^{2 + } ]_{i} }},$$

$$E_{\text{mh}} = \frac{R \cdot T}{F} \cdot \ln \frac{{[{\text{Na}}^{ + } ]_{\text{o}} + 0.12 \cdot [{\text{K}}^{ + } ]_{\text{o}} }}{{[{\text{Na}}^{ + } ]_{\text{i}} + 0.12 \cdot [{\text{K}}^{ + } ]_{i} }}.$$

Intracellular and extracellular Na⁺ and K⁺ concentrations, [Na⁺] _i , [K⁺] _i , [K⁺]_o

$${\frac{{{\text{d}}[{\text{Na}}^{ + } ]_{i} }}{{{\text{d}}t}} = \frac{ - 1}{{V_{\text{Cell}} \cdot V_{{i\_{\text{ratio}}}} \cdot F}} \cdot \left( {i_{\text{Na}} + i_{\text{pNa}} + i_{\text{bNa}} + i_{\text{CaLNa}} + 3 \cdot i_{\text{NaK}} + 3 \cdot i_{\text{NaCa}} } \right),}$$

$$\frac{{{\text{d}}[{\text{K}}^{ + } ]_{i} }}{{{\text{d}}t}} = - \frac{1}{{V_{\text{Cell}} \cdot V_{{i\_{\text{ratio}}}} \cdot F}} \cdot \left( {i_{{{\text{K}}1}} + i_{\text{Kr}} + i_{\text{Ks}} + i_{\text{CaLK}} + i_{\text{to}} - 2 \cdot i_{\text{NaK}} } \right).$$

$${\frac{{{\text{d}}[{\text{K}}^{ + } ]_{\text{o}} }}{{{\text{d}}t}} = \frac{1}{{V_{\text{Cell}} \cdot V_{{e\_{\text{ratio}}}} \cdot F}} \cdot \left( {i_{{{\text{K}}1}} + i_{\text{Kr}} + i_{\text{Ks}} + i_{\text{CaLK}} + i_{\text{to}} - 2 \cdot i_{\text{NaK}} } \right) - 0.7 \cdot \left( {[{\text{K}}^{ + } ]_{\text{o}} - [{\text{K}}^{ + } ]_{\text{b}} } \right).}$$

Membrane potential, V

$${\frac{{{\text{d}}V}}{{{\text{d}}t}} = \frac{ - 1}{{C_{m} }} \cdot \left( {i_{\text{stim}} + i_{\text{Na}} + i_{\text{pNa}} + i_{\text{bNa}} + i_{\text{NaCa}} + i_{\text{CaL}} + i_{\text{bCa}} + i_{\text{to}} + i_{\text{Kr}} + i_{\text{Ks}} + i_{{{\text{K}}1}} + i_{\text{NaK}} } \right.),}$$

$$i_{\text{stim}} = - 3.0{\text{ nA }}\left( {{\text{on at 6}}0{\text{ ms}};{\text{ off at 62}}. 5 {\text{ ms}}} \right).$$

Calcium handling block

Ca²⁺ release from the SR, i _rel

$$i_{\text{rel}} = \left( {10000 \cdot \left( {\frac{{{\text{Fr}}_{\text{Act}} }}{{{\text{Fr}}_{\text{Act}} + 0.25}}} \right)^{2} + 0.05} \right) \cdot [{\text{Ca}}^{2 + } ]_{\text{rel}} ,$$

$$K_{\text{act}} = 500 \cdot {\text{RegBindSite}},$$

$$K_{\text{inact}} = 60 + 500 \cdot {\text{RegBindSite}},$$

$${\text{RegBindSite}} = \left( {\frac{{[{\text{Ca}}^{2 + } ]_{i} }}{{[{\text{Ca}}^{2 + } ]_{i} + K_{{m,{\text{Cacyt}}}} }} + \left( {1 - \frac{{[{\text{Ca}}^{2 + } ]_{i} }}{{[{\text{Ca}}^{2 + } ]_{i} + K_{{m,{\text{Cacyt}}}} }}} \right) \cdot \frac{{[{\text{Ca}}^{2 + } ]_{\text{ds}} }}{{[{\text{Ca}}^{2 + } ]_{\text{ds}} + 0.01}}} \right)^{2} ,$$

$${\text{SpeedRel}} = \left\{ {\begin{array}{*{20}l} 5 \hfill & {{\text{if}}\;V < - 50,} \hfill \\ 1 \hfill & {{\text{otherwise}},} \hfill \\ \end{array} } \right.$$

$${\text{Fr}}_{\text{Prec}} = 1 - {\text{Fr}}_{\text{Act}} - {\text{Fr}}_{\text{Prod}} ,$$

$$\frac{{{\text{dFr}}_{\text{Act}} }}{{{\text{d}}t}} = {\text{SpeedRel}} \cdot ({\text{Fr}}_{\text{Prec}} \cdot K_{\text{Act}} - {\text{Fr}}_{\text{Act}} \cdot K_{\text{inact}} ),$$

$$\frac{{{\text{dFr}}_{\text{Prod}} }}{{{\text{d}}t}} = {\text{SpeedRel}} \cdot ({\text{Fr}}_{\text{Act}} \cdot K_{\text{Inact}} - {\text{Fr}}_{\text{Prod}} ).$$

Ca²⁺ uptake by the SR-pump, i _up

$$K_{1} = 0.00012,$$

$$K_{2} = [{\text{Ca}}^{2 + } ]_{i} + [{\text{Ca}}^{2 + } ]_{\text{up}} \cdot K_{1} + 0.00021,$$

$$i_{\text{up}} = \alpha_{\text{up}} \cdot \frac{{[{\text{Ca}}^{2 + } ]_{i} }}{{K_{2} \cdot \left( {1 + \frac{{[{\text{Ca}}^{2 + } ]_{\text{up}} }}{4}} \right)}} - \beta_{\text{up}} \cdot \frac{{[{\text{Ca}}^{2 + } ]_{\text{up}} \cdot K_{1} }}{{K_{2} }}.$$

Ca²⁺ kinetics within the SR, i _trans

$$i_{\text{trans}} = 15 \cdot \left( {[{\text{Ca}}^{2 + } ]_{\text{up}} - [{\text{Ca}}^{2 + } ]_{\text{rel}} } \right),$$

Ca²⁺ concentration within the SR compartments, [Ca²⁺]_up, [Ca²⁺]_rel

$$\frac{{{\text{d}}[{\text{Ca}}^{2 + } ]_{\text{up}} }}{{{\text{d}}t}} = \frac{{V_{{i\_{\text{ratio}}}} }}{{V_{{{\text{up}}\_{\text{ratio}}}} }} \cdot i_{\text{up}} - i_{\text{trans}} ,$$

$$\frac{{{\text{d}}[{\text{Ca}}^{2 + } ]_{\text{rel}} }}{{{\text{d}}t}} = \frac{{\frac{{V_{{{\text{up}}\_{\text{ratio}}}} }}{{V_{{rel\_{\text{ratio}}}} }} \cdot i_{\text{trans}} - i_{\text{rel}} }}{{1 + \frac{{0.65 \cdot [{\text{CaS}}]_{\text{tot}} }}{{[{\text{Ca}}^{2 + } ]_{\text{rel}} + 0.65}}}}.$$

Ca²⁺ concentration within the cytosol and the DS, [Ca²⁺] _i , [Ca²⁺]_ds

$$\begin{aligned} \frac{{{\text{d}}[{\text{Ca}}^{2 + } ]_{i} }}{{{\text{d}}t}} & = \frac{ - 1}{{2 \cdot V_{\text{Cell}} \cdot V_{{i\_{\text{ratio}}}} \cdot F}} \cdot \left( {i_{{{\text{CaLCa}}_{\text{cyt}} }} + i_{\text{bCa}} - 2 \cdot i_{{{\text{NaCa}}_{\text{cyt}} }} } \right) + [{\text{Ca}}^{2 + } ]_{\text{ds}} \cdot 10 \cdot V_{{{\text{ds}}\_{\text{ratio}}}} \\ & \quad + \frac{{V_{{{\text{rel}}\_{\text{ratio}}}} }}{{V_{{i\_{\text{ratio}}}} }} \cdot i_{\text{rel}} - \frac{{d[{\text{CaTnC}}]}}{dt} - \frac{{d[B_{1} ]}}{dt} - \frac{{d[B_{2} ]}}{dt} - i_{\text{up}} , \\ \end{aligned}$$

$$\frac{{{\text{d}}[{\text{Ca}}^{2 + } ]_{\text{ds}} }}{{{\text{d}}t}} = - \frac{1}{{2 \cdot V_{{{\text{ds}}\_ratio}} \cdot V_{\text{Cell}} \cdot V_{{i\_{\text{ratio}}}} \cdot F}} \cdot (i_{{{\text{CaLCa}}_{\text{ds}} }} - 2 \cdot i_{{{\text{NaCa}}_{\text{ds}} }} ) - 10 \cdot [{\text{Ca}}^{2 + } ]_{\text{ds}} .$$

Kinetics of CaTnC complexes, [CaTnC]

$$\frac{{{\text{d}}[{\text{CaTnC}}]}}{{{\text{d}}t}} = k_{\text{on}} \cdot \left( {[{\text{TnC}}]_{\text{tot}} - [{\text{CaTnC}}]} \right) \cdot [{\text{Ca}}^{2 + } ]_{i} - k_{\text{off}} \cdot e^{{ - k_{A} \cdot [{\text{CaTnC}}]}} \cdot \pi_{{N_{\text{A}} }} \cdot [{\text{CaTnC}}].$$

A key point of this equation is the cooperativity of ${\text{Ca}}^{2 + }$ affinity to TnC (see the main text). In accordance with Xbs–CaTnC cooperativity (Fig. 2), we consider CaTnC off-rate to be decreasing with increases in the number of Xbs (N). The dependence $\pi_{{N_{\text{A}} }}$ expresses this cooperativity in the above equation, and its particular form may be written as follows:

$$\pi_{{N_{\text{A}} }} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {{\text{if}}\;N_{\text{A}} \le 0,} \hfill \\ {0.03^{{N_{\text{A}} }} } \hfill & {{\text{if}}\;0 < N_{\text{A}} \le 1,} \hfill \\ {0.03} \hfill & {{\text{otherwise}},} \hfill \\ \end{array} } \right.$$

where $N_{\text{A}}$ means an average fraction of the attached Xbs falling on one CaTnC complex; i.e. $N_{\text{A}} = \frac{N}{{[{\text{CaTnC}}]}}$.

In accordance with CaTnC–CaTnC cooperativity (Fig. 2), we consider CaTnC off-rate to be decreasing with increases in [CaTnC]. Dependence $e^{{ - k_{\text{A}} \cdot [{\text{CaTnC}}]}}$ defines this cooperativity in the above equation, where $k_{\text{A}}$ is the model parameter (B1, Table 4).

Ca²⁺ buffering system, B ₁, B ₂

$$\frac{{{\text{d}}[B_{1} ]}}{{{\text{d}}t}} = b_{{1_{\text{on}} }} \cdot \left( {[B_{1} ]_{\text{tot}} - [B_{1} ]} \right) \cdot [{\text{Ca}}^{2 + } ]_{i} - b_{{1_{\text{off}} }} \cdot [B_{1} ],$$

$$\frac{{{\text{d}}[B_{2} ]}}{{{\text{d}}t}} = b_{{2_{\text{on}} }} \cdot \left( {[B_{2} ]_{\text{tot}} - [B_{2} ]} \right) \cdot [{\text{Ca}}^{2 + } ]_{i} - b_{{2_{\text{off}} }} \cdot [B_{2} ].$$

Mechanical block of the model

Kinetics of force generating Xbs, N

$$\frac{{{\text{d}}N}}{{{\text{d}}t}} = k_{{p_{v} }} \cdot M([{\text{CaTnC}}]) \cdot n_{1} (l_{1} ) \cdot L_{\text{oz}} \cdot (1 - N) - k_{{m_{v} }} \cdot N,$$

$M([{\text{CaTnC}}]) = \frac{{\left( {\frac{{[{\text{CaTnC}}]}}{{[{\text{TnC}}]_{\text{tot}} }}} \right)^{\mu } \cdot \left( {1 + 0. 6^{\mu } } \right)}}{{\left( {\frac{{[{\text{CaTnC}}]}}{{[{\text{TnC}}]_{\text{tot}} }}} \right)^{\mu } +\, 0. 6^{\mu } }}$ means RU end-to-end interaction between adjacent tropomyosin segments in the case where both of them are affected by the formation of respective ${\text{CaTnC}}$ complexes (Eq. 5, Fig. 2), where µ is the cooperativity degree (B1., Table 4).

$$n_{1} (l_{1} ) = \left\{ {\begin{array}{ll} 0 &\quad {{\text{if}}\;\;0.6 \cdot l_{1} + g_{2} \le 0,} \\ {0.6 \cdot l_{1} + g_{2} } &\quad {{\text{if}}\;\;0 < 0.6 \cdot l_{1} + g_{2} < 1,} \\ 1 &\quad {{\text{otherwise}} .} \\ \end{array} } \right.$$

is the probability of a myosin cross-bridge “finding” a vacant site on the actin filament; in other words, the dependence of $n_{1}$ on $l_{1}$ means sarcomere lattice spacing.

$$L_{\text{oz}} (l_{1} ) = \left\{ {\begin{array}{*{20}c} {\frac{{l_{1} + 1.14}}{1.6}} & {{\text{if}}\;\;l_{1} \le 0.55,} \\ {\frac{1.69}{1.6}} & {{\text{otherwise}}.} \\ \end{array} } \right.$$

is a normalized linear dependence of the sarcomere overlap zone on the SL; $k_{{p_{v} }}$ and $k_{{m_{v} }}$ are velocity-dependent rate constants of cross-bridge attachment and detachment, respectively:

$$k_{{p_{v} }} = 0.6822 \cdot \chi_{0} \cdot q_{v} \cdot G^{*} ,$$

$$k_{{m_{v} }} = \chi_{0} \cdot q_{v} \cdot (1 - 0.6822 \cdot G^{*} ),$$

$$q_{v} = \left\{ {\begin{array}{ll} {17 - \frac{259 \cdot v}{{v_{ \hbox{max} } }}} &\quad {{\text{if}}\;\;v \le 0,} \\ {\frac{ - 2.3 \cdot v}{{0.79 \cdot v_{ \hbox{max} } }} + 17.3} &\quad {{\text{if}}\;\left( {0 < v \le 0.79 \cdot v_{ \hbox{max} } } \right),} \\ {\frac{15}{{\left( {1 + \frac{{ 5\left( {v - 0.79 \cdot v_{ \hbox{max} } } \right)}}{{v_{ \hbox{max} } }}} \right)^{10} }}} &\quad {\text{otherwise}} \\ \end{array} } \right.$$

Mechanical variables, F _CE, F _VS, F _SE, F _PE, F _XSE, l, l ₁, l _ex

The following formula gives the force developed by the activated contractile element CE in the model (Fig. 1):

$$F_{\text{CE}} = \lambda \cdot N \cdot p(v),$$

Function $p(v)$ means the dependence of the average Xb force on the sarcomere shortening/lengthening velocity (v):

$$p(v) = \frac{{P^{*} (v)}}{{G^{*} (v)}},$$

where $P^{*} (v)$ is the dependence of the steady-state sarcomere force on the shortening/lengthening velocity, and $G^{*} (v)$ is the dependence of the steady-state sarcomere stiffness on the velocity. Both functions are normalized so that $P^{*} (0) = 1$, $G^{*} (0) = 1$.

$$P^{*} = \left\{ {\begin{array}{ll} {\frac{{0.25 \cdot \left( {1 + \frac{v}{{v_{\hbox{max} } }}} \right)}}{{0.25 - \frac{v}{{v_{\hbox{max} } }}}}} &\quad {{\text{if}}\;v \le 0,} \\ {1.5 - \frac{0.0625}{{25 \cdot \left( {\frac{v}{{v_{\hbox{max} } }}} \right)^{2} + \frac{1.25 \cdot v}{{v_{\hbox{max} } }} + 0.125}}} &\quad {{\text{otherwise}};} \\ \end{array} } \right.$$

$$G^{*} = \left\{ {\begin{array}{ll} {1 + \frac{0.6 \cdot v}{{v_{\hbox{max} } }}} \hfill &\quad {{\text{if}}\left( { - v_{\hbox{max} } \le v} \right)\;\rm and\;\left( {v \le 0} \right),} \hfill \\ {\frac{{P^{*} }}{{\frac{1.1 \cdot v}{{0.25 \cdot v_{\hbox{max} } }} + 1}}} \hfill &\quad {{\text{if}}\left( {0 < v \le 0.1 \cdot v_{\hbox{max} } } \right),} \hfill \\ {\frac{{P^{*} \cdot e^{{ - \left( {v - 0.1} \right)^{4} }} }}{{\frac{1.1 \cdot v}{{0.25 \cdot v_{\hbox{max} } }} + 1}}} \hfill &\quad {{\text{otherwise}}.} \hfill \\ \end{array} } \right.$$

The force generated by the damper VS during cell deformation (Fig. 1):

$$F_{\text{VS}} = \beta_{v} \cdot e^{16 \cdot l(t)} \cdot \frac{{{\text{d}}l}}{{{\text{d}}t}}(t).$$

The following equations define the SE and PE elements of the rheological schemes underlying the mechanical module of the model:

$$F_{\text{SE}} = \beta_{1} \cdot (e^{{25 \cdot (l(t) - l_{1} (t))}} - 1),$$

$$F_{\text{PE}} = \beta_{2} \cdot (e^{{\alpha_{2} \cdot l(t)}} - 1).$$

The following equations define the force F that is developed by the cardiomyocyte according to the rheological scheme:

$$F_{\text{CE}} = F_{\text{SE}} ,$$

$$F = F_{\text{CE}} + F_{\text{PE}} + F_{\text{VS}} .$$

The force of extra-series element XSE is defined as:

$$F_{\text{XSE}} = \beta_{3} \cdot (e^{{48 \cdot l_{\text{ex}} (t)}} - 1).$$

Thus, at each t, the functions l(t), l ₁(t) and v(t) are solutions of the system of differential–algebraic equations:

$$\lambda \cdot N \cdot p(v) = \beta_{1} \cdot (e^{{25 \cdot (l - l_{1} )}} - 1),$$

$$\frac{{{\text{d}}l}}{{{\text{d}}t}} = \frac{{\beta_{3} \cdot (e^{{48 \cdot l_{\text{ex}} }} - 1) - \beta_{1} \cdot (e^{{25 \cdot (l - l_{1} )}} - 1) - \beta_{2} \cdot (e^{{\alpha_{2} \cdot l}} - 1)}}{{\beta_{v} \cdot e^{16 \cdot l} }},$$

$$\frac{{{\text{d}}l_{1} }}{{{\text{d}}t}} = v.$$

For the single cardiomyocyte model, these equations are complemented with the isometric condition:

$$l + l_{\text{ex}} \equiv {\text{const}}.$$

ENDO and EPI cell model parameters (Tables 5, 6, 7, 8, 9)

All the initial values for the ENDO and EPI phase variables are shown in their steady state, achieved by allowing the models to run for 100 s at a 1-Hz pacing rate.

Table 5 Cellular model parameters

Table 6 Parameters for the electrical block

Table 7 Parameters for the Ca²⁺ handling

Table 8 Parameters for the mechanical block

Table 9 Phase variables