Interaction of the Mechano-Electrical Feedback With Passive Mechanical Models on a 3D Rat Left Ventricle: A Computational Study

In this paper, we are investigating the interaction between different passive material models and the mechano-electrical feedback (MEF) in cardiac modeling. Various types of passive mechanical laws (nearly incompressible/compressible, polynomial/exponential-type, transversally isotropic/orthotropic material models) are integrated in a fully coupled electromechanical model in order to study their specific influence on the overall MEF behavior. Our computational model is based on a three-dimensional (3D) geometry of a healthy rat left ventricle reconstructed from magnetic resonance imaging (MRI). The electromechanically coupled problem is solved using a fully implicit finite element-based approach. The effects of different passive material models on the MEF are studied with the help of numerical examples. It turns out that there is a significant difference between the behavior of the MEF for compressible and incompressible material models. Numerical results for the incompressible models exhibit that a change in the electrophysiology can be observed such that the transmembrane potential (TP) is unable to reach the resting state in the repolarization phase, and this leads to non-zero relaxation deformations. The most significant and strongest effects of the MEF on the rat cardiac muscle response are observed for the exponential passive material law.


INTRODUCTION
Cardiovascular diseases remain the leading causes of death, e.g., 30% in the US and 45% in Europe (Wilkins et al., 2017), even though the cardiovascular system has been extensively studied. A great hope to reduce the mortality of cardiovascular diseases in the future lies in computational models. These models can be an effective tool to study and understand the cardiovascular system and related pathologies in a new fashion. This includes, e.g., the early recognition of heart failure, better understanding of the cardiac function under normal and pathological conditions (Trayanova, 2011;Gao et al., 2015), patient-specific diagnostics and treatments (Asner et al., 2017), as well as the development of cardiac devices (Baillargeon et al., 2015).
Computational models of the heart include an excitationcontraction coupling representing the physiological course of converting an electrical stimulus into an active muscle contraction. The electrical excitation, which induces depolarization and contraction of cardiac cells, is widely represented by the ionic Hodgkin-Huxley model for neurons (Hodgkin and Huxley, 1952). Based on the Hodgkin-Huxley description, a large number of models have been derived, which can be subdivided into physiological and phenomenological models. On the one hand, physiological models are used from cell to tissue level, for example drug testing of a cell drum compromising three different cardiac stem cells (Frotscher et al., 2016) as well as for drug investigations on a complete human heart (Costabal et al., 2018). On the other hand, phenomenological approaches like the Aliev-Panfilov model and the FitzHugh-Nagumo model are also widely used to model cardiac excitability, mostly at tissue level (Fitzhugh, 1961;Rogers and McCulloch, 1994;Aliev and Panfilov, 1996). Compared to physiological models, the phenomenological approach describes the TP evolution, using a significantly smaller number of internal variables, while still capturing the main characteristics of the cardiac electrophysiology. Thus, it is less complex, easier to implement, reduces computational costs (Cherubini et al., 2008) and is often used to investigate arrhythmia.
Apart from the excitation-induced contraction of cardiac cells, the depolarization and subsequently the contraction of cardiomyocytes can be also initiated through the stretchinduced opening of ion channels, commonly referred to as the mechano-electrical feedback (MEF) (Kohl et al., 1999;Keldermann et al., 2009). More specifically, the MEF is capable of modifying the electrophysiology (Kamkin et al., 2005) and might promote stretch-activated arrhythmias, which have been commonly identified as a result of electrical disorders in the cardiac muscle (Franz, 1996). Moreover, the MEF can shorten the action potential (AP) duration or shorten the time course of repolarization (Franz, 1996) and break up spiral waves . Furthermore, it has been shown that high strain and stretch rates can cause premature ventricular excitation (Franz et al., 1992;Hu and Sachs, 1997). In addition to the extensive investigations on humans, similar findings have been observed for arterial cardiomyocytes in rats (Kamkin et al., 2000). In particular, stretch can result in atrial fibrillation after ventricular infarction (Kamkin et al., 2000). Also in other mammals like rabbits, the MEF can give rise to spontaneous arrhythmia in an acute local ischemia (Jie et al., 2010). This phenomenon plays a key role in interpreting the interplay between the electrophysiology and mechanics of cardiomyocytes . In a comprehensive computational modeling approach for cardiac electromechanics, it is essential to account not only for the excitation-triggered contraction of cardiac myocytes but also for the stretch-activated excitation. To investigate the effect of the MEF, strongly coupled electromechanical models are used (Costabal et al., 2017). However, depending on each specific application with different imposed goals in heart research, weakly coupled problems are also considered and sequentially solved (Frotscher et al., 2016;Duong et al., 2018). On the other hand, decoupled or staggered approaches are suitable for one-way coupling formulation (Usyk et al., 2002;Nash and Panfilov, 2004;Baillargeon et al., 2014;Frotscher et al., 2016). While the effects of the MEF have been widely studied in experiments, the effects in connection with various computational models, especially different passive material models, remain relatively unclear. In this study, the MEF is investigated in combination with various types of passive mechanical laws (nearly incompressible/compressible, polynomial/exponential-type, transversally isotropic/orthotropic material models) in order to study their interaction with the MEF.
Specifically, we want to focus on the interaction of the passive mechanical model with the MEF and how different types of material laws (compressible, incompressible, polynomial, exponential, transversely isotropic, orthotropic) are influencing the overall MEF characteristics. We employ a transversely isotropic and nearly incompressible model, a transversely isotropic compressible model (Göktepe and Kuhl, 2010) and the orthotropic and nearly incompressible Holzapfel-Odgen model (Holzapfel and Ogden, 2009). Our study is based on a left ventricle of a rat heart, whose geometry is reconstructed from high resolution MRI at the University of Erlangen-Nuremberg (Duong et al., 2018). The phenomenological model of Aliev-Panfilov is used to represent the electrical excitation (Aliev and Panfilov, 1996). We investigate the variation in TP evolution and mechanical deformation due to the interaction between the passive mechanical models and the MEF.

