Mechano-electric feedback effects in a three-dimensional (3D) model of the contracting cardiac ventricle

Mechano-electric feedback affects the electrophysiological and mechanical function of the heart and the cellular, tissue, and organ properties. To determine the main factors that contribute to this effect, this study investigated the changes in the action potential characteristics of the ventricle during contraction. A model of stretch-activated channels was incorporated into a three-dimensional multiscale model of the contracting ventricle to assess the effect of different preload lengths on the electrophysiological behavior. The model describes the initiation and propagation of the electrical impulse, as well as the passive (stretch) and active (contraction) changes in the cardiac mechanics. Simulations were performed to quantify the relationship between the cellular activation and recovery patterns as well as the action potential durations at different preload lengths in normal and heart failure pathological conditions. The simulation results showed that heart failure significantly affected the excitation propagation parameters compared to normal condition. The results showed that the mechano-electrical feedback effects appear to be most important in failing hearts with low ejection fraction.


Introduction
Investigating the effects of intracellular electromechanical coupling and mechano-electric feedback (MEF) on myocardial function and its contribution to arrhythmogenesis is of great importance. Changes in the electrophysiological behavior explained by stretch-activated channels (SACs) were observed experimentally [1][2][3]. The impact of mechanical events on cardiac electrophysiology was studied in several cellular and tissue models [4][5][6][7] by including stretch activation of ion channels and mechanical modulation of cellular Ca 2+ handling as the major mechanisms of the MEF. Several strongly coupled electromechanical models of the ventricles described mechanical deformation triggered by electrical activation with consideration of the MEF [8][9][10][11]. Such studies typically demonstrate the effect of mechano-electric coupling via the SACs in terms of arrhythmogenesis without taking into account load induced electrophysiological changes. However, the alterations of ventricular loading conditions may provide a basis for initiation of arrhythmia as observed in experimental studies [12,13]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The purpose of this study was to investigate the effects of mechano-electric coupling via stretch-activated ion channels in the contracting ventricle. A previously developed threedimensional (3D) multiscale model of the contracting ventricle [14] was used to determine the influence of mechanical perturbations on electrophysiological parameters. Thus, an investigation of the effects of mechano-electric coupling was possible that could not be assessed experimentally. The model includes a considerable level of physiological detail and a noncontinuum method to represent the myocardium.
The electrophysiological effects of myocardial stretching were evaluated in the contracting ventricle model by varying the preload length and quantitatively analyzing the changes in action potential durations, activation, and recovery patterns. The model describes the initiation and propagation of the electrical impulse and the passive (stretch) and active (contraction) changes in cardiac mechanics. In addition, the model was used to evaluate the effects of different preload lengths on the cardiac electrophysiology under normal and heart failure (HF) pathological conditions.

Materials and methods
The simulation of MEF in this study consists of four major components: the cellular electrophysiology, force development, excitation propagation, and anatomical structure (Fig 1). These components were combined to form the 3D multiscale ventricular model for evaluating the influence of initial length of cardiac tissue (preload) on cardiac electrophysiology. The 3D model of the contracting ventricle was described earlier [14]. In this study, a brief overview of the contracting ventricle model was provided.

