## Abstract

Computational modeling of cardiac function has gradually progressed during the past four decades and now beginning to translate toward clinical use as a noninvasive mean of optimizing cardiac treatment options. Recent experimental findings and numerical investigations have suggested an important role of mechanical and intrinsic properties of cardiac tissues in overall electromechanical dynamics of the heart. The inertial effects, which were usually neglected in earlier computational studies, have now been found to alter cardiac dynamics through stretch activated channels (SAC) and can lead to cardiac disorders under specific set of physiological conditions. Considering possible role of inertia in cardiac dynamics, we have modeled electromechanical activity of the heart with inertia terms for computing pressure volume relation and action potentials over a complete cardiac cycle. To this end, we use the continuum balance laws to capture physiological function of the human left ventricle (LV) on an idealized geometry and solve the resulting equations using a python-based finite element platform. For the same set of pressure boundary conditions, the finite element models for quasi-static (less inertia) and dynamic (with inertia terms) formulation yielded a difference of 4.2% end diastolic volume (EDV), 3.1% ejection fraction, and variations in fiber strain pattern. The mechano-electric transduction channels sensitive to small mechanical perturbations in combination with changes in electrical conductivity due to deformation caused quantitative variations over cardiac electrical activity up to 2.75–5% reduction in action potential duration (APD) at 50% repolarization (APD50) and 3.5–5.75% reduction at 90% repolarization (APD90). Catering the effect of inertia can help the research community to improve future computational models in investigating the electromechanics of the heart.

## 1 Introduction

Cardiovascular disease is the leading cause of death, resulting in 17.5 × 10^{6} global deaths annually [1]. It has a huge impact in terms of mortality, morbidity, and health care cost. Accordingly, substantial effort is being put into understanding the underlying mechanisms of the complex cardiac function. With the increase in computational power, the field of cardiac computational modeling has matured in both scope and methodology. Mathematical models of the heart can now be used to investigate conditions that are otherwise impossible to artificially create in a laboratory setting. In this way, cardiac modeling can be used to simulate a healthy heart and draw comparisons with computational model of an ischemic one. In addition, these computer simulations can also be used to explore the effect of drugs, defibrillation, and pacemakers, in order to optimize sophisticated treatment options at a patient specific level.

To better understand excitation–contraction behavior of the heart, a large number of cardiac electromechanics models have been developed. However, the role of inertial forces has not been investigated in detail as yet. While quasi-static assumption allows the ease of numerical scheme, the actual physiological pump function cannot be assumed as inertia independent. Although some studies include the influence of mass in simulations of the heart, a large majority neglects its influence [2]. In this work, we model the cardiac cycle of left ventricle (LV) with inertia term in the balance of linear momentum which is used to model the cardiac mechanics. To this end, we utilize a python-based finite element framework FEniCS 1.7.0 [3] for modeling the contraction–relaxation cycle of an idealized left ventricle using an invariant-based hyperelastic constitutive relation. Our fully coupled model includes both way coupling, i.e., excitation–contraction as well as mechano-electric feedback, so that the mechanical effects of inertia can be translated into corresponding changes on electrophysiology. FEniCS is a collection of free, open source, software components for solving partial–differential equations that can be used with high-level scripting in python. To the best of our knowledge, this is the first study to model complete cardiac cycle of fully coupled electromechanical activity with inertial effects using an orthotropic constitutive law [4].

## 2 Computational Models of Heart Myocardium Tissue