METHODS
We briefly introduce the basic electromechanical model and numerical methods used to formulate and solve the boundary value problem of the contracting ventricle. More details on the kinematics and constitutive equations can be found in the Appendix A.

Governing Equations for Cardiac Electromechanics
The electromechanics of a left ventricle (LV) (see Figure 1) can be described by two primary field variables, placement ϕ(X, t) and AP Φ(X, t). Thus, two field equations, which govern the state of the material point X at time t, t ∈ [t 0 , t f ] (t 0 and t f are the initial and final time, respectively), can be formulated. The mechanical field equation is the balance of linear momentum where F is the deformation gradient, S is the second Piola-Kirchoff stress tensor and F Φ is the external mechanical body force. The other differential equation describes the spatiotemporal evolution of the AP Φ and can be written asΦ FIGURE 1 | Boundary surface Ŵ 0 decomposed into Ŵ ϕ and Ŵ T for the mechanical model and Ŵ Φ and Ŵ Q for the electrical model.
The non-linear current term F Φ takes into account the electrical excitation,Φ denotes the material time derivative of the AP field, Div[Q] represents the diffusion of the AP, whereby Q describes the electrical flux vector. Together with the boundary conditions, the mechanical and electrical balance equations constitute an initial boundary value problem in the strong form for unknown placement and AP, {ϕ(X, t), Φ(X, t)}. The boundary conditions are Dirichlet and Neumann conditions as ϕ(X, t) =φ on Ŵ ϕ , surface traction vector T(X, t) =T on Ŵ T for all t ∈ [t 0 , t f ] and ϕ(X, t 0 ) = 0 in 0 . For the electrical model, the boundary conditions are Φ(X, t 0 ) =Φ on Ŵ Φ and Q(X, t) · N =Q on Ŵ Q with the unit normal vector N pointing outwards of surface Ŵ Q (see Figure 1). By the Cauchy stress theorem, we get T = F · S · N.

Mechanical Constitutive Models
By the active stress approach (see e.g., Smith et al., 2004;Göktepe and Kuhl, 2010;Pathmanathan et al., 2010), the stress response of cardiac muscles can be decomposed into a passive and an active part. While the passive part is only governed by the mechanical deformation, the active part takes into account the excitation-induced contraction. Hence, the total second Piola-Kirchhoff stress tensor S is written as where C denotes the right Cauchy-Green deformation tensor, f 0 the fiber direction and s 0 the sheet direction in the material configuration. Now, three different models for the passive mechanical response are presented which are utilized to investigate their interaction with the MEF. From the strain energy function , which is used to describe the mechanics of soft tissues, the associated mechanical constitutive equation reads All further derived quantities can be decomposed into their passive ( ) pas and active part ( ) act .

Transversely Isotropic Compressible Model (TIC)
The first passive material model we are utilizing is a transversally isotropic compressible model, which is represented in the form of a polynomial function (see Göktepe and Kuhl, 2010). The passive material stress response is described by a modified neo-Hookean constitutive model in which the basic neo-Hookean model is extended by a transversely isotropic part taking into account the dependency of the material properties in the fiber direction f 0 at each point X in the material configuration. The strain energy function can be written as in terms of the principal strain invariants of the symmetric right Cauchy-Green tensor C as I 1 = tr(C), det(C) = J 2 , and I 4f (C) = f 0 · (Cf 0 ), where the Jacobian J is defined as J = det(F). In Equation (5), the two constants and µ in the isotropic part iso are the Lamé parameters, and the passive stiffness of myofibers is denoted by η in the transversely isotropic (so-called anisotropic) part ani .