Electrophysiological behavior
The electrophysiological model was based on the ten Tusscher-Noble-Noble-Panfilov (TNNP) heterogeneous model of a human ventricular myocyte [15]. It described 12 ionic currents, variations in the intracellular concentration of Na + , K + , and Ca 2+ ions, and the transmembrane voltage. The time-dependent changes in the transmembrane voltage V m [mV] were described by the following differential equation: where t [ms] was the time, and C m [μF/cm 2 ] was the cell capacitance per unit surface area. The parameter I ion [μA/μF] was the sum of the membrane ionic currents through the ionic channels located on the sarcolemma, and I stim [μA/ μF] was the external stimulus current. In this study, thet TNNP model was extended by adding a stretch-activated channel (I SAC ) that corresponded to the stretch effects on the action potential (AP). The following sections describe the details. The sum of the transmembrane ionic currents in the electrophysiological model was described by Eq 2.
The fast Na + current was I Na . The Na + /K + pump current was I NaK , and the background Na + current was I bNa . The parameter I K1 was the inward rectifier K + current, and the transient outward current was I to . The rapid delayed rectifier current was I Kr , and the parameter I Ks was the slow delayed rectifier current. The plateau K + current was I pK , and the L-type Ca 2+ current was I CaL . The parameter I NaCa was the Na + /Ca 2+ exchanger, and the plateau Ca 2+ current was I pCa .
The background Ca 2+ leak current was I bCa ., For all simulations the TNNP model was applied with the parameter set of a human mid-myocardial ventricular cell. Eq 1 was solved using the Euler's method with a time step of Δt = 0.02 ms. This time step ensured a stable converging solution for the TNNP model equations.

Heart failure
To represent electrophysiological changes due to HF, the modifications proposed by Zlochiver were incorporated [16]. Briefly, HF condition in the model was introduced by varying maximum conductances of the transient outward current (I to ), the slow delayed rectifier current (I Ks ), the inward rectifier K + current (I K1 ), the fast Na + current (I Na ), the Na + /Ca 2+ exchanger (I NaCa ) and L-type Ca 2+ current (I CaL ) with values listed in S1 Table in the Supporting information.
In the model, the propagation of electrical excitation was represented by the reaction-diffusion equation for the isotropic medium, as described previously [14].

Mechano-electricfeedback
Stretch-activated channels. The effect of stretch-activated channels on the cellular electrophysiological state was implemented by integrating the formulation of the SACs into the electrophysiological model based on the formulation by Zabel et al. [2] with some modifications.
The parameter V m [mV] was the transmembrane voltage, and V rev [mV] was the SAC reversal potential. The parameter G SAC [nS/pF] was the maximum SAC conductance, the parameter K defined the amount of current when the cell was not stretched, and α was the parameter that controlled the sensitivity to stretch. The sarcomere length in Eq 3 was replaced by the stretch ratio λ that was calculated as the ratio between the current and the resting half sarcomere length (0.97 μm) obtained from the myofilament model [17]: All the parameter values discussed above are listed in Table 1. Due to the absence of consistent data on the SACs, a reversal potential V rev of -20 mV [2] and a maximum conductance G SAC of 0.025 nS/pF (a similar range was observed by Sachs et al. [18]) were chosen.

Force development model
The intracellular calcium released during an electrical activation coupled the electrical and mechanical components in the 3D model simulations. This [Ca 2+ ] i was used as an input parameter for the myofilament model of Negroni and Lascano (NL) [17] previously described [14]. The NL model was incorporated into the simulation to represent the force development and sarcomere shortening within each myocyte.

Fiber-based geometrical model
The fiber-based model of the contracting ventricle [19] formed the geometrical basis for the evaluation of the MEF effects. The 3D structure of the model and the modeling method of the cardiac mechanics were described in earlier research [14]. A twitch obtained from the sarcomere model of force generation [17] (described in the previous section) was used as an input function for the simulation of twitch propagation in the model. Mechanical deformation. The layer of the ventricle was described by the number of muscle units connected in series. As mentioned previously [14], the volume of the ventricle model was estimated for each time step of the analysis. Thus, the length of the outer boundaries L V (circumference) of the ventricle (Fig 2) could be calculated. Using the approach of Shim et. al [20], changes in the volume of the ventricles could be translated into the changes in the halfsarcomere length by dividing the length L V by the number of half sarcomeres (N) in the The updated half-sarcomere length was used to calculate the I SAC current that was added to the TNNP electrophysiological model and, in turn affected the AP characteristics.