Diverse constitutive models have been developed for computational modeling of ventricles during the past four decades [5–11]. A successful constitutive model of cardiac myocardium elasticity was proposed by Holzapfel and Ogden in 2009 [4]. The model is constructed in terms of invariants which depict structural properties of the myocardium and treats the material as orthotropic. In 2011, Göktepe et al. [12] adopted this constitutive law and found an optimum parameter set that fits shear experiment of Dokos et al. [13] conducted in 2002. Furthermore, they performed a finite element simulation of the human heart model having consistent myocardium fiber orientation and investigated stress distributions on a generic heart model. In a recent study, Gao et al. [14] fitted parameters of Holzapfel Ogden's constitutive law to a deforming model of human left ventricle myocardium which was constructed using cardiac magnetic resonance imaging. Subtle changes in the myocardium deformation behavior have been reported as both substrate and trigger for arrythymogenesis in a number of studies [15,16] due to the presence of stretch activated channels (SAC). The only study which addresses the effects of inertia on coupled electromechanics of the heart was conducted by Costable et al. [2]. The results of this study suggest 8.12% increase in conduction velocity due to the effect of inertia. Although this is the first study which attempted to quantify effects of inertia on electrophysiology, it did not cater realistic geometry, near physiological fiber orientations, and more importantly, the four physiological phases of cardiac activity.

In this work, we model the complete cardiac cycle of the heart beat using quasi-static and dynamic formulation on the left ventricle geometry inspired from the actual human left ventricle near physiological dimensions (see Fig. 1). Fiber orientations of the cardiac tissues are modeled using a rule-based mathematical model proposed by Bayer et al. [17]. Inertial effects are catered by the inclusion of acceleration term in the balance equation employed to model mechanical deformation of the heart tissue. Coupled simulations are performed for both quasi-static and dynamic formulation so as to identify the effect of inertia term over pressure volume relationship, fiber strain, mechano-electric transduction, and action potential duration (APD) in the mathematical model.

## 3 Continuum Model of Left Ventricle Electromechanics

^{9}cells known as cardiac myocytes. These functional cells produce action potential which travels from one cell to the other generating an active force that deforms heart tissue resulting in a contraction–relaxation cycle. Following the ground breaking work of Hodgin and Huxely [18], ion channel models of cardiac cells have reasonably matured during the last five decades. Still, the computational resources, time and numerical schemes required to model a single heart beat comprising all cells of the heart, can render its clinical use impractical. Hence, we resort to the continuum theory where volume-averaged quantities are used to approximate the observed physical behavior of materials within the acceptable margin of error. To model the deformation of left ventricle cardiac muscle using continuum mechanics, we begin with the laws of classical mechanics. From Cauchy's equation of motion

Equation (1) describes the balance of forces so that the change in linear momentum for a region Ω equals the external forces acting on that region plus the body forces (*B*) and **U** is the displacement that a material point in region Ω experiences due to some external force. Here, *ρ* is the per unit volume mass of myocardium tissue of left ventricle. The magnitude of *ρ* is reported in the literature as 1.1 gm/cm^{3} and the same is also used in earlier computational studies [2,12].

**σ**

*) and active (*

_{p}**σ**

*) components. Therefore, we can split the total stress into passive and active stress such that*

_{a}**C**) given as

**f**), sheet (

**s**), and sheet normal (

**n**) directions. Using the definition of deformation gradient (

**F**), a mathematical representation of tissue deformation has been provided in Eq. (3) in which $f,s,n$ and $f0,s0,n0$ represent the fiber directions in the deformed and the undeformed configuration, respectively, such that $f=Ff0,\u2009s=Fs0$, and $n=Fn0$. Holzapfel and Ogden's constitutive law provides the stress–strain relation using eight parameters

*a*,

*a*,

_{f}*a*,

_{s}*b*,

*b*,

_{f}*b*,

_{s}*a*, and

_{fs}*b*based on uniaxial, biaxial, and shear experiments conducted on samples of the left ventricle myocardium. If

_{fs}**a**

_{0}and

**b**

_{0}represent unit vectors along the fiber and sheet directions, respectively, then the invariants of

**C**used in Eq. (3) are given as

Equation (5) is the monodomain model of cardiac electric propagation in which **u** is the membrane potential, and **D** is the conductivity tensor comprising conductance (mm^{2}/ms) in the fiber *d _{f}* = 0.0952, sheet

*d*= 0.0125, and sheet normal

_{s}*d*= 0.0125 directions. Here,

_{n}*i*

_{ion}is the ionic current generated by cardiac cellular model,

*χ*= 140/mm

^{2}is the surface to volume ratio of the cell membrane, and

*C*= 0.01 m

_{m}*μ*F/mm