Transversely Isotropic and Nearly Incompressible Model (TII)
The second polynomial-type model is basically derived from the above compressible modified neo-Hookean TIC model, but it is rewritten in a form such that it is transversely isotropic and nearly incompressible reading in terms of the principal invariants of the isochoric Cauchy-Green tensorC =F TF . The principal isochoric invariants ofC are defined asĪ 1 (C) = tr(C) andĪ 4f (C) = f 0 · (Cf 0 ). Further, vol = κ(J − 1) 2 is the volumetric energy, where κ tunes the enforcement of incompressibility (we use κ=10 4 kPa). Similar as in TIC, the anisotropic part of the free energy function ani only occurs if the fibers are under tensionλ > 1 withλ = Ī 4f .

Orthotropic and Nearly Incompressible Model (HO)
The third considered model is an exponential-type strain energy function, the Holzapfel-Ogden model (see Holzapfel and Ogden, 2009). In the field of cardiac modeling, the HO model is one of the most suitable choices for describing the passive mechanical response of the myocardial tissues since it can capture the hyperelastic and orthotropic characteristics which have been found in experiments on porcine hearts (see Dokos et al., 2002). In contrast to the polynomial-type models TIC and TII, the HO is based on a strain energy function which is represented by exponentials. The orthotropy at each point X is characterized by a right-handed set of normalized basis vectors which are determined by the fiber direction f 0 , the sheet direction s 0 and the orthogonal vector of the sheet plane n 0 = f 0 ×s 0 . By applying the incompressibility condition to the finite element setting, the HO model is, therefore, split into its isochoric term (in terms of the principal isochoric invariants) and a volumetric part vol (as in TII). The strain energy function reads as where i ∈ {f , s} and the variables a, b, a f , b f , a s , b s , a fs , b fs are material constants. While all a, a f , a s , a fs parameters have the dimension of stress, all b, b f , b s , b fs are dimensionless. The active stress response is described in the Appendix A.2.

Electrophysiological Constitutive Models
The nonlinear current source term F Φ controlling the AP in Equation (2) is split into two parts as where F Φ e expresses the purely electrical part and F Φ m accounts for the MEF, i.e., a mechanically-induced excitation. The excitationinduced purely electrical part F Φ e characterizes the effective current, which is generated from the inward and outward flow of ions across the cardiac cell membrane. Meanwhile the stretchinduced mechano-electrical part F Φ m accounts for the opening of ion channels due to the mechanical deformation. By introducing the non-dimensional and normalized action potential Φ and the non-dimensional timet, we derive the conversion forms as (Aliev and Panfilov, 1996) where k φ and δ φ relate the dimensionless action potential φ to the physical action potential Φ and k t the dimensionless timet to the physical time t. Consequently, the purely electrical term F Φ e is modeled as where the coefficient α controls the oscillation threshold, c is a scaling parameter, I is an external stimulus and r is the recovery variable which controls the repolarization of the cardiac cell.

Potential Flux
The electrical constitutive equations are also formulated as functions of the deformation gradient and the AP, the material electrical flux reads with the conductivity tensor where D iso = Jd iso and D ani = Jd ani . The parameter d iso accounts for the isotropic and d ani for the additional faster conduction along the fiber direction (Costabal et al., 2017).

Mechano-Electrical Feedback
As mentioned before, the electro-mechanical coupling plays an essential role in the cardiac function. A mechanical cell deformation induces electrical current generation. This behavior is observed due to stretch-induced opening of ion channels which induces AP generation. It can be described by the constitutive equation for the electrical source term F Φ m as where G s denotes the maximum conductance, φ is given in Equation (15), φ s is the dimensionless reversal potential at which there is no net ion flux through the stretch-activated channels, ϑ is a switch function turning the feedback on for λ > 1 and off else. λ = I 4f (see Panfilov et al., 2005;Keldermann et al., 2007).

