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 × 106 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 [511]. 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

Electrical propagation in the heart occurs due to coupled activity of almost 2 × 109 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  
B+·σ=ρd2Udt2
(1)

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].

A unique property of cardiac tissue is its active response to electric stimuli. The stress in the case of ventricle myocardium consists of a combination of passive (σp) and active (σa) components. Therefore, we can split the total stress into passive and active stress such that  
σ=σp+σa
(2)
To model the passive part (σp), we refer to Holzapfel and Ogden's constitutive law [4] which provides relation for the Cauchy stress with incompressibility constraint using invariants of the right Cauchy Green deformation tensor (C) given as  
σp=a2bexp[b(I13)]+2af(I4f1)exp[bf(I4f1)2](ff)+2as(I4s1)exp[bs(I4s1)2](ss)+afsI8fsexp(bfsI8fs2)(fs+sf)pI
(3)
Based on experimental observations, the left ventricle myocardium is an orthotropic material with distinct response along the fiber (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,s=Fs0, and n=Fn0. Holzapfel and Ogden's constitutive law provides the stress–strain relation using eight parameters a, af, as, b, bf, bs, afs, and bfs based on uniaxial, biaxial, and shear experiments conducted on samples of the left ventricle myocardium. If a0 and b0 represent unit vectors along the fiber and sheet directions, respectively, then the invariants of C used in Eq. (3) are given as  
I1=tr(C),I3=det(C),I4a=a0·(Ca0),I8ab=a0·(Cb0)
(4)
In this work, we use a parameter set (see Table 1) identified by Gao et al. [14] in a study on human left ventricle myocardium  
χ(Cmut+iion)=α1+α·(Du)
(5)

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, De=α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].

We model the active stress (σa) using formulation suggested in Refs. [12] and [20] for a fully coupled left ventricle study  
σa˙=ϵ(u)[kf(uur)σa]
(6)
 