Simulation protocol
Several sets of simulations were performed to compare the electrophysiological behavior of the stretched tissue with the SACs recruitment. The first simulation runs were performed to investigate the influence of the SACs on the AP at the single cell level by applying a steady stretch. The cell was excited by applying a stimulus current (52 pAÁpF -1 ) at t = 0 ms with a length of 1 ms. The stretch was induced by changing the resting half-sarcomere length. The effects of the stretch were simulated when the cell was not stretched (λ = 1) and at different degrees of stretch (λ = 1.10, 1.20, 1.25, and 1.30). The second simulation was performed with the 3D model of the contracting ventricle. In these simulations, the AP propagation was compared in the tissue with different half-sarcomere length values (0, 10, 15, 20, and 25% of the reference length) and associated preload values. The stimulus was applied for 6 ms at the predetermined element of the ventricle apex. The local activation time (AT) was determined as the time when the membrane potential upstroke reached -60 mV. The local recovery time (RT) was defined as the time when the repolarized membrane potential reached -60 mV (typically 80% of repolarization). The duration of the AP (APD) of each individual myocardial fiber was calculated as a time interval between the initial activation and recovery of the fiber.
For the next simulations HF pathological condition was incorporated into the TNNP model (a detailed description of the integrated modifications is presented in S1 Table in the Supporting information). The cell was excited under the same conditions as those in the first simulation. No stretch was applied during these simulation runs. At the system level, the corresponding pressure-volume (PV) loops were generated with the fixed preload and the control set of the loading parameters (see [14] for details). The excitation propagation in the 3D model of the contracting ventricle under HF condition was carried out in the same manner as the second simulation.

Results
The contribution of the SACs to the stretch effects on the AP is shown in Fig 3. The effects of the steady stretch were simulated at different degrees of the stretch ratio λ. When the stretch increased, the early phase of repolarization of the AP was accelerated, while late repolarization was prolonged for the stretches of 10%, 20%, and 25% (λ = 1.10, 1.20, and 1.25, respectively). A large stretch (30% of the reference length) depolarized the membrane and prevented repolarization of the AP. The stretches of 10-20% led to an increase in the resting potential by approximately 10 mV. The crossover of the repolarized APs occurred near -20 mV. The peaks of the AP amplitudes obtained during the simulations were independent of the stretch.
The second simulation examined the effects of different preload lengths on the excitation propagation through recruitment of the SACs. Several AP characteristics were compared following the ventricle stimulation. The first AT for λ = 1 (resting length) was observed at the stimulus site at 4.4 ms (Fig 4). No effect was observed on the earliest ATs (the first ATs compared between the fibers) and the latest ATs the last ATs compared between the fibers) when the preload length increased to 20%. A larger increase (25%) produced an acceleration of the latest AT. Fig 5 shows the RTs obtained during the simulation. After the apex stimulation, the first RT occurred at 374.9 ms from the onset of the applied stimulus. Despite the increase of λ, the RTs

Excitation propagation during heart failure
The next set of simulations was performed to examine the influence of different preload lengths on electrophysiological parameters under heart failure condition. Fig 6A showed the time course of the transmembrane voltage at 1 Hz steady state pacing in control and HF conditions. The figure showed the prolongation of simulated HF action potential compared to control (the original TNNP model). The corresponding PV loops were illustrated in Fig 6B. It can be seen that at a constant preload there was a decrease in ventricular pressure and stroke volume in HF PV loop compared to the normal PV loop.
In order to make comparisons, the APDs generated in the tissue with different preload lengths under normal and HF conditions were analyzed. The calculated APDs are listed in  Table 2. The obtained APD 80 values (the action potential duration at 80% repolarization) for varying preload lengths (0%, 10%, 15%, 20%, and 25% of the reference length) were compared. No changes were found in APDs under control condition, while HF resulted in significantly prolonged APDs for preload lengths up to 20% of the tissue reference length. For a 25% increase in the preload length at control condition, the APD 80 was increased and was approximately 50 ms longer than that of the smaller length changes. A 25% increase in the preload length at HF prevented repolarization of the APs.