Excitation Stimulated by Deformation in a Plate
In this section, we refer to the benchmark problem performed on a plate using the compressible TIC model to study the effect of the MEF in which excitation wavefronts are observed caused by the deformation-induced excitation of cardiac tissue (Göktepe and Kuhl, 2010;Dal et al., 2013;Cansiz et al., 2015). The orthotropic conductivity is employed in the plate with d iso = 1.0 mm 2 ms − 1 and d ani = 0.1 mm 2 ms − 1, where according to section 2.3.1, the latter accounts for the additional conduction in fiber direction increased by 10% with respect to the other directions (Göktepe and Kuhl, 2009).
To illustrate the deformation-induced excitation in the plate, a mechanical load p(t) is applied during t ∈ [0, 10] ms as shown in Figure 2A. The plate with the dimensions of 100 × 100 × 12 mm is meshed by 21 × 21 × 2 eight-node brick elements. The fiber orientation f 0 and the sheet plane direction s 0 are defined in x-and y-direction, respectively. The compressible model TIC is employed with parameters given in Table A1. The maximum conductance G s = 15 is utilized to account for the MEF effect. Edge nodes of the plane in the middle at z = 6 mm are fixed in z-direction. Furthermore, the node at (0, 0, 0) is fixed in xand y-directions and the node at (100, 0, 0) in y-direction. At the beginning of the simulation, the cellular TP in the plate is set to the resting value Φ = −80 mV. To trigger the MEF, the nodes of a parallelepiped located in the center region of the plate with a dimension of 20 × 20 × 12 mm are subject to an impulsive cyclic loading p(t) from t = 0 to t = 10 ms. It linearly reaches its maximum of 0.3 N at t = 5 ms and returns to 0 N at t = 10 ms.
This loading results in a stretched region in the middle of the plate and hence gives rise to a local depolarization, which then initiates an excitation wave traveling first elliptically (due to orthotropic conductivity, see t = 12, 21, and 33 ms) and then unidirectionally along the y-direction. At the same time, the plate shortens in the fiber direction x and elongates along the other directions orthogonal to the fibers (t = 61 and 74 ms). Since the model is compressible, the whole plate contracts and reduces its volume. The action potential impulse travels to the plate sides and the plate starts to repolarize with a lower potential wave starting in the middle as an elliptical shape (t = 128 and 163 ms). While repolarizing, the plate recovers the initial volume due to the relaxation phase of the cardiac muscles. This illustrates the change of TP due to the MEF, which will be discussed in more detail in the next section for different types of mechanical models. It is worth noting that for the TIC model, the applied mechanical load needs to be relatively large in order to cause a visible MEF effect. The same test is performed on a plate of the rat heart dimension (10 × 10 × 1.2 mm) with parameters used in the later simulations (see Table A1). We observed similar behavior; however due to the small dimension of the plate, complete depolarization is reached at t = 47 ms ( Figure 2B).

Rat Heart Measurements
To generate a 3D left ventricle model of a rat using MRI data, we conducted several experiments on living and healthy rats. More importantly, we also confirm that the ethics committee (Tierschutzgesetz Regierung Unterfranken) approved our experimental protocol and procedures. Further, only the best measurement result was used to reconstruct the 3D rat left ventricle model in the diastolic phase.

RESULTS
In this section, we investigate the effect of different passive material laws in a strongly coupled electromechanical model of a 3D rat left ventricle.

Parameters of Cardiac Muscles of Rat Heart
The material parameters for our simulation are obtained by curve fitting to experimental data of a porcine heart by Dokos et al. (2002) (see Figure 3; using a scaling factor of 0.5 for rats). The resulting material parameters are displayed in Table A1 (Appendix A.5). The parameters for the electrical and the active model originate partially from our parameter study and partially from work for healthy human hearts (Aliev and Panfilov, 1996;Göktepe and Kuhl, 2010;Baillargeon et al., 2015). We also use k T = 0.49 kPa mV −1 in Equation (24), resulting in T act max = 49 kPa, which is sufficiently close to the maximum tension value of 45 kPa as obtained in rat experiments (see Niederer et al., 2009). The reversal potential φ s in Equation (15) is set to 0.6, corresponding to −20 mV, which is in agreement with the physiological value by Kohl et al. (1999).

Rat Left Ventricle
In this section, the interaction of the three different passive material models with the MEF is investigated with regard to the influence of the maximum conductance with G s = 10 and G s = 0. For the electromechanical simulation, the parameters given in Table A1 (Appendix A.5) are used. The base of the ventricle is mechanically fixed (see Figure 1). However, different boundary conditions (Baillargeon et al., 2015) can be applied such that the base of the ventricle can slightly translate or twist when it contracts and interacts with surrounding tissues, which greatly support and stabilize the whole heart. The mesh of the 3D solid left ventricle model is generated using 4-node tetrahedral elements based on a high resolution MRI from the Universitätsklinikum Erlangen-Nürnberg. The basic workflow from the MRI images to the finite element mesh can be seen in Figure 4.
The fiber and sheet orientations, which are crucially attributing to the mechanics and electrical conduction system, are assigned for the LV by interpolating the local fiber and sheet directions from the endocardial and epicardial surface such that the model can account for the transmural fiber and sheet directions (Vetter and McCulloch, 1998;Wong and Kuhl, 2014). Although the fiber angles on the endocardium and the epicardium vary largely between different rats (Chen et al., 2003; Mekkaoui et al., 2012), the fiber angles are chosen as +80 • on the endocardium and −70 • on the epicardium (see Figure 5). For simplicity, we assume a chamber pressure of zero, no electrical flux on the boundary as well as zero surface traction. The computational mesh of the 3D LV consists of 49,769 tetrahedral elements corresponding to the global mesh size of 0.35 mm. Initially, all 10,237 nodes are set to the resting potential Φ r = −80 mV. Seven nodes at the base are used to trigger the depolarization in the surrounding cardiomyocytes. They are constrained at Φ = −20 mV for a duration of 40 ms. Subsequently, we solve six fully coupled electromechanical problems for the three different passive material models with and without the MEF (TIC G s =10 , TIC G s =0 , TII G s =10 , TII G s =0 , HO G s =10 , HO G s =0 ). We performed simulation on refined meshes up to a mesh size of 0.15 mm. The plot in Figure 6 illustrates that neither the maximal, nor the residual TP show significant changes in regime of smaller meshes. In particular, the resting potential of −80 mV is not reached for all tested mesh sizes. Figure 6 exemplarily shows the numerical result for the HO model with the MEF (HO G s =10 )-the location of the stimulus can be seen at the base at t = 2 ms. From there, the electrical impulse propagates through the ventricle, see snapshots at t = 60 and 70 ms. The entire myocardium contracts at t = 90 ms when all ventricular cells are depolarized (in red). After the depolarization, the repolarization follows (t = 120, 150, 176, and 300 ms), but the ventricle is unable to return to the initial shape (light blue shadow) even for t ≥ 300 ms. This behavior corresponds to a second peak in the AP curve for the HO model (peak 2 above −80 mV) after peak 1 (see Figure 7, right).

