Transmural cellular heterogeneity in myocardial electromechanics.

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 dif- ferences in Ca 2? 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.

Cardiac transmural heterogeneity /C1 Electromechanical coupling /C1 Mechano-calcium-electric feedback /C1 Cardiac modeling /C1 Cardiomyocyte

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 ven- tricular 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 car- diomyocytes responsible for their ability to shorten and generate the force, together assemble the cellular mechanical heterogeneity. The distinctions include differ- ent 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 2? sensitivity of myofilaments [ 8] also contribute to the mechanical heterogeneity, since Ca 2? activation of the crossbridges directly affects cardiomyocyte contractility. All such specific regional features of the cellular elec- tromechanical coupling in cardiomyocytes together con- tribute to distinctions in integrative functional characteristics of cardiomyocytes, such as AP morphology and duration [ 4], ‘force–velocity’ and ‘force–length’ rela- tionships [ 5, 8, 9], and so on. Importantly, not only does electrical intracellular heterogeneity directly contribute to the AP distinctions between the cells, but mechanical & Anastasia Khokhlova [email protected]; [email protected] 1 Ural Federal University, Ekaterinburg, Russia 2 Institute of Immunology and Physiology, Russian Academy of Sciences, 106 Pervomayskaya, Ekaterinburg 620049, Russia 3 Okayama University, Graduate School of Medicine, Dentistry and Pharmaceutical Sciences, 2-5-1 Shikata-cho, Kita-ku, Okayama 700-8558, Japan 123 J Physiol Sci (2018) 68:387–413 https://doi.org/10.1007/s12576-017-0541-0 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 subepi- cardial (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-fre- quency 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 signifi- cant 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 prop- erties of cardiomyocytes, much less is known about the heterogeneity in the mechanical function of these cells and the underlying mechanisms of excitation–contraction cou- pling. 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. Sev- eral mathematical models take into account the transmural gradient in the electrical properties of the myocardium at different levels of myocardial tissue organization [ 24–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–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 2? release from the sarcoplasmic reticulum [ 3], a greater myofilament Ca2? 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 con- trast to other models simulating either unloaded cell shortening [ 29], or isometric twitches at a constant sar- comere length [ 28], we have studied effects of intracellular heterogeneity on the specific features of cellular contrac- tion under more physiological conditions of imposed mechanical load, i.e., simulating afterloaded cellular con- tractions. 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 devel- oping 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 ?/Ca2? 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 compo- nents of delayed K ? current ( iKr and iKs) 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 iNaCa in isolated ENDO and EPI myocytes [36]. Following the experimental findings, we suggested a higher iNaK amplitude in the ENDO model and showed a resulting smaller density of iNaCa 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 388 J Physiol Sci (2018) 68:387–413 123 Table 1 Experimental data on the transmural differences in intracellular electrical and mechanical properties between the ENDO and EPI cells taken into acc ount in the cell subtype models ENDO EPI Animal species T, /C176C Source Target model parameter Parameter in the Appendix Nav protein expression Higher Lower Guinea pig 4 [ 1] Maximum membrane ipNa conductance Table 6 ipNa Higher Lower Dog 37 [ 43] iKr,i Ks Lower Higher Guinea pig 35 ± 1[ 11] iKr,i Ks current–voltage parameters Table 6 iNaCa Lower Higher Dog 37 [ 44] Maximum iNaK, current Table 6 Mouse 37 [ 18] DHPa1 Higher Lower Guinea pig 4 [ 3] Membrane permeability for Ca 2? through L-type channels, maximum membrane ibCa conductance Table 6 Ca2? release from the SR Lower Higher Guinea pig 4 [ 3] Sensitivity of the binding sites of ryanodine receptors to the cytosolic Ca 2? Table 7 SERCA2 Lower Higher Dog 4 [ 38] Rate constant of SERCA Ca 2? uptake Table 7 Guinea pig 4 [ 3] Human 4 [ 13] Passive stiffness Higher Lower Guinea pig 26–27 [ 9] Scaling coefficient for passive force Table 8 Rat, ferret 24–27 [ 5] Rate constant of Xb cycling Lower Higher Guinea pig 22 [ 8] Rate constant for Xb attachment, Xb attachment-detachment ratio coefficient Table 8 Pig 22 [ 7] Ca2? sensitivity of force Higher Lower Guinea pig 22 [ 8] Degree of the RU end-to-end cooperativity Table 8 J Physiol Sci (2018) 68:387–413 389 123 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 cardiomy- ocytes 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 experi- mental findings on the transmural mechanical heterogene- ity between the cells (Table 1). Our study is aimed at the development of the integrative electromechanical models of the ENDO and EPI car- diomyocytes 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.

