A simple phenomenological approach for myocardial contraction: formulation, parameter sensitivity study and applications in organ level simulations

Contraction in myocardial tissue is the result of a complex process through which chemical energy on the cellular level is converted into the mechanical energy needed to circulate blood throughout the body. Due to its vital role for the organism, myocardial contractility is one of the most intensively investigated subjects in medical research. In this contribution, we suggest a novel phenomenological approach for myocardial contraction that is capable of producing realistic intracellular calcium concentration (ICC) and myocyte shortening graphs, can be easily calibrated to capture different ICC and contraction characteristics and, at the same time, is straightforward to implement and ensures efficient computer simulations. This study is inspired by the fact that existing models for myocardial contractility either contain a number of complex equations and material parameters, which reduce their feasibility, or are very simple and cannot accurately mimic reality, which eventually influences the realm of computer simulations. The proposed model in this manuscript considers first the evolution of the ICC through a logarithmic-type ordinary differential equation (ODE) having the normalized transmembrane potential as the argument. The ICC is further put into an exponential-type ODE which determines the shortening of the myocyte (active stretch). The developed approach can be incorporated with phenomenological or biophysically based models of cardiac electrophysiology. Through examples on the material level, we demonstrate that the shape of the ICC and myocardial shortening curves can be easily modified and accurately fitted to experimental data obtained from rat and mouse hearts. Moreover, the performance of the model in organ level simulations is illustrated through several multi-field initial-boundary value problems in which we show variations in volume-time relations, heterogeneous characteristics of myocardial contraction and application of a drug in a virtual left ventricle model.

which in turn may lead to heart failure. Revealing the complex mechanisms behind the myocardial contraction in health and disease is one of the most studied aspects of the heart. In this regard, a big number of experimental and clinical studies have been conducted for both humans and animals. Apart from that, the last decades have witnessed extensive utilization of numerical modelling approaches inquiring into myocardial contraction on the cellular and organ level. Undoubtedly, numerical modelling can be considered a powerful tool for investigating myocardial contraction and its implications on heart function in health and disease. In this paper, we contribute to the literature by elaborating novel phenomenological equations for intracellular calcium concentration (ICC) and myocardial shortening which are simple and straightforward to implement and, at the same time, enable us to mimic reality in an accurate way and perform efficient computer simulations of the virtual heart models. It should be noted that the equations for the ICC and myocardial contraction were first used by Cansız et al. [5] and further adopted by Cansız et al. [6,7] in order to describe myocardial contraction. However, a detailed investigation of the equations at the material and organ level will be carried out for the first time in this manuscript.

Myocardial contraction
Myocardial contraction stems from the tensile force generated by sarcomeres after a series of events. Essentially, the coupling between the electrophysiology and mechanical activity in the myocytes is linked through the ICC. When a myocyte becomes depolarized, Ca 2+ in the extracellular space enters the cell through voltage gated channels, where it further causes the release of Ca 2+ from the sarcoplasmic reticulum (SR). This release results in a sharp increase of the ICC available to bind to a multi-protein structure known as the troponin complex. This process causes the troponin complex to detach and no longer block the binding sites on a long protein strand called the actin filament. These unblocked binding sites allow for crossbridge elements belonging to the protein strand known as the myosin filament to attach. The cross-bridges can change their configuration through the use of adenosine triphosphate forcing the actin and myosin filaments to slide past each other. The actin filaments are attached to transverse structures called Z-discs placed at longitudinal intervals. A chain-like grouping of these systems is in place so that, when the filaments slide past each other, myocyte shortening occurs and tension is generated in this chain. After a myocyte reaches its full tension, the relaxation process begins, in which Ca 2+ travels back to the extracellular space and SR. The ICC decays to its resting value slower than the initial sharp release and the filaments return to their initial configuration thereby the myocyte elongates and the tension is released. Figure 1 illustrates typical behaviour of the ICC and the corresponding muscle contraction-relaxation in mice myocytes [34] for different stimulating frequencies.
Note that the experimental data are normalized so that the peak value of the ICC with 2 Hz stimulating frequency,c, is 1 [-] and the muscle shortening is given in terms of stretch λ a f . Moreover, the ICC transient exhibits regional variations in the ventricles and, naturally, the characteristics of the myocardial contraction vary nonuniformly throughout the ventricles. The mechanical activation of the ventricles occurs in tens of milliseconds upon the local electrical depolarization, namely, a delay between electrical and mechanical activation exists that is often called electromechanical delay (EMD). In a healthy heart, the EMD is shorter at the endocardium compared to the delay at the epicardium and longer near the base than at the apex. Besides, a notable myocardial prestretch is observed at late-depolarized regions due to the contraction of early-depolarized regions, which detains the beginning of myofiber shortening and leads to a longer EMD [23]. Furthermore, the myocardial contraction and fibre strain have longer duration and larger magnitude in the endocardium compared to the epicardium [3,4]. In addition, it has been reported that  [34] in mice myocytes. In the diagrams, a the calcium transients and b the corresponding cell shortenings are shown for different pacing frequencies (2 Hz, 4 Hz and 7 Hz) which increase along the arrow direction a b timing of muscle relaxation onset and its propagation direction differ at the inner and outer layers of the ventricle [46]. Such variations have been observed in many species and are likely to exist in order to provide an efficient pumping function of the ventricles.

Modelling approaches for myocardial contraction
Among heart related modelling research studies, myocardial contraction is perhaps one of the most investigated topics. In this context, Hill's study underlies most of the modelling attempts for skeletal and myocardial muscles [25]. He experimentally studied the energy liberation of isolated skeletal muscles in terms of heat. Based on these observations, he described the tension-velocity relation through an empirical equation. Moreover, he suggested a rheological setting in which the passive response is represented by a purely elastic spring in parallel connected to a branch having a contractile element and an elastic spring that imitates the active response (contraction) of the muscle. Since then, this rheological setting has been referred to as Hill's three-element model in the literature. Tezören [51] proposed a generalized Hill model for isometric and isotonic phases of the contraction in which the muscle tension is formulated as a function of muscle length, strain rate and activation degree. In a similar way, Pinto [43] suggested an exponential-type equation that relates the muscle force to active contraction and identified the material parameters based on experiments performed on papillary muscles. Apart from empirical modelling attempts, more biophysically based models have been suggested. In this context, the work of Huxley [29] is one of the significant studies motivated by the sliding filament theory (interaction of myosin and actin filaments) of muscle contraction [28]. Inspired by the sliding filament theory, Fung [20] proposed a sinusoidal formulation and adopted it into the three-element Hill model to study the active response of the cardiac muscle. In a similar way, Wong [56,57] furnished the three-element Hill structure with the Huxley model [29], where the role of calcium binding was additionally taken into account. Consideration of intracellular calcium dynamics is another way of modelling the relation between the force generation and sarcomere shortening, which is employed by many authors [22,27,35,44,49]. The biophysically based models have also been utilized for organ scale simulations in various studies [19,30,32,33,45]. For a wide review of biophysically based myocardial contraction studies, the reader is referred to Trayanova et al. [52]. Although such models can to some extent mimic the cellular events during contraction, their application in organ scale simulations is not straightforward. These models consist of a large number of ordinary differential equations (ODEs), which are not straightforward to implement and might cause stability issues [37], and possess many material parameters, which might be difficult to adjust for different simulation conditions. Moreover, large computation requirements can limit the resolution or range of application with current computing capabilities. Alternatively, numerical frameworks employing purely phenomenological approaches for myocardial contraction can provide efficient computer simulations without having any direct physical interpretation of the events occurring on the cellular and tissue level. In that sense, Pelce et al. [42] proposed a practical contraction model which relates the transmembrane potential to the shortening of myocytes. The model comprises a simple evolution equation for the ICC which is governed by the transmembrane potential. The ICC is used as an input to a phenomenological equation to compute the myocyte contraction. Cherubini et al. [9] performed finite element (FE) simulations of the myocardium with a rheological model having an elastic spring and a contractile element connected in serial, where the deformation of the contractile element is described by the model of Pelce et al. Göktepe et al. [21] generalized Hill's classical three-element rheology in a threedimensional setting and utilized the same simple model for the contractile unit in the rheology. Moreover, Nobile et al. [38] modified the Pelce et al. contraction model with an additional term that considers the direct influence of the transmembrane potential on the contraction. Although the equations of Pelce et al. are adequate to represent general characteristics of an excitation-contraction problem, they cannot produce very realistic ICC and myocardial contraction curves. We refer to the work of Göktepe et al. [21] for a detailed parameter sensitivity analysis of the Pelce et al. model.