Hales
Moreover, the AP curve for the TII model with G s = 10 in the repolarization phase reaches a minimum value of about −78.5 mV for t ≥ 300 ms, while the AP curve for the TIC model is able to retrieve the resting potential of −80 mV (see Figure 7, right). The same characteristic behavior can  be observed by plotting the apex displacement u (Figure 7, right). On the contrary, Figure 7 (left) exhibits no peak 2 in the AP and displacement curves for all three models without MEF (TIC G s =0 , TII G s =0 , HO G s =0 ). The LV consequently relaxes to its initial shape. The stretch-induced excitation changes the AP duration, leads to a delayed repolarization and thus the resting potential is reached more slowly. In Figure 7 (left), the AP evolution profiles for the three models are very similar; however, the displacement predicted by the compressible TIC model is larger than the displacements predicted by the TII and HO models. Compared to the AP curves, the displacement of the apex starts before the depolarization. This is caused by the fact that the cardiac cells between the initiated nodes at the base and the considered node at the apex already depolarize and start to contract and hence cause the pre-displacement of the apex.
In summary, it can be stated that both incompressible models yield a different electrical and mechanical response (residual deformation and resting potential not reached at the end of the cardiac cycle) of the LV compared to the compressible TIC model. This raises the question, what exactly causes the difference in the results among the three passive mechanical models. On the one hand, the HO represents an orthotropic material law which is represented by exponentials, while the other two models are transversely isotropic and formulated in terms of polynomials. On the other hand, the HO and TII are both incompressible models, while TII is compressible. To study this more comprehensively, we extended the numerical investigations by performing a parameter study on the maximum conductance G s , which is commonly known as the maximum stretch-activated ion channel conductance or the sensitivity of the electrical current to deformation. Note that diseased hearts can cause abnormal changes in the maximum conductance (Zhang et al., 2014). Thus, this value can be considered as a key variable to study the MEF effects.
We performed 18 simulations, which differ in the passive material model (TIC, TII, HO) and the maximum conductance (G s = 0, 10, 20, 30, 40, and 50). Instead of performing a benchmark problem study using a simple cube, bar or beam, the effect is directly studied in a 3D LV model. The result of the AP evolution and the displacement are shown in Figure 8. The TIC model shows an almost independent behavior with respect to the strength of the MEF [see AP and displacement evolution Figure 8 (top, left, and right)]. The polynomial, transversely isotropic and incompressible material law TII (middle, left, and right) and the exponential, orthotropic and incompressible material law HO (bottom, left, and right) exhibit a high sensitivity to the strength of the MEF. The observed characteristics in Figure 7 (early depolarization, delay in the late repolarization phase, decrease of the maximum AP) are intensified with a higher maximum conductance.  Table A1. Excitation by setting seven nodes at the base to Φ = −20 mV for 40 ms. At the end of cardiac cycle (around t = 300 ms), the repolarised (resting) LV (blue) is unable to return to its initial shape (light blue region) since the MEF affects the AP and deformation of the LV. Plot at the right bottom shows the relation of maximal and residual potential versus the global mesh size (HO model with MEF and G s = 50).
FIGURE 7 | Evolution of AP Φ and displacement u of apical node A in a rat left ventricle for three models, without MEF (Left) and with MEF using G s = 10 (Right).

