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.
Cardiovascular disease is the leading cause of death, resulting in 17.5 × 106 global deaths annually . 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 . 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  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 .
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 . 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.  adopted this constitutive law and found an optimum parameter set that fits shear experiment of Dokos et al.  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.  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. . 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. . 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
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/cm3 and the same is also used in earlier computational studies [2,12].
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 (mm2/ms) in the fiber df = 0.0952, sheet ds = 0.0125, and sheet normal dn = 0.0125 directions. Here, iion is the ionic current generated by cardiac cellular model, χ = 140/mm2 is the surface to volume ratio of the cell membrane, and Cm = 0.01 mμF/mm2 is the cell membrane capacitance. Also, 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 .
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 . 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 , the inner boundary, referred to as the endocardium, by , and the base as . 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 , zero traction (T) at the outer boundary of left ventricle , and pressure boundary condition on the inner boundary ().
4.2 Numerical Scheme.
The details of the steps involved in the numerical scheme are explained as follows:
Step 1. Solve the electrophysiology subproblem :
If we assume that and 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:
Step 2. Solve the mechanical deformation subproblem
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 (isac).
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.
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
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 . 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.
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 A Electrophysiology
Numerical accuracy of electrophysiology solution is established through comparison with FEniCS solution available in supplementary material of N version electrophysiology benchmark . 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
and Γi is a plane defined by , i = 1, 2, 3. The deformation in orthotropic cardiac tissue cube is solved using FEniCS with Q1/P0 finite elements by varying active stress Sa (Eq. (B5)) and the resulting stretches are compared with analytical solution obtained through matlab function 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 (P2) 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 (L2 norm) is proportional to (δt) in the case of time discretization and to 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.  and .