Aim of the manuscript
In this study, we propose two ODEs that (i) can precisely mimic the general characteristics of experimentally observed ICC and myocardial shortening behaviour, (ii) can be swiftly adjusted to capture different contraction characteristics on the cellular and organ level and (iii) are straightforward to be implemented within a numerical framework providing efficient computer simulations. To this end, one ODE is developed for the ICC having a logarithmic-type evolution and another ODE is related to the myocardial shortening (active stretch) having an exponential-type evolution in a phenomenological sense. Firstly, we calculate the ICC whose argument is specified as the normalized transmembrane potential. Then, the ICC gives rise to the evolution of the active stretch which causes a muscle shortening on the global scale. The model has only seven material parameters in total and both the ICC and active stretch are treated as internal variables and computed in a local sense.
In order to show the applicability of the proposed equations, we perform a detailed material parameter sensitivity analysis and investigate their fitting capacity to different experimental data of the myocardium obtained from already published works in the literature. The equations are further embedded into our fully implicit staggered solution algorithm for cardiac electromechanics [5] in which the modified Hill model is used to describe the electro-visco-active characteristics of the myocardium. Furthermore, several FE simulations in a personalized left ventricle (LV) model are performed pointing out the usability of the new contraction model. In the first part of the numerical examples, we simulate volume-time (v-t) relations of the LV model having different characteristics with respect to the same electrical activity by simply modifying the material parameters in the active stretch equation. We show that the descending and ascending slope and duration at the peak of the v-t curve can be easily altered. In the second part of the numerical examples, we try to mimic the mechanical activation sequence in the light of clinical studies. For this purpose, we consider the heterogeneous characteristics of the contraction throughout the LV by interpolating some of the material parameters related to the active stretch from the endocardium to the epicardium and from the base to the apex. In the simulations on both the material and organ level mentioned above, the phenomenological Aliev-Panfilov model [2] is used for the description of the cardiac electrophysiology. In the last part of the numerical examples, however, we demonstrate that the proposed contraction model can be coupled to biophysically based models of cardiac electrophysiology, e.g. the model of ten Tusscher et al. [50] which will be referred to as the TNNP (abbreviated based on the names of the authors: ten Tusscher-Noble-Noble-Panfilov) model in the rest of the manuscript. To this end, we keep the active stretch equation unchanged and the ICC, which governs the active stretch, is not calculated through the proposed ICC equation, instead, it is generated by the TNNP model. In this setting, we apply a drug called verapamil, which has effects on the electrophysiology and contractility, to the LV model and study its implications in terms of v-t relations and electrocardiograms (ECGs) and compare them to the control (normal) case.
This paper is organized as follows. In Section 2, we introduce the phenomenological equations of the ICC and active stretch and discuss their performance through a detailed parameter sensitivity analysis. The kinematical setting and governing equations of the electromechanical coupling in the myocardial tissue in a monodomain setting are given in Section 3. The constitutive relations are outlined in Section 4. Section 5 demonstrates the applicability of the proposed model through numerical examples which consist of the fitting to experimental data and several initial-boundary value simulations in a personalized LV model are carried out. Finally, we provide concluding remarks in Section 6.

Phenomenological formulation of myocardial contraction
In this section, we introduce our novel phenomenological equations that mimic the myocardial contraction and evolution of its driving ion concentration, i.e. the ICC. Moreover, a detailed parameter sensitivity analysis demonstrating the capacity of the proposed equations is given, where the Aliev-Panfilov model [2], later introduced in Section 4.3.1, is utilized in order to generate the myocyte transmembrane potential with the parameters α = 0.

Intracellular calcium concentration
Motivated by the resemblance of the negative of the logarithmic function to the descending slope of the ICC observed in experiments, we propose the following logarithm-type evolution equation for the normalized ICC,c, where φ is the normalized transmembrane potential value and k 1 , ζ , p 1 are the material parameters. Before starting with the parameter sensitivity analysis, it is worth mentioning that all the parameters intrinsically affect the magnitude of the ICC in Eq. (1). Therefore, to be able to clearly demonstrate features other than adjusting the maximum, the magnitudes of the ICC curves, except in Fig. 2b, are equalized through the parameter ζ . By equalized, it is meant that the magnitude of the peak is normalized to the same value via parameter ζ . Firstly, shifting the peak point of the ICC curve in time without changing the magnitude is illustrated (see Fig. 2a), where k 1 is varied. It is worth noting that the subtracting termck 1 in Eq (1) is the cause for the descending slope of the ICC. Thus, when k 1 is zero, the ICC reaches its peak and stays there as a plateau (not shown). This is the reason why a decrease in k 1 causes the ICC curves to have a wider nature. Moreover, the magnitude of the ICC is changed independently through the parameter ζ , where an increase in ζ moves the magnitude of the peak up without altering the timing of the peak point (see Fig. 2b). Finally, the effects of the parameter p 1 are investigated. Figure 2c demonstrates that decreasing p 1 values lead to a shift of the location of the peak point in time. At first sight, it seems that the parameter p 1 has a similar effect as k 1 . However, the behaviour in the change of the descending limb slope is different when the parameter p 1 is altered. As p 1 is decreased, we observe only a slight flattening of the descending limb and further reductions do not cause the descending limb to form a plateau as k 1 does. Moreover, for p 1 < 0.5 [-], unrealistic ICC curves are produced that almost resemble to the action potential curve (not shown). A summary of the ICC parameters is given in the following.
-Parameter k 1 changes the magnitude and location of the peak in time. This in turn influences the descending and ascending slope. -Parameter ζ changes the magnitude without changing the location of the peak point in time. The descending and ascending slopes are naturally affected due to the change in magnitude of the peak point. -Parameter p 1 changes the magnitude and timing of the peak point with a slighter change in the descending slope compared to the change when parameter k 1 is varied.
Furthermore, the combination of the parameters k 1 and p 1 gives us the chance to alter the descending slope of the ICC curve while keeping the ascending slope and peak point unchanged (see Fig. 2d). It should be further noted that apart from the material parameters controlling the ICC transient, the way the transmembrane potential curve evolves is closely related to the shape of the ICC curve, since the transmembrane potential is considered as the argument of the ICC equation. For instance, a prolongation of the action potential duration without changing the upstroke slope (by increasing the parameter γ ) results in an elevation in the magnitude of the ICC peak (see Fig. 3a), although none of the parameters in the ICC equation are changed. Likewise, a decrease in the upstroke slope of the transmembrane potential curve (by decreasing the parameter b) leads to a shift of the peak point in time and reduction in the magnitude of the ICC curve (see Fig. 3b).