DISCUSSION
As observed in the last section, the influence of the MEF on the overall electrophysiological and mechanical behavior significantly depends on the passive material model. In the following section, we want to discuss the differences and possible sources or causes. The main differences of the three different passive material models can be summarized as follows: (i) compressible vs. incompressible, (ii) polynomial vs. exponential, (iii) orthotropic vs. transversely isotropic. Their respective results will be evaluated and compared here. Referring to the following explication, we introduce the notation: max−a and rest−a for values of the material a at peak 1 (e.g., the maximal value of the AP or displacement) and at peak 2 close to the resting potential −80 mV and zero deformation (see Figure 7, right), respectively. We define the difference of a value at G s = x to that at G s = 0 as △ First we note that the MEF is capable of speeding up the electrical impulse propagating in the LV (Figure 9). For the compressible TIC model, it is clearly observed that the MEF has negligible effects on the electrophysiology and mechanics.
The largest change in time for the AP peak is △ t 50 max−TIC = 5 ms while there is no peak 2 and thus △ t x rest−TIC = 0 for all values x ∈ {0, 10, 20, 30, 40, 50} (see also peak values in Figure 9). In contrast to that, we observe a significant effect of the MEF (or maximum conductance) for the incompressible models TII and HO. Figure 8 indicates that for higher value of G s , larger changes in AP and displacement curves can be observed. The largest changes in peak time for AP are △ t 50 max−TII = 18 ms and △ t 50 max−HO = 23 ms for the TII and HO models, respectively. The conduction velocity is apparently increased with G s which is in good agreement with the findings in a study by Costabal et al. (2017). Similar observations have also been made in the numerical results by Amar et al. (2018), in which the AP was influenced by the MEF for different stretch levels for examples of a ventricle and a single cardiomyocyte.
Second, depending on the material model, the MEF influences AP Φ and displacement u in the late repolarization phase. For the compressible model TIC, no change in peak 2 for the AP and displacement curves (△ Φ x rest−TIC = 0 mV and FIGURE 9 | Time points for corresponding two AP peaks for the three models. △ u x rest−TIC = 0 mm) for all values x ∈ {0, 10, 20, 30, 40, 50} can be observed. At the same time, there are significant changes in peak 2, which can be observed in the curves for the TII and HO model. For example, the largest changes are introduced by G s = 50 such as △ Φ 50 rest−TII = 16 mV and △ Φ 50 rest−HO = 22 mV, whereas △ u 50 rest−TII = 0.45 mm and △ u 50 rest−HO = 0.79 mm (see Figure 10, right). The maximum value of the AP Φ and displacement u for peak 1 for (G s = x) for all values x ∈ {0, 10, 20, 30, 40, 50} are depicted in Figure 10 (left). While varying the maximum conductance G s , the compressible model leads to an almost unchanged electrophysiology and mechanics, the nearly incompressible models show a stronger effect of the MEF on the electrophysiology (early depolarization, reduced maximum AP, increased conduction velocity, delayed repolarization) and consequently result in a significantly higher relaxation displacement of the LV.

Compressible vs. Incompressible
One important characteristic which influences the impact of the MEF on the overall model behavior is the level of compressibility. The compressible TIC model undergoes a significant volume change (compressible medium) during the active contraction phase, which in turn leads to a reduced overall strain/stretch. This is visualized in Figure 11 showing the isotonic contraction of a cube for the three different material models.
As expected, when considering the volume ratio curves for the three models in isotonic contraction, the volume ratio during the time course of the cardiac cycle for the compressible TIC model decreases ineluctably, while the volume for the HO and TII models remains constant. By plotting the fiber stretch λ of the same apical node of the LV for the HO, TII and TIC model, we observe the fiber stretch λ < 1 (compression) for the TIC model and the fiber stretch λ > 1 for the HO and TII models (see Figure 12).
This leads to the fact that in this particular region, the MEF is active for the TII and HO models and inactive for the TIC model (see Equation 18, λ has to be >1 for the switch function ϑ to be non-zero). Furthermore, the fiber angles transmurally vary and thus the local change in fiber direction leads to an FIGURE 10 | Φ max and u max of peak 1 (Left), Φ rest and u rest of peak 2 (Right) for the three models (Φ rest−TIC curves invisible since u rest−TIC lies on top of it).
Frontiers in Physiology | www.frontiersin.org  increased fiber stretch through the wall (the expansion of the material in sheet and normal direction close to the boundaries leads to a fiber stretch in the middle layers; see Figure 13, left, center).
This behavior is characteristic for the incompressible TII and HO models. The TII model compensates this phenomenon due to the compressibility (volume change), and thus the expansion of the material in the sheet and normal direction is smaller, which finally induces no or insignificant stretches to the fibers in the middle layer (see Figure 13, right). For the TII and HO model, the high fiber stretch in the middle layers during the repolarization phase effectively causes a strong MEF current. In a real healthy rat left ventricle, the MEF effect might be caused by a strong mechanical stimulus when the AP reaches the threshold value. However, it is observed mainly in diseased cases (Kamkin et al., 2000). The MEF can alter the AP or shorten refractory periods in canine left ventricle and atrium, as reported in Franz (1996). Stretchactivated polarization and fibrillation can happen at a low stretch in a rat atrium after the infarction of the left ventricle (Kamkin et al., 2000). Moreover, it has been observed in our simulation of the LV as well as for cubes, plates and other 2D examples in Costabal et al. (2017), that the MEF is capable of altering the AP and increasing the conduction velocity. However, the phenomenon of the residual deformation (nonrelaxed equilibrium between the AP generation and the fiber stretches in the late repolarization phase) is rather unrealistic in the healthy left ventricle simulation.