^{2}is the cell membrane capacitance. Also, $De=\alpha Di$ defines the monodomain assumption of relation between the internal and the external domain conductivity. The magnitudes of these parameters are obtained from the monodomain model based electrophysiology benchmark [19].

*k*, is 0.005 MPa/mV,

_{f}**u**is the membrane potential, and $ur$ is −85.7 mV, which is the resting potential in Tentusscher and Panfilov ionic cellular model [21] for human epicardial cells. In Eq. (7), the switch function $\u03f5(u)$ determines how fast the active fiber tension will change with respect to transmembrane potential

**u**. The parameter values adopted from the fully coupled model [20] are

*ϵ*= 0.1/mV,

_{o}*ϵ*= 1/mV,

_{m}*l*=

*1/mV, and $ua$ = 0 mV. The effect of mechanical deformation on the electric activity is modeled using stretch activated channels that generate mechano-electric currents. Fiber strain opens the stretch activated channels resulting in the generation of an additional component of current*

*i*

_{sac}in addition to ionic transmembrane current

*i*

_{ion}. Different formulations [15] are reported in the literature for incorporating the effect of deformation induced currents in fully coupled models. In this work, we model mechano-electric effect using Eq. (8) by resorting to a formulation used in Refs. [20] and [23]

*G*(0.1) is the maximum conductance, $us$ (20 mV) is the resting potential of the stretch activated channel, and $|b|$ is the stretch at a given material point. Range of constants used in Eq. (8) are obtained from the coupled formulation used in Ref. [20]. As reported in Refs. [15] and [16], mechano-electric transduction through stretch activated channels is only sensitive to acute stretch and not to contraction hence constant

_{s}*δ*acts as a switch to control myofiber stiffness as outlined in the following equation:

## 4 Simulation of a Complete Cardiac Cycle

To generate near physiological results, we adopt the shape of a thickwalled, truncated, axisymmetric ellipsoid as our model of the left ventricle [24]. A large number of software are available for the generation of geometry and meshing. In our work, the left ventricle geometry has been created with open source meshing tool Gmsh. If the complete left ventricle domain (see Fig. 1) is represented by Ω* _{o}*, then the outer boundary of the left ventricle, known as the epicardium, is denoted as $\u2202\Omega epi$, the inner boundary, referred to as the endocardium, by $\u2202\Omega endo$, and the base as $\u2202\Omega base$. The mesh consists of tetrahedra, and to approximate the solution, we use the piecewise Lagrange polynomials.

### 4.1 Boundary Value Problem.

Equations (11)–(16) constitute the mathematical representation of the dynamic boundary value problem in which *T* is the traction force, *p* is the pressure applied by blood on the endocardium, and **N** is the direction normal. We apply zero displacement (**U**) boundary condition at $\u2202\Omega base$, zero traction (*T*) at the outer boundary of left ventricle $\u2202\Omega epi$, and pressure boundary condition on the inner boundary ($\u2202\Omega endo$).

### 4.2 Numerical Scheme.

The details of the steps involved in the numerical scheme are explained as follows:

*Step 1*. Solve the electrophysiology subproblem [25]:

If we assume that $un=u(tn)$ and $sn=s(tn)$ are already known as the initial conditions, where *s* represents the state variables of the ionic cellular model, a single step of the splitting algorithm consists of the following three substeps:

*θ*= 1 corresponding to the Godunov splitting which is first order accurate in time. Rush Larsen scheme is used for solving ordinary differential equation step 1, whereas for the PDE in step 2 we use a time discretization based on a

*θ*-rule. The

*θ*-rule is a commonly used technique for time discretization of PDEs in which the terms in the equation are computed as weighted averages of the values from the current and the next time-step.

*Step 2*. Solve the mechanical deformation subproblem

- (a)
Update the active stress using the membrane potential

**u**(output of the previous step) coming as the input to this subproblem by solving the ordinary differential equation in Eq. (6). Compute SAC currents based on deformation using Eq. (8) and scale conductance using Eq. (10). Update conductance tensor (**D**) and SAC currents (*i*_{sac}). - (b)
Solve the equilibrium Eq. (1) by converting the balance of linear momentum problem to a variational form. Then, solve the system using Newton's iterations at each time-step.