Active stretch
Evolution ofc further gives rise to myocardial contraction represented by the active stretch λ a f that has an exponential-type evolution equatioṅ in terms of the material parameters k 2 , ξ , q and p 2 . Before proceeding with the material parameter sensitivity analysis, we note that as with the parameters for the ICC, all parameters in Eq. (2) naturally influence the magnitude of the active stretch. For that reason, the minimum values of the active stretches, except in Fig. 4b, are brought to the same value through the parameter ξ for a clear illustration of the features other than magnitude adjustability. The parameters for the ICC are determined as k 1 = 0.05 [-], p 1 = 3 [-] and ζ = 0.000222 [-].
Variations in the timing of the minimum value are performed by changing k 2 (see Fig. 4a). For a small k 2 value, the timing of the minimum is delayed and small variations in k 2 can cause a large change in the timing of the minimum. Whereas, when k 2 is large, the timing of the minimum is early in time and k 2 must be largely increased to generate movements towards earlier timings. Another important feature is the ability to independently change the magnitude of the contraction by changing only ξ . This is extremely practical, since it allows ξ to be treated as a parameter that is exclusively tied to the magnitude of the myocardial contraction. Figure 4b displays the movement of the active stretch minimum in magnitude without changing the location in time by a simple adjustment of ξ . In addition, the ascending and descending slopes of the active stretch can be adjusted through q. Figure 4c displays an increase in the slope of both branches of the active stretch along with a shift in time as the parameter q increases. This behaviour can also be interpreted as a prolongation in the duration of the peak contraction. Moreover, the parameter p 2 has a slightly similar effect as the parameter q. As p 2 is decreased in Fig. 4d the slope of ascending limb is reduced and the peak point is shifted in time without a significant alteration in the descending limb slope. Hence, the duration of the active stretch remains almost unchanged in contrast to the effect of the parameter q. A summary of the active stretch parameters is given in the following.
-Parameter k 2 changes the magnitude and location of the minimum in time. This in turn influences the descending and ascending slope. -Parameter ξ changes the magnitude without changing the location of the minimum in time. The descending and ascending slopes are naturally affected due to the change in magnitude of the minimum. -Parameter q changes the magnitude and the ascending and descending slope resulting in a change in the duration and time location of the peak point. -Parameter p 2 predominantly changes the magnitude and the ascending slope while causing a slight change in the the descending slope and the location of the peak point in time.
Furthermore, by a combination of the parameters k 2 and p 2 one can keep the descending slope constant and vary only the ascending slope as demonstrated in Fig. 5a. Another important combination of the parameters k 2 and q is illustrated in Fig. 5b, where a certain amount of increase in the parameters alters the duration of the contraction without changing the location of the peak point.
Moreover, the shape of the ICC (explicitly) and transmembrane potential (implicitly) affects the active stretch evolution. For instance, an increase in the magnitude of the ICC will naturally elevate the magnitude of the active stretch curve (see

Excitation-contraction coupling in myocardial tissue
In this section, we introduce the fundamental variables and kinematics that are utilized to characterize the electro-viscoactive response of the myocardium. Moreover, the coupled field equations governing the deformation and electrophysiology of the myocardium in a monodomain setting are outlined.

Primary field variables and kinematics
The electromechanical state in the myocardial tissue is described by four independent sets of field variables Among the global field variables, the deformation map ϕ(X, t) and the transmembrane potential (X, t) represent the macroscopic motion and the electrical state of the myocardium, respectively. In addition, the micro-motion emerging from the visco-active deformation I va ( , λ a f } and the micro-diffusion I e (X, t) arising from ion concentration differences (e.g. [K] + , [Na] + , [Ca] 2+ ) across the cell membrane are considered as local fields (see Fig. 6). In the numerical framework, the global field variables are discretized spatially whereas the quantities controlling the micro-motion and microdiffusion are treated in a local sense. The myocardial tissue is assumed to be a continuum body B ∈ R 3 and B t ∈ R 3 at initial time t 0 ∈ R and current time t ∈ R, respectively, in which excitation-contraction coupling occurs. The hyperelastic nature of the myocardium is formulated through the nonlinear deformation map that associates material points X ∈ B to spatial points x = ϕ t (X) at time t ∈ R thereby constituting the current state of the body B t . The deformation gradient F := ∇ X ϕ t (X) maps the unit tangent of the reference or the Lagrangian point onto its counterpart in the current or the Eulerian configuration. The non-penetrable deformations are ensured by the constraint J > 0, J := det F . Moreover, the symmetric right Cauchy-Green tensor and left Cauchy-Green tensor are utilized as a strain measure respectively, in terms of the covariant Eulerian metric tensor g ab = δ ab and contravariant Lagrangian metric tensor G AB = δ AB that are used for index raising and lowering purposes in the Cartesian framework. The contractile nature of the myocardium is embodied in the active part of the deformation gradient F a that emerges from the multiplicative decomposition of the total deformation gradient F Furthermore, the mechanical part F m is multiplicatively split into an elastic part F e and a viscous part F v to be able to consider rate-dependent effects arising from the tissue viscosity. The active and viscous deformations of the orthotropic myocardial tissue, which constitute the inelastic part of the deformation gradient, are devised along fibre, sheet and normal directions that are presumed to be orthogonal to each other in the reference configuration. In light of this, the active part of the deformation gradient is prescribed along the fibre direction The multi-field problem: a, b mechanical deformation and c, d electrophysiology of the myocardium. a The macro-motion field ϕ is constrained by essential ϕ =φ and natural boundary conditions t = J −1τ ·n, c the electrical conductance of the action potential is constrained by essential =¯ and natural boundary conditions q = J −1q · n, b micro-motion field I va and d micro-diffusion field I e are locally defined and do not possess any boundary constraints Moreover, the viscous part of the deformation gradient is prescribed in the following orthotropic structure in terms of viscous stretches λ v f , λ v s and λ v n corresponding to three mutually orthogonal orientations f 0 , s 0 and n 0 , respectively. Finally, the prescribed forms of the active and viscous part of the deformation gradient let us attain the explicit expression of the elastic part of the deformation gradient as with the mechanical part of the deformation gradient Note that a fictitious intermediate state between the reference and current configuration is constituted by the rotation free viscous and active part of the deformation gradient in the sense of Cherubini et al. [9] and their orthogonal prescription makes F e insensitive to the order of F v and F a in Eq (6). The stretches along the mutually orthogonal orientations are calculated via the push-forward operation of the Lagrangian unit direction vectors As a result of the orthogonality of the viscous and active part of the deformation gradient, we further multiplicatively split the stretches into active, viscous and elastic parts in the fibre direction and into viscous and elastic parts in the sheet and normal directions, respectively. The rheology for the electro-visco-elastic response of the myocardium and the multiplicative split of the stretches in three mutually orthogonal directions are demonstrated in Fig. 7. Further simplification can be made by transforming the multiplicative decomposition of the stretches to an additive format in the logarithmic strain space as follows This simplified expression of the stretches allows us to construct uncoupled one-dimensional evolution equations for the internal variables, which can be easily updated in a time-discrete setting.

Fundamental equations
In the rest of the manuscript, Cardiac mechanics The quasi-static stress equilibrium is expressed by the balance of linear momentum in its local spatial form in terms of the Kirchhoff stress tensorτ and the volume-specific body forces b in the reference configuration. The transient electrical state in the myocardial tissue leads to a mechanical deformation (i.e. contraction), therefore, the Kirchhoff stress tensor is formulated not only by means of the deformation field but also by taking the electrical field into account The balance of linear momentum (Eq. (15)) is subjected to the prescribed essential boundary condition Additionally, the following natural boundary conditions are introduced with the surface stress traction vector t that is expressed through the Cauchy stress theorem as t := J −1τ · n. The former traction term is prescribed while the latter one is a function of the deformation field and represents the load on the ventricular endocardial surface S v (S v ⊂ S) emerging from ever-changing blood pressure in the ventricular cavity during the cardiac cycle. Note that the expressions S = S ϕ ∪ S t and S ϕ ∩ S t = ∅ have to be satisfied by the surface subdomains.
Cardiac electrophysiology Furthermore, the principal equation of the cardiac electrophysiology in a monodomain setting is expressed aṡ in terms of the flux vectorq describing the propagating electrical waves and the nonlinear source termF φ governing the evolution of the transmembrane potential in time. As in the Kirchhoff stress tensor (Eq. (16)), the flux vector and source term depend on both the electrical and mechanical fields The governing equation of the cardiac electrophysiology (Eq. (19)) is then exposed to the following essential and natural boundary conditions respectively, along with complementary S = S φ ∪ S q and disjoint characteristics S φ ∩ S q = ∅ of the subdomains. The relation between the electrical surface flux term q and the spatial flux vectorq can be written through the Cauchy-type formula q := J −1q · n. Moreover, a specified transmembrane potential value of the material points at an initial state t = t 0 is needed for the solution of the ODE in Eq.

Constitutive equations
This section briefly presents the constitutive equations characterizing the stress responseτ , the ventricular blood pressurê p v , the current sourceF φ and the conductivity tensor D that is necessary to obtain the flux vectorq. In order to formulate the mutual interaction of the mechanical and electrical fields, we equipped the constitutive equations, exceptp v , with additional terms that enable a mechanical deformation resulting from the tissue depolarization, i.e. contraction, and vice versa, a change in the electrical state upon a mechanical stretch, i.e. mechano-electric feedback. All evolution equations are integrated via an implicit Euler backward integration scheme which allows us to use larger time steps compared to frameworks where an explicit or a semi-implicit time integration scheme is used. The reader is referred to Cansız et al. [5] and Wong et al. [58] for a more comprehensive description the constitutive equations, the evolution equations and their update algorithms.

Stress response
The stress response of the living myocardial tissue is formulated through the modified Hill model that incorporates the rheologies proposed by Cansız et al. [8] and Göktepe et al. [21] (see Fig. 7). The significant aspect of the resulting rheological model is the multiplicative split of the deformation gradient into electrical, elastic and viscous parts, while the free energy function is additively decomposed into passive and visco-active parts In the three-dimensional framework, the passive partψ p (g, s; F ) is a function of the total deformation gradient and the visco-active partψ va (g, s; F e ) has a strict dependency on the elastic part of the deformation gradient Moreover, s denotes the set of structural tensors describing the orthotropic material microstructure By making use of the Coleman and Noll method and the Doyle-Ericksen formula, the Kirchhoff stress tensor is written in an additive format with the passiveτ p and visco-activeτ va stress contributions. The well-known orthotropic hyperelastic material model of Holzapfel and Ogden [26] is adopted for the passive part of the free energy function and, accordingly, the passive part of the Kirchhoff stress tensor readŝ by means of the deformation dependent scalar coefficients where a, b, a f , b f , a s , b s , a fs , b fs are material constants describing the deformation state of the isotropic and orthotropic microstructure of the myocardium. For the visco-active part, however, we adopt the simple quadratic free energy function as we previously proposed [5], which leads to the following visco-active part of the Kirchhoff stresseŝ in terms of the nonlinear visco-active Kirchhoff stress coefficients where μ i represents the shear moduli. Moreover, the invariants describing the isotropic and orthotropic deformation state are defined as Evolution equations for internal variables In this part, we briefly refer to the evolution equations that are required to compute the visco-active part of the stress tensor (Eq. (28)). The deformation driven viscous stretches are evolved according to the equatioṅ which is based on the Biot equation, in terms of the material parameters m i , η i andβ i in the fibre, sheet and normal directions. Parameter m i controls the transition from purely linear viscoelastic-to-nonlinear viscoplastic flow. Parameter η i is the viscosity parameter andβ i indicates the activation stress, which is taken as unity to ensure the consistent units.
Note that a closed form solution exists for v i with m i = 0, whereas for m i > 0, v i turns into a nonlinear function of the thermodynamical force β i necessitating the update of the internal variables v i through local Newton-Raphson iterations in a time-discrete setting. Note that the evolution equations forc and λ a f are already given in Eqs. (1) and (2), respectively.

Ventricular blood pressure
Since a quasi-linear relation exists between the ventricular pressure and ventricle cavity volume V during isovolumetric phases and filling, the pressure development in those phases is expressed through a linear evolution laẇ whereκ is a material parameter having distinct values for different phases and imposes the conditionV ≈ 0 during isovolumetric phases andṗ v ≈ 0 during filling. Moreover,θ, θ,Ṽ and V 0 stand for the dilatation at the end of the previous phase, the current dilation, the volume at the end of the previous phase and the initial volume of the ventricle cavity, respectively. For instance,θ andṼ values are measured at the end point of the ejection phase during the isovolumetric relaxation. Furthermore, the pressure development during the ejection phase is described through the three-element Windkessel model where the material parameters C ap , R c and R p determine the resistance and compliance properties of the blood flow and q represents the outflow from the ventricle.

Current source
When we look at the source of the excitation process of the myocardium, it is possible to consider two distinct parts: a purely electrical partF φ e and a mechanical partF where I e is the set of internal variables. The ion exchange between the intracellular and extracellular media, which is exclusively induced by the excitation of nearby cells above a threshold potential value, is described by the purely electrical part, whereas the mechanical part describes the ion exchange through the cell membrane resulting from the stretch of the tissue and is often labelled as the mechano-electric feedback. In the literature, the characterization of the mechanical part is commonly made through the following simple equation where λ f is the stretch along the fibre direction, G s is the maximum conductance and s is the resting potential [41]. Nevertheless, this contribution is only activated if the myocardium is under tension, i.e. λ f > 1. For the purely electrical part, however, two distinct approaches exist: phenomenological and biophysically based models. The former approach requires relatively low computational cost and is capable of reproducing the main features of the transmembrane potential in a myocardial tissue without physical interpretation. In contrast, the latter is based on experimental observations and considers the local evolution of individual ion species at the cost of a demanding computational treatment. In the following, we give a brief introduction of two frequently employed approaches in the context of phenomenological and ionic models.

Phenomenological description: Aliev-Panfilov model
For the phenomenological description of the cardiac electrophysiology, a dimensionless representation is often preferred in the literature by deriving the unitless transmembrane potential φ and unitless time τ as in terms of the scaling factors β φ , β t and the potential difference δ φ . In the Aliev-Panfilov model [2], the source term is expressed aŝ in terms of the material parameters c, α, γ , μ 1 , μ 2 , b and the externally applied stimulus I . The fast depolarization behaviour of the transmembrane potential is formulated through Eq. (37) 1 while the slow repolarization process is described by Eq. (37) 2 in a local sense. Therefore, the so-called recovery variable r is treated as an internal variable in the FE implementation. After the necessary calculations, the physical counterpart of the normalized source term (Eq. (37) 1 ) can be obtained aŝ

Ionic description: Ten Tusscher-Noble-Noble-Panfilov model
The biophysical approach to model the current source considers the dynamics of the main channels, pumps, exchangers and concentrations that influence the flow of ions across the myocyte cell and SR membranes. The flow of ions through the cell membrane leads to the change in the transmembrane potential. One such ionic model for the human ventricular cell was developed by ten Tusscher et al. [50]. The model consists of seventeen ordinary differential equations. One set of the ODEs deals with the gating variables g gate , which are related to how the currents are activated or inactivated. Another set of the ODEs represents the evolution of the concentrations c ion . Therefore, the set of internal variables in the FE implementation is I e = {g gate , c ion }. The evolution and set of the gating variables are expressed aṡ Here, the functions g ∞ gate represent the steady state of the gating variables, while the variables τ gate are related to the time required to reach the steady state [58].
The four concentrations are the free intracellular ion concentrations of sodium (Na + ) c Na , potassium (K + ) c K and calcium (Ca 2+ ) c Ca and the free ion concentration of calcium in the SR (Ca 2+ sr ) c sr Ca . Further, the ionic currents I crt are functions of the transmembrane potential , the gating variables g gate , the concentration variables c ion and the internal variables I e . The ionic current functions along with the set of ionic currents are The currents I crt are defined as the fast sodium current I Na , background sodium current I bNa , sodium potassium pump current I NaK , sodium calcium exchanger current I NaCa , inward rectifier current I K1 , rapid delayed rectifier current I Kr , slow delayed rectifier current I Ks , plateau potassium current I pK , transient outward current I t0 , L-type calcium current I CaL , background calcium current I bCa , plateau calcium current I pCa , leakage current I leak , SR uptake current I up , and the SR release current I rel . Additionally, a stimulus current I can be applied to initiate depolarization. Note that the currents I leak , I up and I rel occur across the SR membrane and not the cell membrane. These currents deal with the storing and release of calcium in the myocyte. Finally, the purely electrical current sourceF φ e can be obtained as the negative of the sum of all the currents that flow across the cell membranê F φ e (I crt ( , I e )) = −(I Na + I bNa + I NaK + I NaCa + I K1 + I Kr + I Ks +I pK + I t0 + I CaL + I bCa + I pCa + I ).
For the sake of brevity, the reader is referred to the work of ten Tusscher et al. [50] and Wong et al. [58] for a more comprehensive discussion of the TNNP model.

Electrical flux
A potential difference between two points in a conductive medium induces an ion flow between those points that is often denoted as electric fluxq and can be expressed aŝ where d iso and d ani are isotropic and anisotropic scalar conductivity parameters, respectively.

Fitting capacity of the proposed model
In order to demonstrate the adequacy and versatility of the proposed equations, we attempt to reproduce several ICC and corresponding cardiomyocyte shortening curves that were obtained from experiments performed on rats [48] and mice [34]. The data concerning the action potential of the cardiomyocytes is not provided in those papers. Therefore, we do not consider real transmembrane potential evolutions which give rise to the ICC and myocardial shortening data in the experiments. Instead, we use reasonable action potential curves that are approximately in line with the duration of the shortening in the experiments. In each of the fittings, the Aliev-Panfilov model is used for the action potential and the excitation is initiated at t = 100 ms by applying a stimulus for 4 ms with I = 1 [-]. The time step is specified as t = 0.5 ms. The ICC data that are extracted from the experiments are shifted in time so that the onset of the ICC in the simulations coincides with the experimental data. For the sake of consistency, the same amount of adjustment in time is applied to the corresponding cardiomyocyte shortening data. The material parameters presented in Table 1 are identified manually without using any optimization procedure. In fact, those parameters can be smoothly determined thanks to the knowledge acquired through the parameter sensitivity analysis in Section 2. The experimental data of Smail et al. [48] obtained from a rat heart and our fits are presented in Fig. 8 for (a) the ICC and (b) the cardiomyocyte shortening. The material parameters of the Aliev-Panfilov model are the same as the ones utilized in Section 2 except for b = 0.2 [-]. The graphs demonstrate a quite accurate fit of the proposed phenomenological equations to the experimental data. The whole ICC data is captured and only slight deviations are observed at the beginning and ending of the myocyte shortening.
Moreover, Lim et al. [34] performed experiments on freshly dissected LV myocytes from young adult mouse hearts with different stimulation frequencies 2, 4, 7 and 9 Hz and measured the ICC along with the corresponding cell shortening. The experimental data and our fits are presented in Fig. 9 with 2 Hz, in Fig. 10 with 4 Hz, in Fig. 11 with 7 Hz and in Fig. 12 with 9 Hz for (a) the ICC and (b) the cardiomyocyte shortening. We do not consider different stimulation frequencies and the material parameters of the Aliev-Panfilov model are the same as the ones utilized in Section 2 except for b = 0.9 [-]. As seen in the figures, we can almost reproduce the entire ICC and muscle shortening data for different frequencies. In some of the fits, a minor inaccuracy appears when the myocyte gets closer to its fully relaxed state. For instance, the slow relaxation in Fig. 9 and rapid relaxation in Fig. 12 cannot be entirely captured by the proposed equations.  [48]. The ICC and cardiomyocyte shortening are shown in a and b, respectively. The ICC data is presented as a normalization to the peak concentration and the cell shortening as the ratio of the cell length to the initial cell length. The initial cell length was assumed to be the mean of the control epicardial myocyte cell length, reported as approximately 130 μm a b

Finite element simulations in a personalized LV model
In this subsection, the performance of the proposed contraction model is demonstrated through several FE simulations in a personalized LV geometry that was generated from 4D echocardiography data [7]. The LV geometry is created at enddiastole and discretized by 18,809 four-node tetrahedral elements over 4116 nodes. The discretization corresponds to an average element size of 4.5 mm. Unless differently specified, a constant time step t = 2 ms is used in the simulations. The same spatial and temporal resolution are used for both the electrophysiology and mechanical problem. The fibre angle is assigned to material points changing linearly from 70 • on the endocardium to −70 • on the epicardium as described by Cansız et al. [6]. Regarding the boundary constraints, linear springs, which impose directional constraints, are attached to the nodes on the basal and epicardial surfaces. The stiffness values of the springs are adopted from our previous work [7] with a slight change in the k z value at the basal surface. Those stiffness values are applied to mimic physiological LV deformation patterns such as rotation, twist and longitudinal shortening. The discretized geometry, fibre orientation and applied boundary constrains are depicted in Fig. 13. Similar to our previous works [5][6][7], before the final results are obtained, subsequent cycles are carried out until the value of the primary field variables becomes almost constant over the subsequent cycles. In all the simulations, the influence of the mechanical deformation on the source term is neglected withF φ m = 0 mV/ms for the sake of simplicity.

Variations in volume-time relations
The v-t relation is a widely used measure by cardiologists, especially because it can be easily obtained through imaging techniques (e.g. echocardiography) in a non-invasive way. These curves indicate valuable data such as the stroke volume (SV) and ejection fraction (EF) that can be used to evaluate cardiac function at first sight. Additionally, the shapes of Fig. 9 Fitting of the proposed phenomenological equations to the experimental data at 2 Hz by Lim et al. [34]. The ICC and cardiomyocyte shortening are shown in a and b, respectively. The data for the ICC is normalized to the maximum value of the 2 Hz graph and shifted downward so that the diastolic ICC is 0. The cell shortening, presented in the form of active stretch, is shown as the ratio of the cell length to the diastolic length a b Fig. 10 Fitting of the proposed phenomenological equations to the experimental data at 4 Hz by Lim et al. [34]. The ICC and cardiomyocyte shortening are shown in a and b, respectively. The data for the ICC is normalized to the maximum value of the 2 Hz graph and shifted downward so that the diastolic ICC is 0. The cell shortening, presented in the form of active stretch, is shown as the ratio of the cell length to the diastolic length a b Fig. 11 Fitting of the proposed phenomenological equations to the experimental data at 7 Hz by Lim et al. [34]. The ICC and cardiomyocyte shortening are shown in a and b, respectively. The data for the ICC is normalized to the maximum value of the 2 Hz graph and shifted downward so that the diastolic ICC is 0. The cell shortening, presented in the form of active stretch, is shown as the ratio of the cell length to the diastolic length a b Fig. 12 Fitting of the proposed phenomenological equations to the experimental data at 9 Hz by Lim et al. [34]. The ICC and cardiomyocyte shortening are shown in a and b, respectively. The data for the ICC is normalized to the maximum value of the 2 Hz graph and shifted downward so that the diastolic ICC is 0. The cell shortening, presented in the form of active stretch, is shown as the ratio of the cell length to the diastolic length a b Fig. 13 Virtualized LV geometry. From left to right: the discretized geometry along with the dimensions (the lengths are in millimetres.), two different views of the fibre orientation and its angle θ distribution and the springs with different stiffness values (in N/mm) attached to the basal and epicardial surface as boundary constraints the curves, namely, magnitudes of the ascending and descending slopes provide important information such as systolic acceleration and filling rate that can be used for a detailed LV function assessment [15]. A study performed by Zeidan et al. [59] presents several v-t curves for both healthy subjects and patients with coronary and hypertensive heart disease. For ill patients, the study found a loss of symmetry in the v-t curves, a lower early filling rate during diastole and a lack of clear distinction between the early filling period from mid-diastolic diastasis when compared to the healthy subjects. Additionally, the researchers observed that the peak filling rates were lower and the time to peak filling rate was significantly prolonged in the patients with coronary and hypertensive heart disease. Moreover, a study by Duvall et al. [15] showed that the values of systolic emptying rates and diastolic filling rates appear to decrease with age.
In the following, we present various v-t relations with respect to the same electrical activity by only modifying the parameters of the active stretch equation in a similar manner as the examples performed at the material level in Section 2.2. Basically, the aim is to demonstrate how the proposed phenomenological model can be utilized to alter the contraction-relaxation (systole-diastole) characteristics on the organ level.
In order to describe the cardiac electrophysiology, the Aliev-Panfilov model is employed. At time t = 0 ms, the myocardial tissue is assumed to be in a resting state with 0 = −80 mV and the depolarization wave is initiated by applying a small amount of stimulus, I = 3 [-] for 16 ms, to the upper part of the septum corresponding to the atrioventricular node. Moreover, the scaling factor β t in the Aliev-Panfilov model is formulated in terms of the activation time t act whereby early excited regions have a longer action potential duration compared to later activated regions [31]. The activation time t act is specified as the elapsed time from the excitation of the atrioventricular node until the transmembrane potential of a particular material point reaches −40 mV. The cycle length is determined as 800 ms. For simplicity, ventricular blood pressure is not taken into account. The necessary material parameters, other than those related to the active stretch, are given in Table 2.
The v-t relations with different parameter sets are presented in Fig. 14 along with the electrical activity which is the same in all examples and illustrated through the ECG graph (most upper diagram). The parameters of the active stretch model are manually calibrated. It is worth noting that the descending and ascending slope and duration at peak point of the v-t curves are mainly determined by the parameters k 2 , q and p 2 , while the parameter ξ is calibrated simply to obtain the same EF value in all curves.
Using a set of parameters which shifts the peak point of the myocardial contraction in time (similar to Fig. 4a), the v-t curves display a comparable movement of the minimum, as seen in Fig. 14 with material set 1. The shifts of the minimum value in time also affect the rate of the ejection (contraction) and filling (relaxation), where a further shift forward in time results in a flatter shape in both the descending and ascending limb of the v-t curve. The same pattern is observed in the myocardial contraction curves in Fig. 4a. Accordingly, the simulations with the parameter set 2 in Fig. 14 demonstrate a similar pattern as that observed in Fig. 5a, where the slope of the myocardial relaxation (filling) is flattened while keeping the rate of the shortening almost the same, thereby the symmetry of the v-t relation is significantly distorted. Moreover, in Fig. 14, the parameter set 3 achieves a notable increase in the duration of the maximum contraction similar to that observed in Fig. 5b, where the duration of the minimum of the myocardial contraction is stretched in time. It is worth noting that the curves generated with material parameter set 4 and 5 do not exhibit a horizontal line at peak value as in Fig. 5b. Observe Table 2 Utilized parameters in the simulations  that the rate of contraction is much steeper for those curves which accordingly results in large viscous stretches. When the peak point is reached, the viscous stretchtes immediately tend to go back their initial state while the active stretch stays as a plateau for a while (similar as in Fig. 5b). Therefore, the LV volume does not stay constant and starts to increase although the active stretch value remains the same for a while. Lastly, we note that the magnitude of the v-t curve can be altered simply with the parameter ξ (not shown) analogously in Fig. 4b.
The changes in the parameters are easily able to capture differences in the v-t relation characteristics, such as the descending and ascending slope, duration and magnitude of the peak, which are related to systolic and diastolic function. The application of this capacity can be quite broad. One benefit would be gained in fitting a v-t relation obtained from an echocardiography or magneto resonance imaging. Once the v-t data is fitted, one can further perform simulations in order to predict, for instance, the pressure-volume (PV) relation, which is measured invasively in reality, or active and passive contributions to the wall stress, which is not possible to acquire in clinical routine. Moreover, the easily adjustable features of the active stretch equation might be useful when one needs to investigate how the characteristics of v-t relation are interrelated with other aspects of the heart, e.g. mechano-electric feedback, phases in the PV relation, boundary conditions etc.

Regional variations in contraction
Experimental and clinical studies show that there is an elegant electrical and mechanical activation sequence in the heart that achieves an efficient pumping function. In a healthy heart, the depolarization wave coming from the atria enters the ventricles through the atrioventricular node and is then rapidly transmitted to the myocardium by the fast conduction system (the His bundle consisting of the left and right bundle branches and the Purkinje fibres) and the whole ventricles are depolarized in a specific order. In general, elapsed time to the electrical activation is shorter on the subendocardial layers and increases towards the epicardium [14].
Upon ventricular depolarization, the myocytes start to contract with distinct regional variations in terms of timing, magnitude and duration of the contraction [46,47]. Experimental observations have revealed that the LV contraction starts with a brief myocyte shortening in the subendocardium, which coincides with a rapid pressure increase during the isovolumetric contraction phase. Since the LV cavity volume remains almost constant in this phase, the subendocardial shortening results in an elongation of the subepicardial fibres. The subepicardial contraction starts with the onset of systolic ejection and the time gap between the beginning of subendocardial and subepicardial activation is related to the duration of isovolumetric contraction. In the ejection phase, however, all parts of the ventricle shorten together. The magnitude of shortening is larger in the apical segments than in the basal segments, which creates an observable wave of shortening travelling from the apex towards the base. On the other hand, during diastole, ventricular relaxation does not follow exactly the same pattern. Firstly, a relaxation wave occurs before the closure of the aortic valve in the subendocardial region that travels from the apex to the base. The subepicardium, however, starts to relax after aortic valve closure and travels in the opposite direction, namely from the base towards the apex. During the isovolumetric contraction, shortening of the subendocardial layers leads to a small amount of clockwise apical rotation and anti-clockwise basal rotation. In the course of ejection, however, the subepicardial layers dominate the rotation because of their longer moment arm and the direction of the rotation is reversed due to the opponent fibre orientation on the layers closer to the epicardium.
In the following, we aim to mimic the aforementioned regional variations of myocardial contraction throughout the virtual LV model with the proposed equations. The Aliev-Panfilov model is used for the cardiac electrophysiology and, in order to simplify the simulation, a simultaneous electrical activation is considered in the LV model by applying I = 0.1 [-] for 20 ms to all elements. Thereby, we neglect the influence of the electrical activation sequence on the ventricular contraction and focus only on the parameters related to the active stretch. The LV wall is divided into two equal segments along its thickness: an inner layer (from the endocardial surface to the middle of the LV wall) and an outer layer (from the middle of the LV wall to the epicardial surface). For each segment a different parameter set is assigned for the active stretch equation to be able to illustrate varying contraction characteristics. The material parameters k 2 and ξ for each layer and q for the outer layer are linearly interpolated from the base to the apex (see Fig. 15). Moreover, constant values q = 0.015 [-] and p 2 = 3 [-] for the inner layer and p 2 = 7.5 [-] for the outer layer are specified. With this setting, we could capture heterogeneous characteristics of the contraction. The rest of the necessary material parameters are adopted from Table 2, except for μ f which is increased by a factor of 3.0. Note that in this simulation the blood pressure in the LV cavity is taken into account. The details of additional considerations regarding the applications of blood pressure in the simulation such as passive inflation and switch conditions between phases of PV relation can be found in our previous works [5,7]. As in the previous example, we consider a cycle length of 800 ms. In the forthcoming results, the cycle, whose results are presented in the following, starts at time t = 100 ms.
In Fig. 16, the PV relation and rotation of the basal and the apical regions along with twist are depicted. The characterization of LV rotation and twist is performed as described by Cansız et al. [7]. Figure 17 illustrates the snapshots of the active stretch λ a f distribution and mechanical activation distribution that is obtained according to λ a f while the values in Fig. 18 are based on the total stretch along the fibre direction λ f . At this point, it is worth mentioning that for the calculation of the current λ f value, the reference state is not taken as the undeformed LV geometry, instead, the deformation state at the end of previous cardiac cycle is considered as the reference configuration for the current cycle since a considerable amount of residual deformation exists at end-diastole upon subsequent cycles. This residual deformation has two sources as observed in our previous works [5,7]. One source is the viscous effects that are considered in the material model. The other one arises from the fact that the value of the active stretch λ a f at some of the regions in the LV does not totally return to its resting value at the end of the cardiac cycle due to the interpolated parameters in the active stretch equation. Moreover, an element is considered to be mechanically activated when the shortening reaches 10% of its maximum [23]. Bear in mind that EMD and elapsed time to mechanical activation are identical since all elements are instantly stimulated at the same time.
In the simulation, physiologically realistic LV contraction with 61% EF and rotation of the basal and apical segments are obtained. Besides, the magnitude of the myocardial shortening (λ a f ) is mostly larger than 70% throughout the LV model. As a result of the applied material parameters in the active stretch equation, which vary throughout the LV geometry, we could reproduce some important characteristics of the ventricular contraction that are mentioned above. The snapshots of λ a f and the corresponding mechanical activation distribution presented in Fig. 17 clearly show that the contraction starts first at the inner (subendocardial) layers at the apex and propagates towards the basal segments. A significant outcome of the  earlier subendocardial contraction is the manifestation of a moderate positive and negative rotation of the basal and apical regions, respectively, as depicted in Fig. 16b. In addition, the beginning of subendocardial activation coincides with the sudden pressure upstroke in the isovolumetric contraction phase. Afterwards, a similar contraction wave, from the apex to the base, is perceived at the outer (subepicardial) layers. Moreover, higher contractility (smaller λ a f values) is present in the regions close to the apex. However, the direction of the basal and apical rotations is reversed as the subepicardial layers start to contract and the magnitude of the rotation is considerably increased. This is undoubtedly a result of the longer moment arm of the subepicardial layers which then dominates the behaviour of the rotation during the ejection. On the other hand, the snapshots of λ a f in Fig. 17 designate that myocardial relaxation first starts at the apex and propagates towards the base at the inner layers, while the outer layers experience the relaxation wave travelling in the opposite direction, from the base to the apex. These results comply well with the observations stated by Sengupta et al. [46].
Nevertheless, if we examine the contraction behaviour based on the total stretch along the fibre direction λ f (see Fig. 18), the results differ from the evaluation based on λ a f . At first glance, the shortening wave propagates with a similar pattern at the inner layers as in Fig. 17. However, as the contraction wave further expands, smaller total stretch (larger shortening) values are measured at the outer basal region although λ a f is larger (less contraction) at these regions, see the snapshots at time t = 230 ms and t = 260 ms in Fig. 18. Accordingly, the mechanical activation distribution differs from the one in Fig. 17 as we move from the apex to the base. Furthermore, the same simulation is performed by neglecting the blood pressure in the LV model to investigate whether the consideration of blood pressure influences the above discussed aspects of the LV contraction. The used parameters are identical except for ξ which is reduced by a factor of 0.78 at every point in the LV model. This is because we would obtain a very high EF in the absence of blood pressure since there would not be any resisting force acting on the endocardium. In the case without blood pressure, the EF is measured as 57%.  Fig. 18 Demonstration of a λ f distributions at different times in the cycle and b mechanical activation distribution calculated according to λ f . In the simulation, LV blood pressure is taken into account a b As one could expect, the characteristics of the mechanical activation based on λ a f are the same (not shown) since the evolution of λ a f is merely governed by the electrical activity. Besides, we observe similar rotation and twist curves with a slight reduction in magnitude (not shown). The essential difference is noticed in the evolution of λ f and corresponding mechanical activation (see Fig. 19). In the absence of the blood pressure, λ f follows almost the same pattern as λ a f and, accordingly, a similar mechanical activation distribution is observed. This difference with respect to the case where the blood pressure is considered can be explained as follows. The presence of blood pressure on the endocardial surface forces myocardial fibres to elongate, particularly, at the inner layers close to the base, which eventually delays the shortening. Therefore, a quite different distribution of λ f and mechanical activation is observed when the pressure in the LV model is taken into account. However, the sequence of the relaxation has similar characteristics in both cases.
On the other hand, the simulation results have discrepancies with some of the previously mentioned clinical findings. In contrast to the observation that the subendocardium has smaller stretch values, i.e. higher contractility, than the subepicardium, we had to specify the parameters related to λ a f so that the outer layers have slightly larger contractility than the inner layers. Such a parameter setting is necessary to be able to obtain mild rotation humps at the beginning of the cycle, which are created by early activation of the subendocardial layers, followed by inverted rotation directions of the inner and outer layers with higher amplitudes (see Fig. 16b). In the simulations where higher contractility is considered for the inner layers, we observe a quite big positive and negative rotation of the basal and apical regions at the onset of the early endocardial contraction, which were far from physiological observations (not shown). Another disagreement was the absence of a significant elongation of the outer layers in the LV model during the isovolumetric contraction phase. Fig. 19 Demonstration of a λ f distributions at different times in the cycle and b mechanical activation distribution calculated according to λ f . In the simulation, LV blood pressure is not taken into account a b

Compatibility to ionic model: drug application
The numerical examples presented here are concerned with the compatibility of the active stretch equation with the TNNP model. For this purpose, we apply a drug called verapamil to the LV model and observe its effects on the cardiac electromechanical behaviour. Verapamil is a calcium channel blocker, which dilates blood vessels, acts as a cardiac depressant and is often used to treat hypertension and chest pain [16,53]. Furthermore, verapamil can suppress arrhythmias which travel through the atrioventricular node. As observed in some in vitro studies, verapamil can abbreviate the action potential duration [1,39]. Regarding the contraction behaviour, calcium channel blockers (such as verapamil) are known to reduce the ICC and muscle contractility [53].
We imitate the effects of drugs on the LV physiology using the approach presented in the work of Costabal et al. [11], where the following expression is employed to block specific ion channels The fractional block β is a fit of a Hill-type equation where C is the concentration of the drug, h is an exponent and I C 50 is the concentration at which 50% of the current is blocked [12]. These parameters can be obtained via a manual patch clamp technique as performed by Crumb et al. [13].
From the same study, we adopt the parameters blocking the rapid delayed rectifier current I Kr as I C 50 = 499 nM, h = 1.1 [-] and the L-type calcium current I CaL as I C 50 = 202 nM, h = 1.1 [-].
To better accommodate the ionic model with the active stretch equation presented in this paper, we modify the original calcium concentration c Ca curve of the TNNP model. We alter the time to peak of the c Ca through a modification in the gate g g (Equation (39) in the publication [58]) by assigning the time variable τ g as 300 ms in the case where g ∞ g > g g and ≥ −60 mV do not apply τ g = ∞ if g ∞ g > g g and ≥ −60 mV 300 ms otherwise. (46) Furthermore, the maximal I up parameter is increased to I max up = 0.00085 mM/ms and the maximal I rel parameter is decreased to I max rel = 0.001399 mM/ms. All in all, the peak of c Ca is shifted further in time and the descending slope becomes steeper as compared to the original evolution of c Ca . Next, regarding the excitation-contraction coupling, the calcium concentration is modified so that it becomes unitless and starts from zero. For this purpose, we subtract the diastolic concentration c Ca,diastolic from the calcium concentration c Ca and divide the value by 1 μM At this point, it is worth mentioning that the ICC, c, that is obtained out of c Ca in the TNNP model will be governing the active stretch λ a f in the following examples. Unless differently specified, the necessary material parameters are given in Table 2 for the stress response, in Table 3 for the TNNP model, specified for human myocytes [54], and its initial conditions are given in Table 4. Here, we note that the LV wall is segmented in three equal parts along the radial direction. The inner, middle and outer segments are used to assign endocardial, M cell and epicardial specific parameters, respectively. Moreover, the parameters related to the conductivity are specified as d iso = 0.05 mm 2 /ms, d ani = 1.0 mm 2 /ms. The diastolic calcium concentration is taken as c Ca,diastolic = 0.03 μM.
Before proceeding with the FE simulations of the LV model, we first investigate the effect of verapamil on a single epicardial cell at the material point level (see Fig. 20). A time step t = 0.1 ms is used for the entire duration of the presented action potential graphs. The parameters of the active stretch are assigned as k 2 = 0.03 [-], ξ = 0.027 [-], q = 0.4 [-] and p 2 = 1.0 [-]. Note that the action potential in Fig. 20 has a different form compared to the original TNNP epicardial myocyte action potential due to the modified c Ca evolution by the time variable τ g in Eq. (46). Higher concentrations of verapamil result in a reduction in the peak and duration of the action potential along with a decrease in the magnitude and time to peak of the ICC c. Accordingly, the active stretch λ a f decreases in magnitude and its minimum occurs earlier in time. Moreover, the relaxation is slowed down as manifested by the milder slope of the active stretch for higher concentrations of verapamil.
For organ scale simulations, we make adjustments in a few maximal current conductances to obtain a meaningful ECG. Due to the discretization used in the simulations, the value of the maximal I Na conductance is decreased to C max Na = 3.56 nS/pF so that the upstroke velocity is proper enough for a satisfactory conduction velocity. Also, to obtain a positive T wave, it is important to take into account that earlier activated tissue has a longer action potential duration compared to later activated tissue [18]. Motivated by the work of Kotikanyadanam et al. [31], we interpolate the maximal I Kr and g m,0 = 0, g h,0 = 0.75, g j,0 = 0.75, g xr1,0 = 0, g xr2,0 = 0, g xs,0 = 0, g r,0 = 0, g s,0 = 1, g d,0 = 0, g f,0 = 1, g f Ca,0 = 1, g g,0 = 0.92 I Ks conductances (C max Kr and C max Ks , respectively) with respect to the electrical activation time t act , which is defined as in Section 5.2.1, through in terms of the maximal current conductance C max crt and the interpolation factors m int = 0.0217 ms −1 and b int = 0.1 [-]. These parameters are interpolated based on the activation time to obtain a heterogeneous distribution of the action potential durations resulting in a positive T wave in the ECG. In each of the FE simulations, we use a smaller time step t = 0.1 ms for the first 90 ms of the cycles to ensure stability during the sharp myocyte depolarization. For the rest of the simulation, a larger time step t = 2 ms is used since the slope of the action potential in time is less steep. For the chosen time steps, we obtain a quadratic convergence behaviour of the global residual norm through Newton-Raphson iterations and do not encounter any stability problems. A cycle length of 900 ms is considered. The parameters used for the active stretch were k 2 = 0.0072 [-], ξ = 0.00728 [-], q = 0.36 and p 2 = 1.0 [-]. Regarding the initiation of the LV depolarization at the atrioventricular node, a larger stimulus is required for higher concentrations of verapamil. Therefore, we applied a stimulus of I = −100 pA/pF, I = −150 pA/pF and I = −300 pA/pF for the cases of control, c verap. = 150 nM and c verap. = 1000 nM, respectively. Furthermore, to reduce the computational load, the blood pressure in the LV is not considered.
The simulations presented in Fig. 21 show the electromechanical behaviour in the control case and in the cases, where two different concentrations of verapamil are applied. The control simulation has a normal ECG with a satisfactory duration of the QRS and QT intervals. Additionally, an EF of 66% is measured. These markers indicate that the LV behaviour is within a physiological range for the electrophysiology and the mechanics. On the other hand, when a verapamil concentration of 150 nM is applied, repolarization occurred earlier and more uniformly, which leads to a smaller QT interval and a reduction in the magnitude of the T wave. Besides, verapamil suppresses the contractility, i.e. the active stretch, resulting in a reduced EF (54%) and a decrease in the duration of systole.
Moreover, the largest concentration of verapamil with 1000 nM makes the repolarization and end systole occur much earlier. The QT interval decreases severely and the T wave is inverted. The T wave inversion agrees with the snapshot at time t = 210 ms, where the repolarization wave propagates in a similar direction to the depolarization wave. In addition, a low EF of 28% is obtained along with a much reduced systole duration. Furthermore, the duration of the QRS complex remains constant, while its magnitude decreases for higher concentrations of verapamil. The tendency observed in the ECGs and v-t relations can be explained by the action potential and active stretch behaviour seen in Fig. 20. The action potential duration is reduced at higher concentrations of verapamil. This causes an abbreviated QT interval and also affects the action potential duration dispersion leading to a smaller and inverted T wave. Additionally, verapamil lowers the magnitude of the active stretch resulting in a smaller EF. The results from the single cell and organ scale simulations show a qualitative agreement with observations from other studies. Action potential durations and QT intervals have been reported to either decrease [1,17,39] or stay the same [24,55] when verapamil is applied. In a recent in silico study which used the O'Hara-Rudy action potential model [40], Costabal et al. [10] found no significant change in the QT interval when verapamil was applied. However, a decrease in both action potential magnitude and duration were also observed in another study [60] which used a more recent version of the TNNP model [54]. Therefore, when verapamil is applied, the changes in the QT interval and action potential are linked with the action potential model. Regarding the contraction behaviour, calcium channel blockers are known to reduce the ICC and muscle contractility [53]. Also, a decrease in the magnitude of the peak tension and the time to peak when verapamil is applied to dog trabecular muscles has been observed [36].

Concluding remarks
In this work, we have proposed a phenomenologically based myocardial contraction model that comprises only two simple ODEs and seven material parameters in total. The main goal was to accurately reproduce the myocardial contraction behaviour on both cell and organ scale simulations with as little effort as possible. The model primarily constitutes the normalized ICC through a logarithmic-type relation governed by the normalized transmembrane potential. Then, the active stretch that represents the myocardial shortening is computed by an exponential-type relation having the ICC as argument. The developed contraction model was further integrated into our fully implicit staggered solution scheme [7] for cardiac electromechanics, where the ICC and active stretch are treated as internal variables. In the framework, the myocardium was considered as an orthotropic electro-visco-elastic material through the modified Hill model. For the description of cardiac electrophysiology, however, either the phenomenological Aliev-Panfilov model or the biophysically based TNNP model was employed depending on the nature of the problem.
The foremost attractive aspect of the proposed model is its easy calibration to represent different curve characteristics. This aspect of the model was demonstrated through a comprehensive material sensitivity analysis at the material level in which we could easily alter the upstroke and descending slope and magnitude of the ICC curve and the ascending and descending slope, magnitude and duration of the active stretch curve. These abilities further enabled us to fit the model to various experimental data with an excellent accuracy. However, we did not consider real action potential data recorded in the experiments [34,48] because it was not provided. Instead, when fitting the ICC and active stretch curves, we employed a reasonable action potential evolution and neglected different stimulation frequencies, which can be addressed as a limitation. Moreover, we discussed the dependency of the ICC and active stretch curves on different action potential behaviour.
Furthermore, the capacity of the myocardial contraction model was illustrated through several initial-boundary value problems in a virtual LV model that was generated from 4D echocardiography data. Firstly, a set of simulations similar to the material parameter sensitivity analysis at the material level was performed and the ability to produce different contraction characteristics at the organ scale was demonstrated by means of the v-t relation of the virtual LV model. In addition, the heterogeneous nature of the myocardial contraction, such as variations in the onset of mechanical activation, relaxation sequence and magnitude, throughout the LV was mimicked in a different set of simulations. For this purpose, some of the parameters related to the active stretch equation were simply interpolated from the base to the apex and from the endocardium to the epicardium. Incorporation of regional myocardial contraction variations in the LV model eventually lead to realistic rotation graphs of the basal and apical segments. In the last set of examples, the ionic model of ten Tusscher et al. [50], in contrast to the phenomenological description [2] of cardiac electrophysiology in the previous examples, was employed to demonstrate the feasibility of the active stretch equation with a biophysically based model of cardiac electrophysiology. In the example, we imitated a drug application (verapamil) and studied its implications on contractility as well as on the electrophysiology of the LV model. It was shown that the incorporation of the ICC that is produced by the TNNP model, i.e. c Ca , into the proposed active stretch equation enables realistic numerical analysis of the excitation-contraction problem without any numerical issue.
One important drawback in our LV model was the absence of the fast conduction system, i.e. the bundle of His and Purkinje fibres. First of all, this limited the ability to obtain a realistic electrical wave propagation and activation sequence and, thus, for instance, the Q and S waves were not detectable in the ECG in Fig. 14. Besides, heterogeneities in the electrical activity are closely linked to different ICC evolutions and, accordingly, to variations in the myocardial contraction throughout the LV. However, we omitted the effect of heterogeneities in the electrical activity on myocardial contraction. Particularly, in Section 5.2.2, we considered the same electrical activation and action potential duration all over the LV model and captured the heterogeneities in the myocardial contraction through the calibration of the material parameters in the active stretch equation.
Some further limitations need to be addressed. First, the phenomenological description is not intended to reproduce some fundamental experimental observations, such as the Frank-Starling law, the force-frequency relationship or Hill's forcevelocity relation. These experimental observations could be considered by adding further dependencies and relationships in the definition of the excitation-contraction coupling. Second, the predictive capabilities of the proposed approach might be reduced without a direct connection to the underlying microstructural mechanisms controlling excitation-contraction coupling. Third, a diminished accuracy of the organ level simulations might be present, since a relatively coarse discretization in space and time was used.
All in all, the proposed contraction model has a huge application potential in computer simulations due to its easy adjustability to different conditions, ability to accurately capture reality, straightforward implementation and compatibility to frameworks having a phenomenological or ionic description of cardiac electrophysiology.