Polynomial vs. Exponential
In addition to the discussed differences due to the level of compressibility, the exponential and polynomial form of the passive material law plays a significant role concerning the overall MEF behavior. To clarify this, we performed the numerical tests for the HO and TII models with the same bulk modulus of κ = 10 4 kPa, which means that these models have the same level of compressibility. The compressibility value is reasonable in our simulations as it is only slightly larger compared to the in vivo measurements by Hassaballah et al. (2013). To compare the influence of the polynomial material law TII and the exponential material law HO, we evaluate the differences in Figures 8, 10, 12, 13.
In the depolarization phase, the HO and TII models show an almost similar behavior concerning the AP and displacement curves even if the observed phenomenon is slightly more prominent for the HO model (see Figures 10, 8, left). Nevertheless, in the late repolarization phase, the AP and displacement curves significantly differ between the HO and TII models (see Figure 10, right). It is worth noting that both models undergo a comparable maximum stretch level (see Figure 13, left, center). Thus the question arises regarding where the significant differences in the late repolarization come from.
Since the same bulk modulus is used, both models generate the same amount of active stress and result in relatively comparable displacements (see Figure 13, left, center). Therefore, the difference can only be explained by the fact that the HO model induces larger strain/stretch in the late repolarization phase and in turn produces a higher MEF current compared to the TII model. Figure 12 exposes that the HO model has a steeper increase in the fiber stretch in the depolarization phase (which explains the faster depolarization compared to the TII model) and a slower change of the fiber stretch in the repolarization phase (which explains the higher stretch in the late repolarization phase). In other words, starting from the same deformation level during the contraction for the HO and TII models, the cardiac tissue relaxes in the repolarization phase, in which the exponential-type stress-strain relation of the HO model responds to the decreasing active fiber tension with a smaller change in fiber stretch compared to the TII model. Finally, this leads to a high MEF current in certain regions, which almost reached their resting potential.

Transversely Isotropic vs. Orthotropic
As we discussed in the previous paragraph, the main difference in the AP and displacement for the HO and TII models originate in the representation of the material law. Nevertheless, the incompressibility condition with the associated stretches along the sheet and sheet normal direction at the boundaries [see Figure 13 (left, center), fiber stretch λ < 1, and thus the compression in fiber direction has to be compensated by stretch in the sheet and sheet normal direction] leads to a stretch in fiber direction in the middle layers. As the fiber stretch λ depends on the material properties in sheet and sheet normal direction, there also exists a difference between a transversely isotropic and an orthotropic material law concerning the MEF behavior. To further investigate the difference, one would have to vary the passive material properties in sheet or sheet normal direction for the same material model (note that for the special case of transversal isotropy, the properties in sheet and sheet normal direction are the same) in order to eliminate the influence of the type of material model (polynomial and exponential). However, this goes beyond the scope of this study. In our case, the difference between the orthotropic material model HO and the transversely isotropic material model TII seems to be rather small as the overall fiber stretch levels are similar.

Limitations of the Work
In the following, we want to briefly discuss the limitations and assumptions for this study. Firstly, we ignore the possible influence of the right ventricle and other components which also contribute to the overall model behavior. Secondly, the electromechanical model uses several parameters, which are directly adapted or come from other studies on human and porcine hearts. Furthermore, the electrophysiological model is an obvious simplification of the actual electrophysiology of the heart. There exist different kinds of cardiac cells in the heart with different contractility and conductivity. In particular, the fast conducting Purkinje fibers need to be included in order to improve the choreographed depolarization and repolarization of the heart. However, the simple electrophysiological model is to be preferred for our application whereby a complex rat heart support system will be developed in order to compensate the function of the pathological heart. Further, more detailed experimental data, especially for rats, about the electrical and active model including the determination of the active muscle tension, transmembrane potential and maximum conductance, G s are required. Measurements to identify the passive material parameters for rats would be desired to further improve the accuracy of the model. Moreover, the fixed base in the simulation model simplifies the actual support of the heart in the body (e.g., connective tissue, atria) and an advanced boundary condition should be applied for future studies. Further investigations are necessary to evaluate the influence of the specific morphology of the heart (different fiber distributions). Additionally, to more precisely describe the exact influence and differences of transversely isotropic and orthotropic material laws concerning the MEF behavior, a parameter study about the level of orthotropy is necessary.