EO model of cellular electromechanics as a predecessor of ENDO and EPI models Our ENDO and EPI models of the electrical and mechan- ical activity in cardiomyocytes from ENDO and EPI 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 ( ix) involved in the AP generation are shown: iNa fast Na ? current, ipNa persistent Na? current, ibNa background Na? current, iNaCa Na?/Ca2? exchange (NCX) current, iCaL Ca2? current via L-type Ca 2? channels, ibCa background Ca 2? current, ito transient outward K ? current, iKr, iKs rapid and slow components of delayed K ? current, iK1 rectifier K? current, iNaK Na?/K? pump current, DS dyadic space, TC terminal cisterns, LR longitudinal reticulum, irel Ca2? release flow from the sarcoplasmic reticulum (SR), iup SERCA pump flow from cytosol to the SR, itr Ca2? diffusion between the LR and the TC, CaS Ca2? complexes with calsequestrin, CaTnC Ca2? complexes with TnC, B1, B2 Ca2? 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 1 and lex define deformations of PE, CE and XSE, respectively, relative to their slack lengths (see text for details) 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. Cooper- ative mechanisms of calcium activation of myofilament RUs (Xb– CaTnC, CaTnC–CaTnC and RU end-to-end cooperativity) are described in the text 390 J Physiol Sci (2018) 68:387–413 123 regions are based on an Ekaterinburg–Oxford mathemati- cal 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. An electrophysiological block of the EO model is based on the Noble 098 ionic model [ 40]. The transmembrane potential is calculated with the following equation: where Cm is the membrane capacitance, istim is the stimu- lation current, iNa is the fast Na ? current, ipNa is the per- sistent Na ? current, ibNa is the background Na ? current, iNaCa Na?/Ca2? is the exchange current, iCaL Ca2? is the current via L-type Ca 2? channels, ibCa is the background Ca2? current, ito is the transient outward K ? current, iKr, iKs are the rapid and slow components of the delayed K ? current, iK1 is the rectifier K ? current, and iNaK Na?/K? is the pump current (Fig. 1a). Figure 2 shows a scheme of links between mechanisms of the excitation contraction coupling and mechano-cal- cium-electric feedback taken into account in the EO model. AP generation governs cytosolic Ca 2? transient (direct links, solid lines) and in turn depends on Ca 2? kinetics via transmembrane Ca 2? fluxes (feedback links, dashed lines). The equation for the concentration of free cytosolic Ca 2? ([Ca2?]i) describes Ca2? 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 2? binding to myofilament regulatory protein tro- ponin C (TnC) and other Ca2? binding ligands (Fig.1a) [39]: d½Ca2þ/C138 i dt ¼ /C0 1 2 /C1 VCell /C1 Vi ratio /C1 F /C1 iCaLCacyt þ ibCa /C0 2 /C1 iNaCacyt /C0/C1 þ½ Ca2þ/C138 ds /C1 10 /C1 Vds ratio þ Vrel ratio Vi ratio /C1 irel /C0 d½CaTnC/C138 dt /C0 d½B1/C138 dt /C0 d½B2/C138 dt /C0 iup; ð2Þ where VCell is the cell volume, Vi_ratio is the volume frac- tion of cytosolic space, Vds_ratio is the volume ratio of DS, Vrel_ratio is the volume ratio of the SR release site, F is Faraday’s constant, iCaLCacyt is the Ca 2? current via L-type Ca2? channels, ibCa is the background Ca 2? current, iNaCacyt is the Na ?/Ca2? exchange current, [Ca 2?]ds is the Ca2? concentration within the DS, irel is the Ca 2? release flux from the SR, [CaTnC] is the concentration of Ca 2? complexes with TnC, [ B1], [ B2] are the concentrations of Ca2? complexes with fast and slow binding ligands, and iup is the sarcoplasmic reticulum (SR) Ca 2? ATPase (SERCA) pump flux from the cytosol to the SR (Fig. 1a). Note that [Ca2?]ds,10,Vds_ratio in Eq. ( 2) is a term used in the model Noble098 [40] to describe Ca2? diffusion from DS to cytosol. Strictly speaking, the formula D ([Ca2?]ds - [- Ca2?]i) would be a more correct description of the diffusion between DS and the cytosol (with diffusion coefficient D)i n this equation, allowing for a reverse flow if Ca 2? concen- tration 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 098, probably owing to the passing of the input calcium flow into the cytosol only through the DS. Therefore, ‘ ‘uni- directional diffusion’ ’ of Ca2? from the DS to the cytosol in Eq. ( 2) seems to be a permissible simplification inherent from the model Noble 098, 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 jus- tified 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 disso- ciation 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 2? binding by TnC, located within one regulatory unit (RU) of actin–tropo- myosin–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– dV dt ¼ /C0 1 Cm istim þ iNa þ ipNa þ ibNa þ iNaCa þ iCaL þ ibCa þ ito þ iKr þ iKs þ iK1 þ iNaK /C0 Þ; istim ¼/C0 3:0 nA every cardiac cycle on at 60 ms ; off at 62 :5m sðÞ : ð1Þ J Physiol Sci (2018) 68:387–413 391 123 CaTnC cooperative mechanism), making Ca 2? kinetics mechano-sensitive as well [ 39]: d½CaTnC/C138 dt ¼ kon ½TnC/C138 tot /C0½ CaTnC/C138 /C0/C1 /C1½ Ca2þ/C138 i /C0 koff ð½CaTnC/C138 ; NÞ/C1½ CaTnC/C138 ; ð3Þ where kon and koff ð½CaTnC/C138 ; NÞ are on- and off-rate ‘ ‘con- stants’ ’ 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 cardiomy- ocytes (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 mecha- nisms 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 sar- comere ensemble of the cardiomyocyte. Sarcomeres gen- erate 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 2? 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 proper- ties 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 1 and shortening velocity v][ 39]: dN dt ¼ konð½CaTnC/C138 ; l1; vÞ/C1ð 1 /C0 NÞ/C0 koff ðvÞ/C1 N; ð4Þ where konð½CaTnC/C138 ; l1; vÞ and koff ðvÞ¼ kmv 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ð½Ca TnC/C138Þ ¼ ½CaTnC/C138 ½TnC/C138 tot /C16/C17 l /C1 1 þ 0:6lðÞ ½CaTnC/C138 ½TnC/C138 tot /C16/C17 l þ 0:6l ; ð5Þ where [CaTnC] is the concentration of CaTnC complexes, [CaTnC]tot is the total TnC concentration, and l 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 2? activation in the model. These mechanisms underlie a number of essential properties of cardiac muscle, such as the depen- dence of contraction and relaxation on myofilament length and load. Suppose that l 1 is the relative change in the SL of the cell against its slack length (deformations of the CE, nor- malized 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 ( FCE) is defined by the fraction of force-generating Xbs N and depends on the velocity of the CE shortening/stretching ( v ¼ dl1 dt ) (Fig. 2)[ 39]: FCE ¼ k /C1 N /C1 pðvÞ; ð6Þ where k 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 PE ¼ b2 /C1 ea2 /C1 l /C0 1 /C0/C1 ; ð7Þ where b2, a2 are model parameters (for more details, see the ‘ ‘Appendix’ ’). The elastic element SE and the damper VS also modu- late the passive and active mechanical properties [ 39]. The element SE generates a force proportional to the difference between l and l 1. 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 reg- istered experimentally (Table 1). In the next Section, ‘ ‘ENDO and EPI cellular model development’ ’, we describe selection of model parameters, which differ in the 392 J Physiol Sci (2018) 68:387–413 123 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 expres- sion 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 suben- docardium than in the subepicardium in both heart ventri- cles 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 subendo- cardium than in the subepicardium in canine hearts [ 43]. Based on these data, the maximum membrane ipNa con- ductance 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 0m V [ 43] (see the Appendix ‘ ‘ Persistent Na ? current, ipNa’ ’, Appendix Table 6). Fast and slow components of delayed potassium current, iKr and iKs 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 iKr and iKs. They found that the density of iKr and iKs 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 cur- rents in the guinea pig, reproducing the current–voltage characteristics of i Kr and iKs as registered experimentally in the EPI and ENDO cells [ 34]. Here, we used these findings in the cellular models (see Appendix ‘ ‘ Rapid delayed rec- tifier potassium current, iKr’’ , ‘‘Slow delayed rectifier potassium current, iKs’ ’, Appendix Table 6). L-type calcium current, iCaL No differences were found in the current–voltage charac- teristics of iCaL between ENDO and EPI cells in the mouse [18] and dog [ 45]. Nonetheless, the expression of the L-type calcium channel protein DHP a1 was greater in the subendocardium than in the subepicardium in the guinea pig LV [ 3]. According to these data, the membrane per- meability for Ca 2? through L-type channels was set as larger in the ENDO model than in the EPI model (see Appendix ‘ ‘L-type Ca2? current, iCaL’ ’, Appendix Table 6). The slightly higher maximum membrane ibCa conductance chosen for the ENDO cell maintains a greater Ca 2? current into the cell (see Appendix ‘ ‘ Background Ca 2? current, ibCa’ ’, Appendix Table 6). Na1/Ca21 exchange and Na 1/K1 pump currents, iNaCa and iNaK Experimental data on the regional differences in iNaCa and iNaK in LV are still controversial. Laurita and co-workers found no significant differences in the expression of the dominant isoform of NCX 1 in the dog [ 38], while other authors registered iNaCa as smallest in ENDO cells in canine and mouse hearts [ 18, 44]. Quinn and co-workers have demonstrated a higher iNaCa 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 gra- dient in the corresponding Na ? concentration [Na ?]i [36] Fig. 3 Ratio of maximal rate constant of force redevelopment Ktr in EPI to ENDO cells registered in experiments [ 8] and produced by the models at SL 0 = 2.3 lm. Ktr was determined by monoexponential fitting to the maximal tension for each pCa (= -log[Ca2?]). Ratio of mean values of Ktr from experimental data [ 8] is used. The latter data we calculated as a ratio of the Ktr values corresponding to equal values of pCa [see left and right plots in Fig. 3b (for SL 0 = 2.3 lm) in the cited article] J Physiol Sci (2018) 68:387–413 393 123 and suggested that this transmural gradient in iNaCa 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 iNaK underlies the heterogeneity in iNaCa in isolated ENDO and EPI myocytes. We assigned a higher iNaK amplitude in the ENDO model and showed a resulting smaller density of iNaCa 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, iNaK’’ , Appendix Table 6). This mechanism is utilized in the cellular models presented here. Heterogeneity in Ca 21 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 2? uptake was given a smaller value in the ENDO model than in the EPI model (see Appendix ‘ ‘Ca2? uptake by the SR-pump, iup’’ , Appendix Table 7). Wan and co-workers suggested that less Ca2? 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 2? in the ENDO model (see Appendix ‘ ‘ Ca2? release from the SR, irel’ ’, 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]. Trans- mural differences in the active mechanical properties of cardiomyocytes may be associated with regional differences in Xb cycling kinetics. Ait Mou and co-workers demon- strated in skinned cardiomyocytes of guinea pig LV at sat- urated Ca 2? concentration that maximal rate constant of force redevelopment Ktr 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]. 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, FVS, FSE, FPE, FXSE, l, l1, lex’ ’, Appendix Table 8) and a higher rate constant of Xb cycling in the EPI model (see Appendix ‘ ‘ Mechanical variables, FCE, FVS, FSE, FPE, FXSE, l, l1, lex’ ’, Appendix Table 8). These modulations in model parameters allowed the models to reproduce qualitatively experimental data on higher Ktr in EPI cardiomyocytes (Fig. 3). Note that quantitative mis- match between experimental data and simulations may be attributed to the differences in parameters of skinned car- diomyocytes used in experiments and that of simulated intact cardiomyocytes. Ait Mou and co-workers showed in skinned cells that Ca 2? sensitivity of force at SL of about 2.3 lm was higher in ENDO cells than in EPI cells, suggesting a greater molecular cooperativity of myofilament Ca 2? activation in the ENDO cells [ 8]. Based on these data, we set a different degree of nonlinearity l for the RU end-to-end coopera- tivity 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 l in Eq. 5 for the ENDO model makes the effects of RU end-to-end inter- action 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 2?-force relationship in the ENDO cells at the physiological range of cytosolic Ca 2? concentrations as shown experimentally. The modeling

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 sev- eral 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 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) 394 J Physiol Sci (2018) 68:387–413 123 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 cardiomy- ocytes from the guinea pig LV [ 9] (Fig. 4 and see Appendix ‘ ‘ Mechanical variables, F CE, FVS, FSE, FPE, FXSE, l, l1, lex’ ’, Appendix Table 8). 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 SLo of 2.0 lm. 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 after- loaded isotonic contractions were examined.

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 Table 2 Contribution of ENDO/EPI differences in individual ionic currents to differences in action potential, Ca 2? transient, and contraction between the ENDO/EPI cellular models Ionic current RP APD 10 APD50 APD90 CaT70 CaT90 Tmax TR75 ipNa 0.99 1.08 1.08 1.08 1.00 0.97 1.08 1.13 iKr 0.99 1.06 1.06 1.06 1.00 0.97 1.07 1.13 iKs 0.99 1.21 1.18 1.17 1.00 0.97 1.09 1.19 iCaL 0.99 1.35 1.13 1.11 1.00 0.97 1.13 1.22 ibCa 1.00 1.05 1.05 1.05 1.00 0.97 1.08 1.13 iNaK/iNaCa 0.99 1.20 1.11 1.10 1.00 0.97 1.03 1.08 Total 0.99 1.31 1.26 1.25 1.00 0.97 1.06 1.14 The ENDO/EPI ratio between the characteristics produced by the models at isometric twitches is shown RP resting membrane potential; APD10,50,90 action potential duration at 10, 50 and 90% repolarization, respectively; CaT70,90 calcium transient duration at 70 and 90% decay, respectively; Tmax time to peak contraction; TR75 time to 75% relaxation (measured as the time interval from Tmax to the time of 75% relaxation) 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 J Physiol Sci (2018) 68:387–413 395 123 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 parame- ters set the same as in the complete ENDO model (see the Appendix Table 6, ENDO parameters). 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 dif- ferences 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 iNaK between the models is viewed as negligible (Fig. 5), this difference led to the regional heterogeneity in iNaCa,a s mentioned above (see the Na ?/K? pump and Na ?/Ca2? exchanger section). The modeling results suggested that the smaller outward iNaCa 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 iKr and iKs and the smaller inward ipNa 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]. To reveal the contribution of the heterogeneity in cal- cium 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 cal- cium handling properties between the cells (see Appendix Tables 6, 7), and (c) the complete ENDO model account- ing 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 2? transient but ensured slower contraction (Table 3). However, the gra- dient 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 2? transient (CaT 90) due to slower calcium uptake from the cytosol, and essentially delayed Tmax and TR75 in the ENDO cell (Table 3). Finally, the change in the mechanical parameters for the complete ENDO model produced the most prolonged APD 90, CaT 90, and Tmax (Table 3). As compared to the ENDO (E ? Ca) model, delayed contraction in the complete ENDO model essen- tially 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 Table 3 Contribution of ENDO/EPI differences in the electrical, calcium handling and mechanical properties to the characteristics of action potential, Ca 2? transient, and contraction at isometric twitches Model APD 90, ms CaT 90,m s Tmax,m s T R 75,m s EPI 204 71 115 131 ENDO (E) 256/1.25 69/0.97 122/1.06 149/1.14 ENDO (E ? Ca) 259/1.27 80/1.12 140/1.22 176/1.34 ENDO 262/1.28 97/1.37 154/1.34 175/1.34 ENDO ( -MEF) 256/1.25 82/1.15 158/1.37 178/1.36 The ENDO/EPI ratio for the characteristics produced by the models is shown after the slash APD 90 action potential duration at 90% repolarization, CaT90 calcium transient duration at 90% decay, Tmax time to peak contraction, TR75 time to 75% relaxation (measured as the time interval from Tmax to the time of 75% relaxation). ENDO (E) model differs from the EPI model in the electrical properties only; ENDO (E ? Ca) model accounts for the heterogeneity in the electrical and calcium handling properties between the cells; ENDO is a complete ENDO model accounting for the heterogeneity in the electrical, calcium handling and mechanical properties; ENDO ( -MEF) model is a complete ENDO model excluding mechano-calcium-electric feedback (MEF) Table 4 Characteristics of action potential and low-loaded shortening in the ENDO and EPI models, compared to experimental data [ 23] APD90 exp [34 ± 1 /C176C], ms APD90 mod, ms Tmax exp [34 ± 1 /C176C], ms Tmax mod, ms TR75 exp [34 ± 1 /C176C], ms TR75 mod, ms ENDO 272 ± 11 [23] 287 158 ± 5 [23] 186 119 ± 8 [23] 102 EPI 221 ± 8 [23] 217 126 ± 3 [23] 133 91 ± 5 [23] 71 D 1.23 1.32 1.25 1.40 1.31 1.44 Delta ( D) shows the ENDO/EPI ratio between characteristics (mean values are used for the experimental data). In square brackets, the temperature at which experiments were carried out is shown 396 J Physiol Sci (2018) 68:387–413 123 feedback in our models [ENDO (–MEF), Table 3]. In this case, delayed Tmax and TR 75 had no effect on CaT 90 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 cal- cium handling parameters is most essential for the dis- tinctions in Ca 2? 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. 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). Ca21 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 2? transient (Table 3) that qualita- tively corresponded to experimental recordings in guinea pig ENDO and EPI cardiomyocytes ( *10% difference in Ca2? transient duration at 90% decay (CaT 90)[ 3]) and canine transmural wedges ( *20% in CaT 90 [3, 13, 38]). In accordance with the experimental data [ 3], the dias- tolic [Ca 2?]i was similar in the ENDO and EPI models (9.18 and 9.48 nM under isometric conditions, respec- tively). The lower maximal velocity of the SR pump led to a lower SR Ca 2? 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 cardiomy- ocytes 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 2?]i (1.78 vs 1.68 lM, respectively). Experi- mental data on the differences in peak [Ca 2?]i are contro- versial. Wan and co-workers demonstrated a lower amplitude of the Ca 2? 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 2? 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 2?]i between the cell subtypes, despite lower expression of SERCA2 in the subendocardium [ 13, 47, 48]. We determined the mechanisms underlying the differ- ences in Ca 2? transients between ENDO and EPI cells by analyzing the amount of Ca 2? (D[Ca2?]i) that enters into and exits from the cytosol over the duration of the calcium transient ( &500 ms) during isometric twitches (Fig. 6). A greater amount of Ca 2? (?34.4%, 9.4 lM in the ENDO model vs 7.0 lM in the EPI model) entered the ENDO cell via iCaL, while a lesser amount of Ca 2? entered the cytosol via iNaCa(r) (-34.9%, 1.7 vs 2.6 lM) and a lesser amount of Ca 2? extruded from the cytosol via iNaCa(f) (?2.2%, 17.6 vs 18.0 lM) in the ENDO cell. Thus, in total, about the same amount of Ca 2? entered the ENDO and EPI cells via transmembrane currents over the duration of cal- cium transient. Despite a smaller amplitude of irel in the ENDO cell, the Ca 2? release from the SR was longer lasting, resulting in the release of a greater amount of Ca 2? from the SR ( ?8.6%, 84.7 vs 78.0 lM). Therefore, the latter mechanism provided for almost identical Ca 2? amplitude, despite the lower SR Ca 2? content, and con- tributed 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 2? transient contributed to a delayed time Fig. 6 Amount of Ca 2? (D[Ca2?]i, lM) 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 2? are shown: iCaL Ca2? current via L-type Ca 2? channels, ibCa background Ca2? current, iNaCa(f) NCX current during the forward mode, iNaCa(r) NCX current during the reverse mode, irel Ca2? release flow from the SR, iup SERCA pump flow from cytosol to the SR J Physiol Sci (2018) 68:387–413 397 123 to peak contraction ( Tmax) and a slower relaxation of con- traction (longer period of time from peak to 75% of relax- ation, TR 75) 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% dis- crepancy between Tmax and about 20% between TR 75 in the cells. Transmural differences in the inherent parameters of intracellular Ca 2? handling essentially increased discrep- ancy between Tmax and TR 75 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 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 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 Fmax) at isotonic mode of contractions. Tension ( F) developed by cells ( top panel ) and cell deformations (bottom panel , DL in % of the initial length L0) during the cardiac cycle in ENDO and EPI models. The stimulus was applied at 60 ms 398 J Physiol Sci (2018) 68:387–413 123 [see ENDO (E ? Ca ? M) in Table 3; Fig. 7], which, via mechano-dependence of Ca 2? kinetics, modulate the Ca 2? transient and Ca2? 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 con- traction 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). Our EPI and ENDO models produced similar contrac- tion amplitudes at an initial SL o = 2.0 lm in both high- loaded (isometric) and low-loaded isotonic contractions at different afterloads (Fig. 8). In both cell types, a decrease in the afterload resulted in a prolongation of Tmax (ENDO: 20% increase in the Tmax during the afterloaded contraction at a low load of 20% F max compared to Tmax 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 90 (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 multicel- lular myocardial preparations [ 49–51]. The mechanisms underlying the load dependence of both T max and TR75 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-depen- dence of their kinetics, which resulted in a prolongation of Ca 2? 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 proper- ties of the cells via mechano-calcium-electric feedback.

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 mech- anisms 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 quantita- tively to available experimental data specific to guinea pig (Table 1). Unlike other models that have described trans- mural heterogeneity in the electrical properties of the myo- cardium [25–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 2? transients and the time course of contraction. The use of different parameters of Ca 2? handling and myofilament mechanics (Xb cycling) for ENDO and EPI models allowed us to reproduce a longer duration of the Ca 2? 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 con- tributes essentially to the distinctions in the time course of the Ca2? transient and, via mechano-calcium-electric feed- back, 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 differ- ences in ion current expression, particularly potassium cur- rents [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, ito is absent [56]. Our mod- eling results suggested that the higher iCaL in guinea pig ENDO cardiomyocytes mainly underlies the longer APD and the longer contraction in ENDO cells (Table 2). Ca21 handling heterogeneity A wide range of studies have described various aspects of spatial heterogeneity in the electrophysiological properties of J Physiol Sci (2018) 68:387–413 399 123 cardiomyocytes and underlying ionic mechanisms. By con- trast, less is known about the heterogeneity in the intracellular Ca2? 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 differ- ences in Ca 2? transients between ENDO and EPI myocytes [3, 18, 26, 38]. Some authors have suggested that the hetero- geneity in Ca 2? 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 2? release from the SR in the ENDO region of the guinea pig [3]. We showed that changes only in the electrophysiologi- cal parameters did not allow the ENDO model to reproduce significantly slower Ca 2? transient against the EPI model (Table 3) as observed experimentally in isolated car- diomyocytes [ 3, 13, 38]. Our models predict that despite the lower SR Ca 2? content in the ENDO cells, a greater amount of Ca 2? released from the SR (Fig. 6) in the cells may provide for similar Ca2? amplitudes in both cell types, which together with slower uptake of Ca2? 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 myo- cardium 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 mechanicalproperties 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 Fr ank–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 theslopes of the ‘ ‘cell length-ten- sion’ ’ 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 ventric- ular 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 2? 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 dif- ferences 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 Tmax and TR 75 between the cells of 13 and 22%, respec- tively. Note that the whole set of ENDO electrical parameters produced less effect on contraction-relaxation (6% for T max and 13% for TR 75) than that of some indi- vidual currents, e.g., iCaL and iKs. However, neither of the changes in the electrical parameters produced the roughly 25–30% difference in contraction characteristics as regis- tered experimentally (Table 4). Then, we compared the relative contribution of the differences in ionic currents with that combined with dif- ferences in calcium dynamics, and then additionally with the distinctions in the mechanical properties as well (Table 3). We revealed that distinctions in calcium han- dling and mechanical parameters slightly enhance the transmural gradient in AP characteristics between the ENDO and EPI cells (Table 3) and produce essential gra- dients in Ca 2? transient and contraction between the cells. Then, by excluding mechano-calcium-electric feedback in the ENDO model [Table 3, ENDO ( -MEF)], we demon- strated explicitly that these feedback mechanisms underlie a significant mechanical effect on the duration of Ca 2? transient in the ENDO cell, which slightly affects the 400 J Physiol Sci (2018) 68:387–413 123 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 trans- mural differences in AP generation [ 25–27] referred to in our work, the incorporation would have led to a higher gradient in Ca 2? 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 electro- physiology 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 repro- duce the combination of experimental data on the gradients in the profiles of AP, Ca 2? transient and contraction observed in isolated ENDO and EPI cells, and to explain underlying mechanisms. There is the lack of experimental data on transmural dif- ferences in mechanical function of ENDO and EPI car- diomyocytes 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 dif- ferently 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 electrome- chanical properties of cells. Under normal physiologic con- ditions, heterogeneous myocardial elements confer system homogeneity to optimize cardiac function [ 53, 59]. Con- versely, 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.

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 /C176C( T a b l e1, Table 4). However, some character- istics, such as ‘SL-passive tension’ (Fig. 4), were specified at room temperature. Temperature can possibly modulate transmural differences in the electrical and mechanical func- tions 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 2? activation of contractile proteins, and have not accounted for mechano- sensitive ionic channels and mechano-dependent SR Ca 2? 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 experi- mental 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 iNaK current (to repro- duce experimental data on transmural differences in iNaCa) 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 bor- rowing 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 iNaCa for a guinea pig model. In particular, it is known that the Na ? channel isoform Na v 1.5 is a major contributor to ipNa [1, 62], and NCX1 is a major contributor to iNaCa [36, 63]i n both guinea pig and canine ventricular cardiomyocytes, suggesting that the role of ipNa and iNaCa 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 2?,o r 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.).

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 J Physiol Sci (2018) 68:387–413 401 123 modeling results predict that heterogeneity in the parame- ters 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 2? 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. Acknowledgments This work was supported by RF Government Resolution #211 of March 16, 2013 and Program of the RAS Pre- sidium #I.33 G. Author contributions AK, NB-V: conception of the mathematical models, computational simulations, design, analysis and interpreta- tion of the computational experiments. LK, OS: conception of the mathematical models, design, analysis and interpretation of the computational experiments. GI: analysis and interpretation of the computational experiments. The manuscript was written by AK and OS, with the assistance of NB-V, LK, GI. All authors approved the final version of the manuscript. Compliance with ethical standards Conflict of interest The authors declare that they have no conflicts of interest related to this study. Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors. 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, iNa iNa ¼ gNa /C1 m3 /C1 h /C1ð V /C0 EmhÞ; dm dt ¼ am /C1ð 1 /C0 mÞ/C0 bm /C1 m; dh dt ¼ ah /C1ð 1 /C0 hÞ/C0 bh /C1 h; am ¼ 2000 if jV þ 41j\dm; 200 /C1ð V þ 41Þ 1 /C0 e/C0 0:1/C1ð Vþ41Þ otherwise ( bm ¼ 8000 /C1 e/C0 0:056/C1ð Vþ66Þ; ah ¼ 20 /C1 e/C0 0:125/C1ð Vþ75Þ; bh ¼ 2000 1 þ 320 /C1 e/C0 0:1/C1ð Vþ75Þ : Persistent Na? current, ipNa ipNa ¼ gpNa 1 þ e /C0ð Vþ56:51Þ 6:62 /C1ð V /C0 ENaÞ: L-type Ca2? current, iCaL iCaL ¼ d /C1 f /C1 f2;ds /C1ð iCaLCa þ iCaLNa þ iCaLKÞ; iCaLCa ¼ PCa /C1 4 /C1 F R /C1 T 1 1 /C0 e /C0 2/C1ð V/C0 50Þ/C1 F R/C1 T ½Ca2þ/C138 i /C1 e 100/C1 F R/C1 T /C0½ Ca2þ/C138 o /C1 e /C0 2/C1ð V/C0 50Þ/C1 F R/C1 T /C16/C17 /C1ð V /C0 50Þ; iCaLNa ¼ 0:01 /C1 PCa /C1 F R /C1 T 1 1 /C0 e /C0ð V/C0 50Þ/C1 F R/C1 T ½Naþ/C138 i /C1 e 50/C1 F R/C1 T /C0½ Naþ/C138 o /C1 e /C0ð V/C0 50Þ/C1 F R/C1 T /C16/C17 /C1ð V /C0 50Þ; iCaLK ¼ 0:002 /C1 PCa /C1 F R /C1 T 1 1 /C0 e /C0ð V/C0 50Þ/C1 F R/C1 T ½Kþ/C138 i /C1 e 50/C1 F R/C1 T /C0½ Kþ/C138 o /C1 e /C0ð V/C0 50Þ/C1 F R/C1 T /C16/C17 /C1ð V /C0 50Þ; dd dt ¼ ad /C1ð 1 /C0 dÞ/C0 bd /C1 d; df dt ¼ af /C1ð 1 /C0 f Þ/C0 bf /C1 f ; E0d ¼ V þ 19; E0f ¼ V þ 34; ad ¼ 360 if jE0dj\0:0001; 90 /C1 E0d 1 /C0 e/C0 E0d 4 otherwise ; 8 < : bd ¼ 360 if jE0dj\0:0001; /C0 36 /C1 E0d 1 /C0 e E0d 10 otherwise ; 8 < : af ¼ 7:5i f jE0f j\0:0001; /C0 1:875 /C1 E0f 1 /C0 e E0f 4 otherwise ; 8 < : bf ¼ 3:6 1 þ e /C0 E0f 4 ; 402 J Physiol Sci (2018) 68:387–413 123 df2;ds dt ¼ 20 /C1 1 /C0 ½Ca2þ/C138 ds 0:001 þ½ Ca2þ/C138 ds þ f2;ds /C18/C19/C18/C19 : Na?–Ca2? exchange current, i NaCa iNaCacmpt ¼ Fraccmpt /C1 iNaCa max /C1 e 0:37/C1 V/C1 F R/C1 T /C1½ Naþ/C138 3 i /C1½ Ca2þ/C138 o /C0 e /C0 0:63/C1 V/C1 F R/C1 T /C1½ Naþ/C138 3 o /C1½ Ca2þ/C138 cmpt 1 þ ½Ca2þ /C138 cmpt 0:0069 ; where cmpt is used to denote either cytosole (cyt) or dyadic space compartments (ds) of the cell; Fraccyt ¼ 0:999; Fracds ¼ 0:001; iNaCa ¼ iNaCacyt þ iNaCads : Na?–K? pump, iNaK iNaK ¼ iNaK max /C1 ½Kþ/C138 o 1 þ½ Kþ/C138 o /C1 ½Naþ/C138 i 30 þ½ Naþ/C138 i : Rapid delayed rectifier potassium current, i Kr iKr ¼ xr /C1 gKr /C1 ffiffiffiffiffiffiffi K0 5:4 r /C1 V /C0 EK 1 þ e V/C0 VoKr kKr ; axr ¼ 10 1 þ e /C0 V/C0 10ðÞ 370 ; bxr ¼ AKr /C1 e /C0 V/C0 20ðÞ k1Kr ; dxr dt ¼ axr /C1 1 /C0 xrðÞ /C0 bxr /C1 xr: Slow delayed rectifier potassium current, i Ks iKs ¼ gKsffiffiffiffiffiffiffiffiffi ½K/C138 0 p /C1 x2 s /C1 V /C0 EKsðÞ ; axs ¼ 14 1 þ AKs /C1 e /C0 V/C0 68ðÞ BKs ; bxs ¼ e /C0 V 88 ; dxs dt ¼ axs 1 /C0 xsðÞ /C0 bxsxs: Inward rectifier K ? current, iK1 iK1 ¼ gK1 /C1 ½Kþ/C138 o ½Kþ/C138 o þ 10 /C1 ðV /C0 EKÞ 1 þ e 1:25/C1 V/C0 EK /C0 10ðÞ /C1 F R/C1 T : Transient outward K ? current, ito ito ¼ gto /C1 s /C1 r /C1ð V /C0 EKÞ; dr dt ¼ 333 /C1 1 1 þ e /C0 Vþ4ðÞ 5 /C0 r ! ; as ¼ 0:033 /C1 e /C0 V 17 ; bs ¼ 33 1 þ e/C0 0:125/C1 Vþ10ðÞ ; ds dt ¼ as /C1 1 /C0 sðÞ /C0 bs /C1 s:

Na? current, ibNa ibNa ¼ gbNa /C1ð V /C0 ENaÞ:

Ca2? current, ibCa ibCa ¼ gbCa /C1ð V /C0 ECaÞ: Reversal potentials, E Na,E K,E Ks,E Ca,E mh ENa ¼ R /C1 T F /C1 ln ½Naþ/C138 o ½Naþ/C138 i ; EK ¼ R /C1 T F /C1 ln ½Kþ/C138 o ½Kþ/C138 i ; EKs ¼ R /C1 T F /C1 ln ½Kþ/C138 o þ PKNa /C1½ Naþ/C138 o ½Kþ/C138 i þ PKNa /C1½ Naþ/C138 i ; ECa ¼ R /C1 T 2 /C1 F /C1 ln ½Ca2þ/C138 o ½Ca2þ/C138 i ; Emh ¼ R /C1 T F /C1 ln ½Naþ/C138 o þ 0:12 /C1½ Kþ/C138 o ½Naþ/C138 i þ 0:12 /C1½ Kþ/C138 i : Intracellular and extracellular Na ? and K? concentrations, [Na ?]i,[ K?]i,[ K?]o d½Naþ/C138 i dt ¼ /C0 1 VCell /C1 Vi ratio /C1 F /C1 iNa þ ipNa þ ibNa þ iCaLNa þ 3 /C1 iNaK þ 3 /C1 iNaCa /C0/C1 ; d½Kþ/C138 i dt ¼/C0 1 VCell /C1 Vi ratio /C1 F /C1 iK1 þ iKr þ iKs þ iCaLK þ ito /C0 2 /C1 iNaKðÞ : d½Kþ/C138 o dt ¼ 1 VCell /C1 Ve ratio /C1 F /C1 iK1 þ iKr þ iKs þ iCaLK þ ito /C0 2 /C1 iNaKðÞ /C0 0:7 /C1½ Kþ/C138 o /C0½ Kþ/C138 b /C0/C1 : J Physiol Sci (2018) 68:387–413 403 123 Membrane potential, V dV dt ¼ /C0 1 Cm /C1 istim þ iNa þ ipNa þ ibNa þ iNaCa þ iCaL þ ibCa þ ito þ iKr þ iKs þ iK1 þ iNaK /C0 Þ; istim ¼/C0 3:0n A o na t6 0m s ; off at 62 :5m sðÞ : Calcium handling block Ca2? release from the SR, i rel irel ¼ 10000 /C1 FrAct FrAct þ 0:25 /C18/C19 2 þ0:05 ! /C1½ Ca2þ/C138 rel; Kact ¼ 500 /C1 RegBindSite; Kinact ¼ 60 þ 500 /C1 RegBindSite; RegBindSite ¼ ½Ca2þ/C138 i ½Ca2þ/C138 i þ Km;Cacyt þ 1 /C0 ½Ca2þ/C138 i ½Ca2þ/C138 i þ Km;Cacyt ! /C1 ½Ca2þ/C138 ds ½Ca2þ/C138 ds þ 0:01 Þ2; SpeedRel ¼ 5i f V\ /C0 50; 1 otherwise ; /C26 FrPrec ¼ 1 /C0 FrAct /C0 FrProd; dFrAct dt ¼ SpeedRel /C1ð FrPrec /C1 KAct /C0 FrAct /C1 KinactÞ; dFrProd dt ¼ SpeedRel /C1ð FrAct /C1 KInact /C0 FrProdÞ: Ca2? uptake by the SR-pump, i up K1 ¼ 0:00012; K2 ¼½ Ca2þ/C138 i þ½ Ca2þ/C138 up /C1 K1 þ 0:00021; iup ¼ aup /C1 ½Ca2þ/C138 i K2 /C1 1 þ ½Ca2þ /C138 up 4 /C16/C17 /C0 bup /C1 ½Ca2þ/C138 up /C1 K1 K2 : Ca2? kinetics within the SR, i trans itrans ¼ 15 /C1½ Ca2þ/C138 up /C0½ Ca2þ/C138 rel /C16/C17 ; Ca2? concentration within the SR compartments, [Ca2?]up, [Ca2?]rel d½Ca2þ/C138 up dt ¼ Vi ratio Vup ratio /C1 iup /C0 itrans; d½Ca2þ/C138 rel dt ¼ Vup ratio Vrel ratio /C1 itrans /C0 irel 1 þ 0:65/C1½ CaS/C138 tot ½Ca 2þ /C138 relþ0:65 : Ca2? concentration within the cytosol and the DS, [Ca2?]i, [Ca2?]ds d½Ca2þ/C138 i dt ¼ /C0 1 2 /C1 VCell /C1 Vi ratio /C1 F /C1 iCaLCacyt þ ibCa /C0 2 /C1 iNaCacyt /C0/C1 þ½ Ca2þ/C138 ds /C1 10 /C1 Vds ratio þ Vrel ratio Vi ratio /C1 irel /C0 d½CaTnC/C138 dt /C0 d½B1/C138 dt /C0 d½B2/C138 dt /C0 iup; d½Ca2þ/C138 ds dt ¼/C0 1 2 /C1 Vds ratio /C1 VCell /C1 Vi ratio /C1 F /C1ð iCaLCads /C0 2 /C1 iNaCads Þ/C0 10 /C1½ Ca2þ/C138 ds: Kinetics of CaTnC complexes, [CaTnC] d½CaTnC/C138 dt ¼ kon /C1½ TnC/C138 tot /C0½ CaTnC/C138 /C0/C1 /C1½ Ca2þ/C138 i /C0 koff /C1 e/C0 kA /C1½ CaTnC/C138 /C1 pNA /C1½ CaTnC/C138 : A key point of this equation is the cooperativity of 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 pNA expresses this cooperativity in the above equation, and its particular form may be written as follows: pNA ¼ 1i f NA /C20 0; 0:03NA if 0 \NA /C20 1; 0:03 otherwise ; 8 < : where NA means an average fraction of the attached Xbs falling on one CaTnC complex; i.e. NA ¼ N ½CaTnC/C138 . In accordance with CaTnC–CaTnC cooperativity (Fig. 2), we consider CaTnC off-rate to be decreasing with increases in [CaTnC]. Dependence e/C0 kA /C1½ CaTnC/C138 defines this cooperativity in the above equation, where kA is the model parameter (B1, Table 4). Ca2? buffering system, B 1,B 2 d½B1/C138 dt ¼ b1on /C1½ B1/C138 tot /C0½ B1/C138 /C0/C1 /C1½ Ca2þ/C138 i /C0 b1off /C1½ B1/C138 ; d½B2/C138 dt ¼ b2on /C1½ B2/C138 tot /C0½ B2/C138 /C0/C1 /C1½ Ca2þ/C138 i /C0 b2off /C1½ B2/C138 : 404 J Physiol Sci (2018) 68:387–413 123 Mechanical block of the model Kinetics of force generating Xbs, N dN dt ¼ kpv /C1 Mð½CaTnC/C138Þ /C1 n1ðl1Þ/C1 Loz /C1ð 1 /C0 NÞ/C0 kmv /C1 N; Mð½CaTnC/C138Þ ¼ ½CaTnC/C138 ½TnC/C138 tot /C16/C17 l /C1 1þ0:6lðÞ ½CaTnC/C138 ½TnC/C138 tot /C16/C17 l þ 0:6l means RU end-to-end interaction between adjacent tropomyosin segments in the case where both of them are affected by the formation of respective CaTnC complexes (Eq. 5, Fig. 2), where l is the cooperativity degree (B1., Table 4). n1ðl1Þ¼ 0i f 0 :6 /C1 l1 þ g2 /C20 0; 0:6 /C1 l1 þ g2 if 0 \0:6 /C1 l1 þ g2\1; 1 otherwise : 8 < : is the probability of a myosin cross-bridge ‘ ‘finding’ ’ a vacant site on the actin filament; in other words, the dependence of n1 on l1 means sarcomere lattice spacing. Lozðl1Þ¼ l1 þ 1:14 1:6 if l1 /C20 0:55; 1:69 1:6 otherwise: 8 >: is a normalized linear dependence of the sarcomere overlap zone on the SL; kpv and kmv are velocity-dependent rate constants of cross-bridge attachment and detachment, respectively: kpv ¼ 0:6822 /C1 v0 /C1 qv /C1 G/C3 ; kmv ¼ v0 /C1 qv /C1ð 1 /C0 0:6822 /C1 G/C3 Þ; qv ¼ 17 /C0 259 /C1 v vmax if v /C20 0; /C0 2:3 /C1 v 0:79 /C1 vmax þ 17:3i f 0 \v /C20 0:79 /C1 vmaxðÞ ; 15 1 þ 5 v /C0 0:79 /C1 vmaxðÞ vmax /C18/C19 10 otherwise 8 >> > >>> > > >>> > > > > : Mechanical variables, F CE,F VS,F SE,F PE,F XSE,l ,l 1,l ex The following formula gives the force developed by the activated contractile element CE in the model (Fig. 1): FCE ¼ k /C1 N /C1 pðvÞ; Function pðvÞ means the dependence of the average Xb force on the sarcomere shortening/lengthening velocity ( v): pðvÞ¼ P/C3 ðvÞ G/C3 ðvÞ ; where P/C3 ðvÞ is the dependence of the steady-state sar- comere force on the shortening/lengthening velocity, and G /C3 ðvÞ is the dependence of the steady-state sarcomere stiffness on the velocity. Both functions are normalized so that P/C3 ð0Þ¼ 1, G/C3 ð0Þ¼ 1. P/C3 ¼ 0:25 /C1 1 þ v vmax /C18/C19 0:25 /C0 v vmax if v /C20 0; 1:5 /C0 0:0625 25 /C1 v vmax /C18/C19 2 þ 1:25 /C1 v vmax þ 0:125 otherwise; 8 >> > >>> > > >>> > > > > : G/C3 ¼ 1 þ 0:6 /C1 v vmax if /C0 vmax /C20 vðÞ and v /C20 0ðÞ ; P/C3 1:1 /C1 v 0:25 /C1 vmax þ 1 if 0 \v /C20 0:1 /C1 vmaxðÞ ; P/C3 /C1 e/C0 v/C0 0:1ðÞ 4 1:1 /C1 v 0:25 /C1 vmax þ 1 otherwise: 8 >> >>> > >>> > >>> > >>> > >: The force generated by the damper VS during cell deformation (Fig. 1): F VS ¼ bv /C1 e16/C1 lðtÞ /C1 dl dt ðtÞ: The following equations define the SE and PE elements of the rheological schemes underlying the mechanical module of the model: FSE ¼ b1 /C1ð e25/C1ð lðtÞ/C0 l1 ðtÞÞ /C0 1Þ; FPE ¼ b2 /C1ð ea2 /C1 lðtÞ /C0 1Þ: The following equations define the force F that is developed by the cardiomyocyte according to the rheo- logical scheme: F CE ¼ FSE; F ¼ FCE þ FPE þ FVS: The force of extra-series element XSE is defined as: FXSE ¼ b3 /C1ð e48/C1 lex ðtÞ /C0 1Þ: Thus, at each t, the functions l(t), l1(t) and v(t) are solutions of the system of differential–algebraic equations: k /C1 N /C1 pðvÞ¼ b1 /C1ð e25/C1ð l/C0 l1 Þ /C0 1Þ; dl dt ¼ b3 /C1ð e48/C1 lex /C0 1Þ/C0 b1 /C1ð e25/C1ð l/C0 l1 Þ /C0 1Þ/C0 b2 /C1ð ea2 /C1 l /C0 1Þ bv /C1 e16/C1 l ; J Physiol Sci (2018) 68:387–413 405 123 dl1 dt ¼ v: For the single cardiomyocyte model, these equations are complemented with the isometric condition: l þ lex /C17 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 Parameter Definition ENDO EPI R [mJ M-1 K-1] Gas constant 8314.472 F [C M-1] Faraday's constant 96485.341 T [K] Temperature 310 Cm [μF] membrane capacitance 12.0 10-5 RCell [μm] cell radius 12 LCell [μm] cell length 74 VCell [μm3] cell volume 10-9π⋅RCell⋅LCell Vi_ratio volume fraction of cytosolic space 0.49 Ve_ratio volume fraction of extracellular space 0.4 Vds_ratio volume fraction of dyadic space 0.1 Vrel_ratio volume fraction of SR release site 0.003 Vup_ratio volume fraction of SR uptake site 0.03 [Ca2+]o [mM] extracellular Ca2+ concentration 2 [Na+]o [mM] extracellular Na+ concentration 140 [K+]b [mM] K+ concentration in the space adjacent to the cell 4 406 J Physiol Sci (2018) 68:387–413 123 Table 6 Parameters for the electrical block Parameter Definition ENDO EPI gNa [μS] maximum membraneiNa conductance 2.5 gK1 [μS] maximum membraneiK1 conductance 0.5 gto [μS] maximum membraneito conductance 0.0045 iNaCa_max [nA] maximum iNaCa 0.00017 gbNa [μS] maximum membraneibNa conductance 0.0006 gKr [μS] maximum membraneiKr conductance 0.00157 0.00183 VoKr [mV] potential at 50% activation of iKr -17 -20 kKr [mV] slope coefficient 84.2 83.0 k1Kr [mV] slope coefficient 7.0 7.2 AKr [s-1] parameter of iKr inactivation 0.008 0.004 gKs [μS] maximum membraneiKr conductance 0.0048 0.0073 AKs [s-1] parameter of iKs activation 1.09 1.00 BKs [s-1] parameter of iKs activation 24 28 iNaK_max [nA] maximum iNaK 0.8 0.7 gpNa [μS] maximum membraneipNa conductance 0.00062 0.00054 PCa [nA mM] membrane permeability for Ca2+ through L-type channels 0.16 0.10 gbCa [μS] maximum membraneibCa conductance 0.00027 0.00025 J Physiol Sci (2018) 68:387–413 407 123 Table 7 Parameters for the Ca 2? handling Parameter Definition ENDO EPI [CaS]tot [mM] total calsequestrin concentration 40 [TnC]tot [mM] total CaTnC concentration 0.07 [B1]tot [mM] total B1 buffer concentration 0.08 [B2]tot [mM] total B2 buffer concentration 0.1 kon [mM-1 s-1] rate constant 70000 koff [s-1] rate constant 200 kA [mM-1] CaTnC-CaTnC cooperativity degree 30 b1on [mM-1 s-1] rate constant 100000 b1off [s-1] rate constant 100000 b2on [mM-1 s-1] rate constant 1000 b2off [s-1] rate constant 3 βup [mM s-1] SERCA Ca2+ reverse flux rate 0.03 αup [mM s-1] SERCA Ca2+ uptake flux rate 0.8 1 Km,CaCyt [mM] Sensitivity of the binding sites of RyR to the cytosolic Ca2+ 0.001 0.0005 408 J Physiol Sci (2018) 68:387–413 123 Parameter Definition ENDO EPI β1[mN mm-2] scaling coefficient for SE force 0.093 β3 [mN mm-2] scaling coefficient for XSE force 0.015 βv [mN s mm-3] scaling coefficient for VS force 0.00016 λ[mN mm-2] scaling coefficient for CE force 40 α2 [μm-1] passive tension–cell length curve parameter 2.92 β2 [mN mm-2] scaling coefficient of PE force 0.2606 0.0838 vmax [μm s-1] maximum velocity of unloaded shortening 5.0 7.0 g2 coefficient of probability for a myosin attachment on an actin filament 0.40 0.52 χ0 Xb attachment- detachment ratio coefficient 3.0 4.0 μ RU end-to-end cooperativity degree 3.0 3.3 All the initial values for the ENDO and EPI pha se variables are shown in their steady state, achieved by allowing the models to run for 100 seconds at 1 Hz pacing rate. Phase Variable Definition ENDO EPI l1 [μm] deformation of CE relative to the slack length 0.2139 0.2262 l [μm] deformation of PE relative to the slack 0.2139 0.2262 length lex [μm] deformation of XSE relative to the slack length 0.0436 0.0313 v [μm/s] sarcomere shortening velocity 6.7E-6 0.38E-6 N fraction of force- generating Xbs 3.3E-7 0.52E-7 Table 8 Parameters for the mechanical block J Physiol Sci (2018) 68:387–413 409 123 Phase Variable Definition ENDO EPI FrAct fraction of activated channels 0.00056 0.00205 FrProd fraction of inactivated channels 0.22573 0.30025 [Ca2+]up [mM] Ca2+ concentration in longitudinal reticulum 1.38 1.50 [Ca2+]rel [mM] Ca2+ concentration in terminal cisterns 1.37 1.49 [Ca2+]i [mM] cytosolic Ca2+ concentration 9.2E-6 9.5E-6 [Ca2+]ds [mM] DS Ca2+ concentration 1.6E-6 2.0E-6 [CaTnC][ m M ] CaTnC concentration 2.3E-4 2.3E-4 [B1][ m M ] B1 buffer concentration 4.0E-4 4.1E-4 [B2][ m M ] B2 buffer concentration 9.0E-4 7.2E-4 Phase Variable Definition ENDO EPI xr coefficient of open probability of the potassium channels 0.51E-4 1.61E-4 xs coefficient of open probability of the potassium channels 0.0659 0.0598 r potassium channel activation coefficient 1.78E-8 1.80E-8 s potassium channel inactivation coefficient 0.995 0.996 m sodium channel activation coefficient 0.0015 0.0016 h sodium channel inactivation coefficient 0.9948 0.9948 d L-channel activation coefficient 2.2E-8 2.2E-8 f L-channel inactivation coefficient 0.9999 0.9999 f2,ds degree of activation of the L-channels in DS 0.9983 0.9980 [Na+]i [mM] intracellular Na+ concentration 4.6423 4.9982 [K+]i [mM] intracellular K+ concentration 138.0930 137.7124 [K+]o [mM] extracellular K+ concentration 3.9997 4.0019 V [mV] membrane potential -93.208 -93.142 Table 8 continued 410 J Physiol Sci (2018) 68:387–413 123

1. Osadchii OE, Soltysinska E, Olesen SP (2011) Na ? channel distribution and electrophysiological heterogeneities in guinea pig ventricular wall. Am J Physiol Heart Circ Physiol 300(3):H989–H1002 2. Nerbonne JM, Guo W (2002) Heterogeneous expression of voltage-gated potassium channels in the heart: roles in normal excitation and arrhythmias. J Cardiovasc Electrophysiol 13(4):406–409 3. Wan X et al (2005) Molecular correlates of repolarization alter- nans in cardiac myocytes. J Mol Cell Cardiol 39(3):419–428 4. Antzelevitch C, Fish J (2001) Electrical heterogeneity within the ventricular wall. Basic Res Cardiol 96(6):517–527 5. Cazorla O et al (2000) Differential expression of cardiac titin isoforms and modulation of cellular stiffness. Circ Res 86(1):59–67 6. Litten RZ et al (1985) Heterogeneity of myosin isozyme content of rabbit heart. Circ Res 57(3):406–414 7. Stelzer JE et al (2008) Transmural variation in myosin heavy chain isoform expression modulates the timing of myocardial force generation in porcine left ventricle. J Physiol 586(Pt 21):5203–5214 8. Ait Mou Y et al (2008) Differential contribution of cardiac sar- comeric proteins in the myofibrillar force response to stretch. Pflugers Arch 457(1):25–36 9. Cazorla O et al (1997) Resting tension participates in the mod- ulation of active tension in isolated guinea pig ventricular myo- cytes. J Mol Cell Cardiol 29(6):1629–1637 10. Solovyova O et al (2003) Mechanical interaction of heteroge- neous cardiac muscle segments in silico: effects on Ca2 ? han- dling and action potential. Int J Bifurc Chaos 13(12):3757–3782 11. Bryant SM et al (1998) Regional differences in the delayed rectifier current (IKr and IKs) contribute to the differences in action potential duration in basal left ventricular myocytes in guinea-pig. Cardiovasc Res 40(2):322–331 12. Cordeiro JM et al (2004) Transmural heterogeneity of calcium activity and mechanical function in the canine left ventricle. Am J Physiol Heart Circ Physiol 286(4):H1471–H1479 13. Lou Q et al (2011) Transmural heterogeneity and remodeling of ventricular excitation–contraction coupling in human heart fail- ure. Circulation 123(17):1881–1890 14. Derumeaux G et al (2000) Assessment of nonuniformity of transmural myocardial velocities by color-coded tissue Doppler imaging: characterization of normal, ischemic, and stunned myocardium. Circulation 101(12):1390–1395 Table 9 Phase variables Phase variable Definition ENDO EPI xr Coefficient of open probability of the potassium channels for i kr 0.51E-4 1.61E -4 xs Coefficient of open probability of the potassium channels for i ks 0.0659 0.0598 r Potassium channel activation coefficient for i to 1.78E-8 1.80E -8 s Potassium channel inactivation coefficient for i to 0.995 0.996 m Sodium channel activation coefficient 0.0015 0.0016 h Sodium channel inactivation coefficient 0.9948 0.9948 d L-channel activation coefficient 2.2E -8 2.2E -8 f L-channel inactivation coefficient 0.9999 0.9999 f2,ds Degree of activation of the L-channels in DS 0.9983 0.9980 [Na?]i (mM) Intracellular Na ? concentration 4.6423 4.9982 [K?]i (mM) Intracellular K ? concentration 138.0930 137.7124 [K?]o (mM) Extracellular K ? concentration 3.9997 4.0019 V (mV) Membrane potential -93.208 -93.142 FrAct Fraction of activated channels 0.00056 0.00205 FrProd Fraction of inactivated channels 0.22573 0.30025 [Ca2?]up (mM) Ca 2? concentration in longitudinal reticulum 1.38 1.50 [Ca2?]rel (mM) Ca 2? concentration in terminal cisterns 1.37 1.49 [Ca2?]i (mM) Cytosolic Ca 2? concentration 9.2E -6 9.5E -6 [Ca2?]ds (mM) DS Ca 2? concentration 1.6E -6 2.0E -6 [CaTnC] (mM) CaTnC concentration 2.3E -4 2.3E -4 [B1] (mM) B1 buffer concentration 4.0E -4 4.1E -4 [B2] (mM) B2 buffer concentration 9.0E -4 7.2E -4 l1 (lm) Deformation of CE relative to the slack length 0.2139 0.2262 l (lm) Deformation of PE relative to the slack length 0.2139 0.2262 lex (lm) Deformation of XSE relative to the slack length 0.0436 0.0313 v (lm/s) Sarcomere shortening velocity 6.7E -6 0.38E -6 N Fraction of force-generating Xbs 3.3E -7 0.52E -7 J Physiol Sci (2018) 68:387–413 411 123 15. Ashikaga H et al (2007) Transmural dispersion of myofiber mechanics: implications for electrical heterogeneity in vivo. J Am Coll Cardiol 49(8):909–916 16. Glukhov AV et al (2010) Transmural dispersion of repolarization in failing and nonfailing human ventricle. Circ Res 106(5):981–991 17. Clark RB et al (1993) Heterogeneity of action potential wave- forms and potassium currents in rat ventricle. Cardiovasc Res 27(10):1795–1799 18. Dilly KW et al (2006) Mechanisms underlying variations in excitation-contraction coupling across the mouse left ventricular free wall. J Physiol 572(Pt 1):227–241 19. Stankovicova T et al (2000) M cells and transmural heterogeneity of action potential configuration in myocytes from the left ven- tricular wall of the pig heart. Cardiovasc Res 45(4):952–960 20. Antzelevitch C (2010) M cells in the human heart. Circ Res 106(5):815–817 21. Quinn FR et al (2003) Myocardial infarction causes increased expression but decreased activity of the myocardial Na ?–Ca2? exchanger in the rabbit. J Physiol 553(Pt 1):229–242 22. Bell SP et al (2000) Alterations in the determinants of diastolic suction during pacing tachycardia. Circ Res 87(3):235–240 23. Wan X, Bryant SM, Hart G (2003) A topographical study of mechanical and electrical properties of single myocytes isolated from normal guinea-pig ventricular muscle. J Anat 202(6):525–536 24. ten Tusscher KH et al (2004) A model for human ventricular tissue. Am J Physiol Heart Circ Physiol 286(4):H1573–H1589 25. Pandit SV et al (2001) A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes. Biophys J 81(6):3029–3051 26. Bondarenko VE, Rasmusson RL (2010) Transmural heterogene- ity of repolarization and Ca 2? handling in a model of mouse ventricular tissue. Am J Physiol Heart Circ Physiol 299(2):H454– H469 27. Benson AP et al (2008) The canine virtual ventricular wall: a platform for dissecting pharmacological effects on propagation and arrhythmogenesis. Prog Biophys Mol Biol 96(1–3):187–208 28. Mullins PD, Bondarenko VE (2013) A mathematical model of the mouse ventricular myocyte contraction. PLoS One 8(5):e63141 29. Campbell SG et al (2008) Mechanisms of transmurally varying myocyte electromechanics in an integrated computational model. Philos Trans Ser A Math Phys Eng Sci 366(1879):3361–3380. https://www.ncbi.nlm.nih.gov/labs/articles/18593662/ 30. Kerckhoffs RC et al (2003) Homogeneity of cardiac contraction despite physiological asynchrony of depolarization: a model study. Ann Biomed Eng 31(5):536–547 31. Nickerson D, Smith N, Hunter P (2005) New developments in a strongly coupled cardiac electromechanical model. Europace 7(Suppl 2):118–127 32. Campbell SG et al (2009) Effect of transmurally heterogeneous myocyte excitation-contraction coupling on canine left ventricu- lar electromechanics. Exp Physiol 94(5):541–552 33. Gurev V et al (2010) Distribution of electromechanical delay in the heart: insights from a three-dimensional electromechanical model. Biophys J 99(3):745–754 34. Vasilyeva AD, Solovyova OE (2012) Electromechanical coupling in cardiomyocytes from transmural layers of guinea pig left ventricle. Biophysics 57(5):661–667 35. Vasilyeva A, Solovyova O (2012) Modeling of heterogeneity in electrical and mechanical properties of guinea pig ventricular myocytes. Comput Cardiol (CinC) 39:453–456 36. Gao J et al (2005) Transmural gradients in Na/K pump activity and [Na ?]I in canine ventricle. Biophys J 89(3):1700–1709 37. Bryant SM, Shipsey SJ, Hart G (1997) Regional differences in electrical and mechanical properties of myocytes from guinea-pig hearts with mild left ventricular hypertrophy. Cardiovasc Res 35(2):315–323 38. Laurita KR et al (2003) Transmural heterogeneity of calcium handling in canine. Circ Res 92(6):668–675 39. Sulman T et al (2008) Mathematical modeling of mechanically modulated rhythm disturbances in homogeneous and heteroge- neous myocardium with attenuated activity of Na ?–K? pump. Bull Math Biol 70(3):910–949 40. Noble D et al (1998) Improved guinea-pig ventricular cell model incorporating a diadic space, IKr and IKs, and length- and ten- sion-dependent processes. Can J Cardiol 14(1):123–134 41. Izakov V et al (1991) Cooperative effects due to calcium binding by troponin and their consequences for contraction and relaxation of cardiac muscle under various conditions of mechanical load- ing. Circ Res 69(5):1171–1184 42. Garny A et al (1895) Cellular open resource (COR): current status and future directions. Philos Trans Ser A Math Phys Eng Sci 2009(367):1885–1905 43. Zygmunt AC et al (2001) Larger late sodium conductance in M cells contributes to electrical heterogeneity in canine ventricle. Am J Physiol Heart Circ Physiol 281(2):H689–H697 44. Zygmunt AC, Goodrow RJ, Antzelevitch C (2000) I(NaCa) contributes to electrical heterogeneity within the canine ventricle. Am J Physiol Heart Circ Physiol 278(5):H1671–H1678 45. Banyasz T et al (2003) Endocardial versus epicardial differences in L-type calcium current in canine ventricular myocytes studied by action potential voltage clamp. Cardiovasc Res 58(1):66–75 46. Khokhlova A, Iribe G, Solovyova O (2015) Load-dependency in mechanical properties of subepicardial and subendocardial car- diomyocytes. Comput Cardiol 42:965–968 47. McIntosh MA, Cobbe SM, Smith GL (2000) Heterogeneous changes in action potential and intracellular Ca 2? in left ven- tricular myocyte sub-types from rabbits with heart failure. Car- diovasc Res 45(2):397–409 48. Cordeiro JM et al (2007) Cellular and subcellular alternans in the canine left ventricle. Am J Physiol Heart Circ Physiol 293(6):H3506–H3516 49. Lab MJ, Allen DG, Orchard CH (1984) The effects of shortening on myoplasmic calcium concentration and on the action potential in mammalian ventricular muscle. Circ Res 55(6):825–829 50. White E, Boyett MR, Orchard CH (1995) The effects of mechanical loading and changes of length on single guinea-pig ventricular myocytes. J Physiol 482(Pt 1):93–107 51. Iribe G, Helmes M, Kohl P (2007) Force-length relations in isolated intact cardiomyocytes subjected to dynamic changes in mechanical load. Am J Physiol Heart Circ Physiol 292(3):H1487–H1497 52. Katsnelson LB et al (2004) Influence of viscosity on myocardium mechanical activity: a mathematical model. J Theor Biol 230(3):385–405 53. Markhasin VS et al (2003) Mechano-electric interactions in heterogeneous myocardium: development of fundamental experimental and theoretical models. Prog Biophys Mol Biol 82(1–3):207–220 54. Katsnelson LB et al (2011) Contribution of mechanical factors to arrhythmogenesis in calcium overloaded cardiomyocytes: model predictions and experiments. Prog Biophys Mol Biol 107(1):81–89 55. Markhasin VS et al (2012) Slow force response and auto-regu- lation of contractility in heterogeneous myocardium. Prog Bio- phys Mol Biol 110(2):305–318 56. Ryder KO, Bryant SM, Hart G (1993) Membrane current changes in left ventricular myocytes isolated from guinea pigs after abdominal aortic coarctation. Cardiovasc Res 27(7):1278–1287 57. Ashikaga H et al (2009) Transmural myocardial mechanics dur- ing isovolumic contraction. JACC Cardiovasc Imaging 2(2):202–211 412 J Physiol Sci (2018) 68:387–413 123 58. Bollensdorff C, Lookin O, Kohl P (2011) Assessment of con- tractility in intact ventricular cardiomyocytes using the dimen- sionless ‘Frank-Starling gain’ index. Pflugers Arch 462(1):39–48 59. Katz AM, Katz PB (1989) Homogeneity out of heterogeneity. Circulation 79(3):712–717 60. Kohl P et al (2001) Sudden cardiac death by commotio cordis: role of mechano-electric feedback. Cardiovasc Res 50(2):280–289 61. Chung CS, Campbell KS (2013) Temperature and transmural region influence functional measurements in unloaded left ven- tricular cardiomyocytes. Physiol Rep 1(6):e00158 62. Maltsev VA et al (2008) Molecular identity of the late sodium current in adult dog cardiomyocytes identified by Nav1. 5 anti- sense inhibition. Am J Physiol Heart Circ Physiol 64(2):H667 63. Fujioka Y, Hiroe K, Matsuoka S (2000) Regulation kinetics of Na?–Ca2? exchange current in guinea-pig ventricular myocytes. J Physiol 529(3):611–623 J Physiol Sci (2018) 68:387–413 413 123