### 4.3 Pressure Calculations.

*p*is computed in four consecutive steps. Phase 1 is passive inflation which is modeled through linear increase of pressure up to 1.2 kPa (end diastolic pressure). This simulates the filling of blood into LV during diastole which lasts for the first 410 ms. Phase 2 entails isovolumic contraction in which volume of the LV is constrained to remain fixed at end diastolic volume (EDV). This constraint models the phase of cardiac cycle when both atrioventricular valves are closed and LV contracts due to electrical stimulus. For this phase, the two-element Windkessel model [26] provides relation between the volume (

*v*) and the pressure (

*p*) by the following equation:

*C*= 0.001 ml/Pa and

_{p}*R*=

*18,000 Pa/ml represent aortic compliance and resistance to outflow of blood. These parameters have been tuned to achieve a near physiological evolution of pressure. The isovolumic step during phase 2 is calculated iteratively for each time-step as described in Ref. [27] to keep volume (*

*v*) in Windkessel model using the following equation:

After isovolumic contraction, ejection is modeled during phase 3. Instead of keeping the volume constant, we allow it to vary using the two-element Windkessel model describing the outflow of blood as a function of pressure and volume. After ejection phase, the isovolumic relaxation (phase 4) is modeled similar to isovolumic contraction by fixing volume of the LV to end systolic volume (49 ml).

## 5 Numerical Accuracy and Convergence

*E*is measured in

_{n}*L*

_{2}and

*L*error norms as $\u03f5=\Vert x\u2212xref\Vert L2$ and $\u03f5\u221e=\Vert x\u2212xref\Vert L\u221e$, respectively, the rate of convergence (

_{∞}*r*) is calculated after each refinement using Eq. (24), and

_{n}*n*is the level of refinement

## 6 Results and Discussion

### 6.1 Contraction Deformation Cycle.

Left ventricle undergoes extensive geometry changes during the cardiac cycle. Computed displacement field provides a measure to assess contraction–relaxation pattern of myocardial wall. To visualize that our model captures these shape changes, we include a time series of the computed geometry. The results featured in Fig. 2 depict major trends in geometry change during the cardiac cycle. The first 410 ms show passive filling or the diastolic phase. Thereafter, the ventricle undergoes isovolumetric contraction under the influence of electrical activity. After isovolumetric contraction, the maximum geometry change can be seen during ejection in which ventricle contracts up to maximum pressure and then relaxes back to original volume. We use these geometric patterns to study strains at different portions of the ventricle in Sec. 6.2. Myocardium wall thickening can be clearly seen in the computed geometry between 300 and 750 ms where equatorial portion of the myocardium wall thickens in wake of incompressibility.

### 6.2 Fiber Strain.

The time series plotted in Fig. 3 show a comparison of strain pattern between the quasi-static and dynamic models at the base, apex, and mid portion of the left ventricle. End diastolic strain is marked in each plot to show the difference in magnitude of strain between quasi-static and dynamic models. The increase of 2–6% fiber strain at the base and mid portion of the ventricle can be observed in dynamic model as compared to quasi-static model as shown in Fig. 3. As compared to the base and mid portion, larger gradients can be seen at the apex of the ventricle which is in agreement with physiological facts and the geometric changes depicted in Fig. 2. Acute stretch during different phases of the cardiac cycle can possibly alter the cardiac dynamics by activating stretch activated channels and inducing mechano-electric feedback. The effect of these strains results in the generation of mechano-electric currents during diastolic phase and changes in the electrical conductivity tensor in systolic phase when electric activity starts to propagate throughout the myocardium causing active contraction.

### 6.3 Action Potential Duration.