ϵ(u)=ϵo+(ϵmϵo)exp[exp(l(uua))]
(7)
In Eq. (6), the maximum active fiber tension, kf, is 0.005 MPa/mV, 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 ϵ(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 ϵo = 0.1/mV, ϵm = 1/mV, 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 isac in addition to ionic transmembrane current iion. 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]  
isac=δGs(|b|1)(uus)
(8)
In Eq. (8), Gs (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 δ acts as a switch to control myofiber stiffness as outlined in the following equation:  
δ={1,for|b|>10,otherwise
(9)
In addition to excitation–contraction coupling and mechano-electric feedback, another important element of fully coupled electromechanics lies in scaling of conductance on account of mechanical deformation [22]. In order to cater the effect of mechanical deformations due to inertia on scaling of electrical conductance, we include this effect over conductivity tensor using the following equation:  
D=df(ff)+ds(ss)+dn(nn)
(10)

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 Ωepi, the inner boundary, referred to as the endocardium, by Ωendo, and the base as Ωbase. The mesh consists of tetrahedra, and to approximate the solution, we use the piecewise Lagrange polynomials.

4.1 Boundary Value Problem.

Using Eqs. (1) and (5), dynamic coupled boundary value problem for the left ventricle can be represented as  
·σ=ρd2Udt2inΩo
(11)
 
χ(Cmut+iion)=α1+α·(Du)inΩo
(12)
 
T=p(t)FTNonΩendo
(13)
 
n·(Du)=0onΩepi
(14)
 
U=0onΩbase
(15)
 
T=0onΩepi
(16)

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 Ωbase, zero traction (T) at the outer boundary of left ventricle Ωepi, and pressure boundary condition on the inner boundary (Ωendo).

Equation (12) along with the human epicardial ionic cellular model proposed by Tentusscher and Panfilov [21] and zero flux boundary condition (Eq. (14)) has been used to formulate boundary value problem for the propagation of electric charges throughout the myocardium.

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:

(a) Solve the system  
ut=iion(u,s),u(tn)=un
(17)
 
st=f(u,s),s(tn)=sn
(18)
for tn<t<tn+θΔt. Here, u and s at t=tn+θΔt are represented by uθn and sθn, respectively.
(b) Solve the linear partial differential equation (PDE)  
ut=α1+α·(Diu),u(tn)=uθn
(19)
for tn<t<tn+Δt, with the boundary condition. The resulting u(tn+Δt) is represented by uθn+1
(c) Solve the system  
ut=iion(u,s),u(tn,θΔt)=uθn+1
(20)
 
st=f(u,s),s(tn,θΔt)=sθn
(21)
for tn<t<tn+Δt, to get the solution un+1 and sn+1 at t=tn+Δt. We choose θ = 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 (isac).

  • (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.

The pressure–volume (PV) relation in cardiac ventricular simulations has been modeled to follow the classic pressure–volume loop. LV PV loop is computed by running the simulation with a cycle time of 800 ms (equivalent to 75 beats per minute). The follower-type pressure load 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:  
Cpdpdt+pR=dvdt
(22)
Values of parameters Cp = 0.001 ml/Pa and 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:  
pn+1=pn+(vn+1vn)/Cp
(23)

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

Validating accuracy of numerical computations of cardiac electromechanics with respect to actual physiological function is still an open question. A number of attempts have been made to evolve N-version benchmark problems separately for electrophysiology [19] and cardiac mechanics [28]. However, to date, no consensus solution for fully coupled cardiac electromechanics have been published. To ascertain numerical accuracy of our model, we begin with isolated testing of cardiac mechanics and electrophysiology code first. For cardiac mechanics, we resort to comparison with an analytical test case described in an earlier study [29]. We also use the method of manufactured solution to ascertain code accuracy and convergence. In the case of electrophysiology, the monodomain benchmark problem is solved and the results are compared with the published FEniCS benchmark solution [19]. Both solutions yielded good match as compared to analytical and benchmark problem. Reader is referred to Appendix for further details. In order to compute convergence rates for our fully coupled electromechanics model with inertial effects, we prepare a test case similar to problem 3 (active contraction case) in a recently published benchmark of cardiac mechanics [28]. We inflate the LV up to end diastolic pressure of 1500 KPa and compute error norms for membrane potential and displacement after contraction is triggered under the effect of electrical activity. We use the solution obtained using temporal resolution of 0.001 ms as the manufactured, reference solution. Here, we do not resort comparison with cardiac mechanics benchmark data for problem 3 because of the absence of inertia term in the benchmark solution and obvious differences between material model and activation means. Since no exact solution for the electromechanics problem is available, the error analysis with respect to time steps is carried out by using the reference solution on the same mesh with a spatial resolution of 0.2 mm. En is measured in L2 and L error norms as ϵ=xxrefL2 and ϵ=xxrefL, respectively, the rate of convergence (rn) is calculated after each refinement using Eq. (24), and n is the level of refinement  
rn=log(En/En1)log(Δtn/Δtn1)
(24)

Tables 2 and 3 show that the expected convergence rates are recovered in both cases and the converging rate is indeed coherent with the first-order backward difference scheme under use.

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

To verify the numerical framework, we resort to analytical test case. We consider a unit cube with fibers, sheets, and sheet normals aligned along the global coordinate axes Xi (i =1, 2, 3). The deformation can be expressed in terms of stretches λj as xi=λjXi such that f=Ff0=[λf,0,0],s=Fs0=[0,λs,0],n=Fn0=[0,0,λn]. We have the Right Cauchy Green deformation tensor  
C=(λf2000λs2000λn2)
(B1)
Using stretches λ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  
I1=λf2+λs2+λn2I4f=λf2I4s=λs2I8fs=0
As the unit cube is free to deform, the stress at equilibrium is zero for all components of the Cauchy stress matrix and the Lagrange multiplier p is readily determined from the relation σ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:  
σ11=2ψ1(λf2λn2)+2ψ4fλf2+Sa=0
(B2)
 
σ22=2ψ1(λs2λn2)+2ψ4sλs2=0
(B3)
 
λfλsλn=1
(B4)
In the above equation from Holzapfel Ogden's orthotropic constitutive law, ψ1=exp(b(I13)) 
ψ4f=af(I4f1)exp(bf(I4f1)2)
 
ψ4s=af(I4s1)exp(bs(I4s1)2)
Constants a,b,af,as,bs,bf are model parameters listed in Table 1. The system of Eqs. (B2)(B4) can be solved for λf, λs, and λn 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  
Sa=kt
(B5)
where 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 Sa 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  
·σ=0Ωo,U1=0onΓ1,U2=0onΓ2,U3=0onΓ3
(B6)

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 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 O(δt) in the case of time discretization and to O(δ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].

References

References
1.
Mendis
,
S.
,
Puska
,
P.
,
Norrving
,
B.
, and World Health Organization,
2011
,
Global Atlas on Cardiovascular Disease Prevention and Control
,
World Health Organization
,
Geneva, Switzerland
.
2.
Costabal
,
F. S.
,
Concha
,
F. A.
,
Hurtado
,
D. E.
, and
Kuhl
,
E.
,
2017
, “
The Importance of Mechano-Electrical Feedback and Inertia in Cardiac Electromechanics
,”
Comput. Methods Appl. Mech. Eng.
,
320
, pp.
352
368
.
3.
Alnæs
,
M. S.
,
Blechta
,
J.
,
Hake
,
J.
,
Johansson
,
A.
,
Kehlet
,
B.
,
Logg
,
A.
,
Richardson
,
C.
,
Ring
,
J.
,
Rognes
,
M. E.
, and
Wells
,
G. N.
,
2015
, “
The FEniCS Project Version 1.5
,”
Archive Numer. Software
,
3
(
100
), pp.
9
23
.http://publications.lib.chalmers.se/records/fulltext/228672/local_228672.pdf
4.
Holzapfel
,
G. A.
, and
Ogden
,
R. W.
,
2009
, “
Constitutive Modelling of Passive Myocardium: A Structurally Based Framework for Material Characterization
,”
Philos. Trans. R. Soc. A
,
367
(
1902
), pp.
3445
3475
.
5.
Weise
,
L. D.
, and
Panfilov
,
A. V.
,
2013
, “
Correction: A Discrete Electromechanical Model for Human Cardiac Tissue: Effects of Stretch-Activated Currents and Stretch Conditions on Restitution Properties and Spiral Wave Dynamics
,”
PloS One
,
8
(
3
),
31
(
5
), p.
e59317
.
6.
Demiray
,
H.
,
1976
, “
Stresses in Ventricular Wall
,”
ASME J. Appl. Mech.
,
43
(
2
), pp.
194
197
.
7.
Yin
,
F. C.
,
Strumpf
,
R. K.
,
Chew
,
P. H.
, and
Zeger
,
S. L.
,
1987
, “
Quantification of the Mechanical Properties of Noncontracting Canine Myocardium Under Simultaneous Biaxial Loading
,”
J. Biomech.
,
20
(
6
), pp.
577
589
.
8.
Fung
,
Y. C.
,
1993
, “
Mechanical Properties and Active Remodeling of Blood Vessels
,”
Biomechanics
,
Springer
,
New York
, pp.
321
391
.
9.
Nash
,
M. P.
, and
Hunter
,
P. J.
,
2000
, “
Computational Mechanics of the Heart
,”
J. Elast. Phys. Sci. Solids
,
61
(
1/3
), pp.
113
141
.
10.
Costa
,
K. D.
,
Holmes
,
J. W.
, and
McCulloch
,
A. D.
,
2001
, “
Modelling Cardiac Mechanical Properties in Three Dimensions
,”
Philos. Trans. R. Soc. London. Ser. A
,
359
(
1783
), pp.
1233
1250
.
11.
Bischoff
,
J. E.
,
Arruda
,
E. A.
, and
Grosh
,
K.
,
2002
, “
A Microstructurally Based Orthotropic Hyperelastic Constitutive Law
,”
ASME J. Appl. Mech.
,
69
(
5
), pp.
570
579
.
12.
Göktepe
,
S.
,
Acharya
,
S. N. S.
,
Wong
,
J.
, and
Kuhl
,
E.
,
2011
, “
Computational Modeling of Passive Myocardium
,”
Int. J. Numer. Methods Biomed. Eng.
,
27
(
1
), pp.
1
12
.
13.
Dokos
,
S.
,
Smaill
,
B. H.
,
Young
,
A. A.
, and
LeGrice
,
I. J.
,
2002
, “
Shear Properties of Passive Ventricular Myocardium
,”
Am. J. Physiol.-Heart Circ. Physiol.
,
283
(
6
), pp.
H2650
H2659
.
14.
Gao
,
H.
,
Li
,
W. G.
,
Cai
,
L.
,
Berry
,
C.
, and
Luo
,
X. Y.
,
2015
, “
Parameter Estimation in a Holzapfel–Ogden Law for Healthy Myocardium
,”
J. Eng. Math.
,
95
(
1
), pp.
231
248
.
15.
Reed
,
A.
,
Kohl
,
P.
, and
Peyronnet
,
R.
,
2014
, “
Molecular Candidates for Cardiac Stretch-Activated Ion Channels
,”
Global Cardiol. Sci. Pract.
,
2
, p.
19
.
16.
Li
,
W.
,
Kohl
,
P.
, and
Trayanova
,
N.
,
2004
, “
Induction of Ventricular Arrhythmias Following Mechanical Impact: A Simulation Study in 3D
,”
J. Mol. Histol.
,
35
(
7
), pp.
679
686
.
17.
Bayer
,
J. D.
,
Blake
,
R. C.
,
Plank
,
G.
, and
Trayanova
,
N. A.
,
2012
, “
A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models
,”
Ann. Biomed. Eng.
,
40
(
10
), pp.
2243
2254
.
18.
Hodgkin
,
A. L.
, and
Huxley
,
A. F.
,
1952
, “
A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve
,”
J. Physiol.
,
117
(
4
), pp.
500
544
.
19.
Niederer
,
S. A.
,
Kerfoot
,
E.
,
Benson
,
A. P.
,
Bernabeu
,
M. O.
,
Bernus
,
O.
,
Bradley
,
C.
,
Cherry
,
E. M.
,
Clayton
,
R.
,
Fenton
,
F. H.
,
Garny
,
A.
,
Heidenreich
,
E.
,
Land
,
S.
,
Maleckar
,
M.
,
Pathmanathan
,
P.
,
Plank
,
G.
,
Rodriguez
,
J. F.
,
Roy
,
I.
,
Sachse
,
F. B.
,
Seemann
,
G.
,
Skavhaug
,
O.
, and
Smith
,
N. P.
,
2011
, “
Verification of Cardiac Tissue Electrophysiology Simulators Using an N-Version Benchmark
,”
Philos. Trans. R. Soc. A
,
369
(
1954
), pp.
4331
4351
.
20.
Xia
,
H.
,
Wong
,
K.
, and
Zhao
,
X.
,
2012
, “
A Fully Coupled Model for Electromechanics of the Heart
,”
Comput. Math. Methods Med.
,
2012
, p.
927279
.
21.
Ten Tusscher
,
K. H.
, and
Panfilov
,
A. V.
,
2006
, “
Alternans and Spiral Breakup in a Human Ventricular Tissue Model
,”
Am. J. Physiol.-Heart Circ. Physiol.
,
291
(
3
), pp.
H1088
H1100
.
22.
Quarteroni
,
A. L. F. I. O.
,
Manzoni
,
A.
, and
Vergara
,
C.
,
2017
, “
The Cardiovascular System: Mathematical Modelling, Numerical Algorithms and Clinical Applications
,”
Acta Numer.
,
26
, pp.
365
590
.
23.
Keldermann
,
R. H.
,
Nash
,
M. P.
, and
Panfilov
,
A. V.
,
2007
, “
Pacemakers in a Reaction-Diffusion Mechanics System
,”
J. Stat. Phys.
,
128
(
1–2
), pp.
375
392
.
24.
Peterson
,
D. R.
, and
Bronzino
,
J. D.
, eds.,
2007
,
Biomechanics: Principles and Applications
,
CRC Press
,
New York
, pp.
8
3
.
25.
Sundnes
,
J.
,
Lines
,
G. T.
,
Cai
,
X.
,
Nielsen
,
B. F.
,
Mardal
,
K. A.
, and
Tveito
,
A.
,
2007
,
Computing the Electrical Activity in the Heart
, Vol.
1
,
Springer Science and Business Media
,
Berlin
, pp.
70
78
.
26.
Tsanas
,
A.
,
Goulermas
,
J. Y.
,
Vartela
,
V.
,
Tsiapras
,
D.
,
Theodorakis
,
G.
,
Fisher
,
A. C.
, and
Sfirakis
,
P.
,
2009
, “
The Windkessel Model Revisited: A Qualitative Analysis of the Circulatory System
,”
Med. Eng. Phys.
,
31
(
5
), pp.
581
588
.
27.
Royse
,
C. F.
, and
Royse
,
A. G.
,
2005
, “
The Myocardial and Vascular Effects of Bupivacaine, Levobupivacaine, and Ropivacaine Using Pressure Volume Loops. Anesthesia and
,”
Analgesia
,
101
(
3
), pp.
679
687
.
28.
Land
,
S.
,
Gurev
,
V.
,
Arens
,
S.
,
Augustin
,
C. M.
,
Baron
,
L.
,
Blake
,
R.
,
Bradley
,
C.
,
Castro
,
S.
,
Crozier
,
A.
,
Favino
,
M.
,
Fastl
,
T. E.
,
Fritz
,
T.
,
Gao
,
H.
,
Gizzi
,
A.
,
Griffith
,
B. E.
,
Hurtado
,
D. E.
,
Krause
,
R.
,
Luo
,
X.
,
Nash
,
M. P.
,
Pezzuto
,
S.
,
Plank
,
G.
,
Rossi
,
S.
,
Ruprecht
,
D.
,
Seemann
,
G.
,
Smith
,
N. P.
,
Sundnes
,
J.
,
Rice
,
J. J.
,
Trayanova
,
N.
,
Wang
,
D.
,
Jenny Wang
,
Z.
, and
Niederer
,
S. A.
,
2015
, “
Verification of Cardiac Mechanics Software: Benchmark Problems and Solutions for Testing Active and Passive Material Behaviour
,”
Proc. R. Soc. A
,
471
(
2184
), p.
20150641
.
29.
Eriksson
,
T. S.
,
Prassl
,
A. J.
,
Plank
,
G.
, and
Holzapfel
,
G. A.
,
2013
, “
Influence of Myocardial Fiber/Sheet Orientations on Left Ventricular Mechanical Contraction
,”
Math. Mech. Solids
,
18
(
6
), pp.
592
606
.
30.
Johnson
,
C.
,
2012
,
Numerical Solution of Partial Differential Equations by the Finite Element Method
,
Courier Corporation
,
Chelmsford, MA
.
31.
Langtangen
,
H. P.
,
2013
,
Computational Partial Differential Equations: Numerical Methods and Diffpack Programming
, Vol.
2
,
Springer Science and Business Media
,
Berlin
.