Discussion
In this study, we investigated the effects of mechano-electric coupling in the heart at the system level using a previously developed 3D multiscale model of the contracting ventricle [14]. To understand the effects, the model was modified by adding the SAC that corresponded to the stretch effects on the AP. The analysis of the MEF effects in the 3D model included simulations with the SACs induced by changing the preload lengths. At the first stage, the effect of a steady stretch on the AP in a single myocyte was assessed. A number of experimental studies have shown that myocardial stretch may shorten or lengthen the APD. A majority of the experiments with rabbit or canine hearts showed a lengthening of the APD during the late repolarization phase [2,12]. Similarly, at a single cell level, this study showed APD shortening at the early repolarization phase and delayed repolarization and crossover of the APs near -20 mV (Fig 3). This was in agreement with the experimental results [2,21].
The second stage of the simulations evaluated the electrophysiological behavior of the 3D model at the system level for different preload lengths with SAC recruitment. The activation and recovery times as well as the minimum and maximum APD for varying preload lengths (0%, 10%, 15%, 20%, and 25%) were estimated and compared. Hansen et al. [22] showed a shortening of the plateau phase of the AP and a delay during phase 3 repolarization owing to the mechanical stretch that occurred from an increase of the ventricular preload. In this study, a 10-20% increase in the preload length had no effect on the AT and RT during excitation propagation. However, the recovery patterns were affected significantly by a 25% increase in the preload length of the tissue reference length. This was in agreement with the results of Hansen et al. [22]. In contrast to the single cell results, during the system level simulations, an increase in the preload length of 20% did not influence the APD 80 values. However, for a 25% increase, the APD80 was enlarged. Similar observations were obtained by Hansen [22], where minor changes in the measured APDs at 90% repolarization while altering the preload values were reported.
The model was used then to assess the influence of different preload lengths on electrophysiological parameters under HF condition. HF was modeled by decreasing the maximum conductance of the potassium currents (I to , I Ks , I K1 ), the sodium current I Na and the calcium current I CaL while increasing the Na + /Ca 2+ exchanger-the mechanisms that contributed to AP prolongation [23]. The simulated action potentials of normal and failing ventricular cells were of similar morphology to measured physiological data (Fig 6A). The corresponding PV loops (Fig 6B) demonstrate an expected decrease in inotropy (less ventricular pressure was generated and stroke volume was decreased) as was observed experimentally [24]. These simulation results showed that the model presented realistic ventricular contractile performance under HF condition.
Results of the subsequent simulations showed that the presence of heart failure had significant effect on excitation propagation. The ATs and RTs, the minimum and maximum APDs for varying preload lengths to 20% were affected. It can be seen that the mechano-electrical feedback effects that we observed appear to be most important in failing hearts with low ejection fraction. When the preload length was above 20% of the tissue reference length, it

Limitations
To ensure the usefulness of the developed model, a number of assumptions were made during design. As there was just one layer of myocardial fibers in the 3D model, transmural heterogeneity was not incorporated in this study. This could potentially affect the activation pattern.
Another limitation refers to the propagation of electrical excitation which was represented by the reaction-diffusion equation for the isotropic medium. Consideration of the anisotropy in anatomically accurate geometries, such as previously developed anatomical 3D model of the contracting heart [25], will allow to describe physiologically accurate electrical propagation in myocardial tissue. This will be done in future.
Because of the absence of homogeneous data for the SACs (e.g. channel conductance, density, and reversal potential), the parameter values were chosen from different data sources. This may have reduced the accuracy of the model.
Mechanical modulation of cellular Ca 2+ was not included in the analysis of the MEF effects. The mechanism by which intracellular Ca 2+ (the affinity of troponin for calcium) is affected by mechanical perturbations, was not within the scope of this study.
Supporting information S1