Cardiac rythm and function is sensitive to APD as well as the shape of action potential. Cardiac myocyte fiber stretch amplitude and timing shows, both experimentally and computationally, to effect action potentials through mechanotransduction in a number of ways [5]. In order to quantify the inertial effects, we perform computations of APD at different portions of LV during the cardiac cycle for quasi-static and dynamic formulations. At base of LV, the inertial effects resulted in a 5.75% reduction in action potential at APD90 and 5% reduction in action potential at APD50 (See Fig. 3). At the mid portion, these changes are comparatively moderate resulting in 3.5 and 2.75% reduction at APD90 and APD50, respectively, whereas at the apex, the changes in APD are almost negligible. These results show that the inertial effects causing variation in strain patterns and mechanical deformation possess the capability to perturb electric activation of the myocardium by altering the action potentials in different portions. These perturbations can vary point to point in LV and also result in changes in the conduction velocity since opening and closing timings of ionic channels as well as other vital indices, such as refractory period, are reflected through the action potential shape and duration.

### 6.4 Pressure–Volume Relationship.

Pressure–volume curves are useful indicators of the health of cardiac pumping function. Simulating one cardiac cycle with and without inertia term in the coupled model yield results shown in Figs. 4 and 5. For quasi-static and dynamic model, the left ventricle undergoes a volume change of 67 ml and 70 ml, respectively, which corresponds to an ejection fraction of 55.8% and 58.3%. Computed values of ventricular volume and ejection fraction for both quasi-static and dynamic model relate closely to the reported normal human physiological range. By looking at the first 400 ms of volume curves, the inertial effects due to the pressure on endocardium can be observed for both the quasi-static and the dynamic model (see Fig. 4).

It is observed that the initial conditions in the dynamic model play an important role at this early stage of simulation. The velocities are initially at zero and a rapidly increasing pressure load is applied. At this point, the magnitude of strain is low and only small amount of stress is required to yield significant deformations in the case of dynamic model (up to 100 ms of solid line plot in Fig. 4) due to inertial effects of mass. However, in the case of quasi-static model (dotted line plot in Fig. 4), the deformation occurs rather gradually since the inertial effect on account of mass is not catered in the balance of linear momentum (see Eq. (1)). Due to inertial effects in the dynamic model, increasing pressure instantly changes the volume such that equilibrium is maintained. As the inclusion of mass adds an acceleration term with no form of damping or energy dissipation, the dynamic model overshoots the equilibrium due to inertial effects. This leads to a slight periodic motion in the volume during mid-diastolic phase around 200 ms. The oscillations in the volume curve die out as the deformation increases and the ventricle gets stiffer. Ventricular pressure also undergoes a near physiological change with an end diastolic pressure of 70 mm of Hg (9.332 kPa) and an end systolic pressure of 110 mm of Hg (15.980 kPa). Difference due to inertial effects is not evident in the pressure evolution between the quasi-static and dynamic models.

## 7 Conclusion

We simulated a complete cardiac cycle of the left ventricle using a coupled dynamic electromechanical model that responds to the inertial effects. In this study, displacements are computed during the four physiological phases of contraction relaxation cycle using a finite element computational platform in order to investigate the effects of mass on the ventricle pressure–volume relation, strain pattern, and APD. Our major finding is that additional deformations are caused during mid- and end diastole under the effect of inertia. These deformations bring quantitative changes over the action potentials in terms of shortening of APD. As part of the strain pattern analysis, we also identified the apex as the region of maximum fiber strain at the end diastole as well as during the ejection phase. Further, numerical investigations could refine the effects of electromechanical changes observed in this study. Comparison between the quasi-static and dynamic models presented here will assist the research community to identify effect of inertia in cardiac computational models for exploring complex interactions that govern the cardiac function.

On a broader spectrum, the computational modeling of cardiac system has not fully entered the clinical usage. There is a huge gap to be bridged until organ level cardiac simulations can be used by clinicians as a regular means of practice. The utilization of this study lies in its efficacy to improve future computational models. Furthermore, it can also be used to study the effects of myocardium deformation on the mechano-electric currents in a near physiological manner as compared to earlier studies. Hence, the addition of inertial effects can be a step to bridge the gap between computational models and the actual physiological contraction relaxation of myocardium. This addition can prove to be a small but valuable step toward the translation of computational models to clinical use.

## Appendix

### Appendix A Electrophysiology