CONCLUSION
In this paper, we focus on the interaction of passive mechanical models with the MEF and how different types of material laws (compressible, incompressible, polynomial, exponential, transversely isotropic, and orthotropic) influence the overall MEF characteristics. We employ a transversely isotropic and nearly incompressible model (TII), a transversely isotropic compressible model (TIC) (Göktepe and Kuhl, 2010) and the orthotropic and nearly incompressible Holzapfel-Odgen model (HO) (Holzapfel and Ogden, 2009). We investigate the variation in AP evolution and mechanical deformation due to the interaction between the passive mechanical models and the MEF.
The interaction between the passive models and the MEF is discussed through a computational study of a rat LV, whereby the following findings are obtained: (i) compressibility: the transversely isotropic material law (TIC) predicts a significantly smaller fiber stretch (compression almost everywhere in the LV) and thus leads to a nearly unrecognizable change in the overall MEF behavior (change in electrophysiology and mechanical contraction); additionally, for the incompressible models, we observe a residual deformation caused by a nonrelaxed equilibrium; (ii) polynomial vs. exponential material laws: due to the exponential strain energy function, the HO model shows a faster temporal change of the fiber stretch in the depolarization phase (leading to a faster depolarization compared to the TII model; higher MEF current) and at the same time slower temporal change of the fiber stretch in the late repolarization phase; (iii) transversely isotropic vs. orthotropic: the incompressibility condition with the associated stretches along the sheet and sheet normal direction on the boundaries [see Figure 13 (left, center), fiber stretch λ < 1 and thus the compression in fiber direction has to be compensated by stretch in the sheet and sheet normal direction] leads to a stretch in fiber direction in the middle layers. As the fiber stretch λ depends on the material properties in the sheet and sheet normal direction, there also exists a difference in the electromechanical behavior and the MEF between a transversely isotropic and an orthotropic material law concerning the MEF behavior.
Obviously, the type of passive material model plays a key role in defining the MEF behavior in a fully coupled electromechanical model. It has to be further investigated which of the considered models reflect the cardiac tissue best concerning the overall MEF behavior.

ETHICS STATEMENT
The experimental protocol and procedures were approved by the ethics committee (Tierschutzgesetz Regierung Unterfranken).

AUTHOR CONTRIBUTIONS
MD mainly contributed to this paper concerning the manuscript, modeling, and simulation. DH contributed through the geometric and fiber modeling, simulation as well as the revision of the paper. MA, SD, and SL contributed through the conception and design of the study. MA provided the MRI data and the segmentation.

FUNDING
The authors thank the Klaus Tschira Stiftung grant 00.289.2016 for funding support.
A. APPENDIX

A.1. Kinematics and Finite Element Approximation
The deformation gradient is defined as: where the determinant of the deformation gradient J is also known as the volume ratio and has to be positive as the body is impenetrable. Further, the right Cauchy-Green tensor C can be introduced as C = F T F. To enforce the material incompressibility conditions (J = 1) in a framework of a finite element setting, the deformation gradient can be split into two parts F = (J 1/3 I)F, where I is the identity tensor. Specifically, J 1/3 I describes purely volumetric deformation whereasF denotes the purely isochoric deformation (J = det(F) = 1). Employing the Galerkin procedure, the residuals (2) and (1) are multiplied with the scalar-and vector-valued test functions δΦ and δϕ which satisfy δΦ = 0 on Ŵ Φ and δϕ = 0 on Ŵ ϕ , respectively. These resulting expressions are further reformulated by using integration by parts and the Gauss theorem over the body 0 such that the weak forms are obtained as: Further, the Dirichlet boundary conditions prescribe the state of the respective surface points to beφ andΦ. The Neumann boundary conditions prescribe the surface tractionT and the surface flux termQ = Q · N. F Φ . All quantities of the boundary conditions are supposed to be given. In addition, in our problems T,Q and F Φ are independent of deformation and AP. In what follows, the notations are introduced: the material time derivative as {˙ } = d{ }/dt and the material divergence and gradient as Div{ } = ∂{ }/∂X : I and ∇{ } = ∂{ }/∂X. The equation of motion expresses the quasi-static force equilibrium and has to hold for all points X in 0 . In the following, we will review how to derive the main terms in the (A2).

A.2.3. Orthotropic and Nearly Incompressible Model (HO)
Based on (11), the passive second Piola-Kirchhoff stress is given as S pas = S iso + S ani , where

A.2.4. Transversely Isotropic Active Stress Response
Taking into account f 0 only, the active second Piola-Kirchhoff stress is defined as follows: The transversely isotropic characteristic is realized such that the magnitude of the active fiber tension T act (Φ), which is driven by the AP, only has an effect along the material fiber orientation f 0 .

A.2.5. Orthotropic Active Stress Response
When both f 0 and s 0 are considered for active contraction, the stress is written as: S act (f 0 , s 0 , Φ) = T act (Φ) ν ff f 0 ⊗ f 0 + ν ss s 0 ⊗ s 0 , where ν ff and ν ss are weighting factors. Furthermore, T act (Φ) is modeled via the evolution equation asṪ act = T(Φ, T act ) in Nash and Panfilov (2004). According to the local ordinary differential equation of the active muscle traction, T act can be treated as an internal variable and locally updated on the integration point level. The evolution equation for the muscle traction is herein used for both active tension models and reads: with its sensitivity with respect to the action potential.