Numerical accuracy of electrophysiology solution is established through comparison with FEniCS solution available in supplementary material of N version electrophysiology benchmark [19]. We resort to spatial resolutions of 0.2 and 0.1 mm (edge length using linear tetrahedral elements) and temporal resolutions of 0.05, 0.01, and 0.005 ms. Comparison of activation times of our numerical solution (M1–M8) with benchmark solution (P1–P8) is shown in Table 4. The difference in activation times at the terminal point (P-8) varies up to maximum 3% error for coarse spatial and temporal resolution, i.e., 0.2 mm edge length and 0.05 ms. Results obtained using mesh with spatial resolution of 0.2 mm converged within 2% error as compared to the results obtained with 0.1 mm spatial resolution.

### Appendix B Deformation of an Orthotropic Tissue Cube

*i*=

*1, 2, 3). The deformation can be expressed in terms of stretches*

*λ*as $xi=\lambda jXi$ such that $f=Ff0=[\lambda f,0,0],\u2009s=Fs0=[0,\lambda s,0],\u2009n=Fn0=[0,0,\lambda n]$. We have the Right Cauchy Green deformation tensor

_{j}*λ*, $j=f,s,n$ in the fiber, sheet, and sheet normal directions, the invariants of the Cauchy–Green deformation tensor are given using Eq. (3) as

_{j}*p*is readily determined from the relation $\sigma 33=0$. Using the condition of a volume preserving incompressible material $(J=Det(F)=1)$, and inserting the above invariants in Eq. (3), the following nonlinear system of equations is obtained:

*λ*,

_{f}*λ*, and

_{s}*λ*for a given value of the active stress $(Sa)$ using matlab nonlinear equation solver $fsolve()$ function. For the purpose of this test case, we prescribe an active stress similar to cardiac depolarization/polarization cycle only in the fiber $X1$ direction using a function of the form

_{n}*k*is a constant representing stress rate of change and

*t*is the time elapsed since the moment of electrical stimulation. In order to draw comparison between the analytical and the FE solution, we also solve the same test case as a boundary value problem using finite element method. We use a set of boundary conditions to get uniform deformation, while ensuring uniqueness of the solution. The cube (Fig. 6) can freely deform, and it is subjected to an active stress

*S*resulting in a contraction in the global $X1$, i.e., the fiber direction. The material is assumed to be incompressible orthotropic and, therefore, the cube must expand in the global $X2$ and $X3$ directions. The boundary value problem can be stated as

_{a}$U=U1i+U2j+U3k$ and Γ* _{i}* is a plane defined by $Xi=0$,

*i*=

*1, 2, 3. The deformation in orthotropic cardiac tissue cube is solved using FEniCS with*

*Q*1/

*P*0 finite elements by varying active stress

*S*(Eq. (B5)) and the resulting stretches are compared with analytical solution obtained through matlab function

_{a}*fsolve*(). Comparison of results is depicted in Fig. 7 which shows a good agreement between the analytical result and the numerical solution. So far, we compared the analytical and FE solution of the physical problem while considering orthotropic material model, incompressibility of material, contraction through active stress, etc. Henceforth, we also include acceleration term and use a manufactured solution to verify our code. We assume a unit elastic cube of 1 cm edge length with Saint Venant–Krichhoff constitutive model. To find residual, we insert the exact solution into Eq. (1). Thereafter, we investigate convergence through changing spatial refinement and evaluating the exact, discrete solution over a given time interval, while refining mesh between each simulation. See Table 5 Lagrange elements with degree 2 (

*P*2) are used for displacement,

*N*is the number of elements,

*E*is the error between exact and discrete solution, and

*r*is the convergence rate as defined in Eq. (24). Convergence due to refinement in time is conducted similarly, except that cell size (

*δx*) is held constant and is replaced by

*δt*. Using the results in Tables 5 and 6, it can be proven that the error (

*L*

_{2}norm) is proportional to $O$(

*δ*t) in the case of time discretization and to $O(\delta x2)$ in the case of spatial discretization. The proof of this error estimate lies beyond the scope of this work and may be found in Refs. [30] and [31].