PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Methods Appl Mech Eng. Author manuscript; available in PMC 2014 January 1.
Published in final edited form as:
Comput Methods Appl Mech Eng. 2013 January 1; 253: 323–336.
Published online 2012 July 24. doi:  10.1016/j.cma.2012.07.004
PMCID: PMC3501134
NIHMSID: NIHMS396431

A fully implicit finite element method for bidomain models of cardiac electromechanics

Abstract

We propose a novel, monolithic, and unconditionally stable finite element algorithm for the bidomain-based approach to cardiac electromechanics. We introduce the transmembrane potential, the extracellular potential, and the displacement field as independent variables, and extend the common two-field bidomain formulation of electrophysiology to a three-field formulation of electromechanics. The intrinsic coupling arises from both excitation-induced contraction of cardiac cells and the deformation-induced generation of intra-cellular currents. The coupled reaction-diffusion equations of the electrical problem and the momentum balance of the mechanical problem are recast into their weak forms through a conventional isoparametric Galerkin approach. As a novel aspect, we propose a monolithic approach to solve the governing equations of excitation-contraction coupling in a fully coupled, implicit sense. We demonstrate the consistent linearization of the resulting set of non-linear residual equations. To assess the algorithmic performance, we illustrate characteristic features by means of representative three-dimensional initial-boundary value problems. The proposed algorithm may open new avenues to patient specific therapy design by circumventing stability and convergence issues inherent to conventional staggered solution schemes.

Keywords: electromechanics, bidomain model, finite element method, coupled problems, monolithic, cardiac mechanics

1. Introduction

The past two decades have seen tremendous progress in the computational modeling of the heart [23, 24]. Efficient computational tools for the assistance of patient specific treatment of cardiac disorders is of great scientific and socio-economical interest [3]. We have come to appreciate that these tools can provide access to regionally varying quantities such as wall strains or stresses, which are virtually impossible to measure in the beating heart [30, 58]. Although tremendous effort has been devoted to understand the coupled electrical and mechanical response of the heart, most existing algorithms solve the electrical and mechanical fields in a decoupled way, typically by using different discretization techniques in time and space for the individual fields [27]. Historically, the biochemical response has been modeled by biophysicists [12, 21, 36], the electrical response by electrical engineers [2, 7], the mechanical response by mechanical engineers [8, 20], and the clinical response by clinicians [5, 25]. The lack of cross-talk between the individual disciplines has hampered the creation of a unified, robust and stable, fully coupled multiscale-multifield solution strategy.

The bidomain equations represent a homogenization of the in-tracellular and extracellular medium [33, 59]. The coupled bidomain equations of electrophysiology have been traditionally solved through staggered solution schemes. The solution of the parabolic and elliptic part of the the bidomain equations via operator splitting follows a common recipe: In the first step, the elliptic part is solved for a constant transmembrane potential [var phi]. In the second step, the parabolic part is solved for a constant external potential [var phi]e [56, 61, 63]. Apparently, the operator splitting simplifies the coupled nonsymmetric set of equations with symmetric and smaller subsystems at the expense of computation time and stability [15, 48]. The strong coupling due to steep excitation wavefront causes significant stability issues and renders the staggered algorithms computationally inefficient and expensive. This fact motivated the intense study of the monodomain equation, which is obtained through the proportionality assumption between the conduction tensors of the intra- and extracellular domains of the bidomain model. This assumption reduces the coupled bidomain equations to a single parabolic reaction-diffusion equation, which a priori satisfies the elliptic part of the bidomain model. It should be mentioned that the simplified monodomain model of electrophysiology is incapable of modeling the externally applied shock-induced polarization due to nonproportional conductances between the intra- and extracellular spaces and typically mispredicts the velocity of the front depolarization wave [44]. This limits use of the monodomain model in defibrillation simulations. An alternative remedy for the solution of the bido-main equations is the improvement of the solution algorithms. To this end, semi-implicit integration methods where the local and global field variables are updated through explicit and semi-implicit methods [26], predictor-corrector type time stepping algorithms [61, 63, 54] and three-step operator splitting techniques are employed in order to improve the stability and accuracy of solution algorithms [62]. Various preconditioning strategies, e.g. symmetric successive over relaxation [41], block Jacobi [14, 61], block triangular monodomain [15], multigrid [40, 41, 56, 48] and multilevel (hybrid) Schwarz methods have been devised for the linearized set of equations resulting from operator splitting algorithms [35, 49, 50, 51]. A hybridomain model, which is adapted based on a posteriori error estimator to either bi- or monodomain, is introduced in order to improve the comptutational efficiency of the semi-implicit integration [34].

Various anisotropic continuum based models [22, 37] and finite element formulations [16] are developed for the passive mechanical response of the heart. The three-dimensional structure of the heart can be constructed from segmented MRI data [43, 45] and the myofiber orientation in the heart can be obtained by anisotropic least square filtering followed by fiber and sheet tracking techniques applied to the Diffusion Tensor Magnetic Resonance Imaging data of the excised human heart [46]. Operator splitting schemes based on monodomain electrophysiology have been proposed for the electromechanical excitation-contraction problem in [39, 47, 60, 38]. These approaches combine a finite difference method to integrate the excitation equations, with a Galerkin finite element method to solve the equations governing tissue mechanics.

Current state-of-the-art electrophysiological and passive mechanical models enable researchers to explore arrhythmogenesis at the organ level. In order to enable the use of these models for patient-specific therapy design, we still lack efficient numerical methods. At this point, it is noteworthy to say that the above mentioned procedures do not heal the inherent instability issues associated with both decoupled semi-explicit solution strategies and nonconsistent linearization of the variational formulation. Recently, we have proposed a unified, unconditionally stable, finite element formulation based on a Galerkin-type variational formulation for the monodomain electrophysiology [17], bidomain electrophysiology [9] and monodomain based two-field electromechanics [18]. The unconditional stability and the quadratic convergence of the formulations results from the fully implicit backward Euler scheme employed for the integration of the global and local field variables and the consistent linearization of the residual terms obtained from the weak form.

In this follow up work, we extend the monolithic schemes proposed for the bidomain electrophysiology [9] and the two-field electromechanics [18] to three-field excitation contraction coupling [10]. The proposed structure is inherently modular and can be easily generalized to the three-field electro-chemistry [64], the four-field photo-electro-chemistry [1], or the four-field chemo-electro-mechanics of the heart. Unlike existing discretization schemes, which are most powerful on regular grids [6], the proposed finite element based discretization can be applied to arbitrary geometries with arbitrary initial and boundary conditions. It is easily applicable to medical-image based patient-specific geometries [29, 42, 66]. The resulting algorithm provides an unconditionally stable and geometrically flexible framework, which opens possibilities for the analysis of defibrillation phenomena and their impact on the electromechanical behavior of heart tissue. We demonstrate the performance of the proposed approach through a coupled electromechanical analysis of spiral wave initiation in a slice of cardiac tissue and excitation contraction of a three-dimensional generic biventricular heart model. The computational cost of the model in comparison to monolithic monodomain electromechanics formulation is commented.

This manuscript is organized as follows: In Section 2, we introduce the governing equations of cardiac electromechanics consisting of the bidomain model of electrophysiology and the quasi-static linear momentum balance. Section 3 is devoted to the derivation of the weak forms of the field equations, their consistent linearization, and their spatio-temporal discretization. Therein, we apply the finite element method in space and a backward Euler type finite difference scheme in time. In Section 4, we specify the constitutive equations for the underlying source and flux terms, and derive the corresponding consistent algorithmic tangent moduli. Section 5 is concerned with numerical examples demonstrating the distinctive performance of the proposed approach. We conclude the manuscript with closing remarks in Section 6.

2. Governing equations of the cardiac electromechanics

This section introduces the field equations and corresponding state variables of the coupled boundary value problem of cardiac electromechanics.

2.1. Geometric mappings and the field variables

A body B is a three-dimensional manifold consisting of material points An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg [set membership] B. The motion of the body is defined by a one-parameter function of time via bijective mappings

equation M1
(1)

The point x = An external file that holds a picture, illustration, etc.
Object name is nihms396431ig2.jpg( An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg, t) denotes the configuration of the particle An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg at time t [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig3.jpg. Let the configuration of the material points at a reference time t0 be denoted by X = An external file that holds a picture, illustration, etc.
Object name is nihms396431ig2.jpg( An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg, t0) [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig4.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig5.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg) = An external file that holds a picture, illustration, etc.
Object name is nihms396431ig2.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms396431ig1.jpg, t) denote the placement map for a frozen time frame t. Then, the deformation map equation M2 with

equation M3
(2)

maps the reference configuration X [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig7.jpg of a material point onto the spatial counterpart x [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg. The deformation gradient

equation M4
(3)

maps the unit tangent of the reference or the Lagrangean configuration onto its counterpart in the current or the Eulerian configuration. The gradient operators [nabla]X [•] and [nabla]x[•] denote the spatial derivative with respect to the reference X and current x coordinates, respectively. Moreover, the Jacobian J:= det F > 0 characterizes the volume map of infinitesimal reference volume elements onto the associated spatial volume elements. Furthermore, the reference An external file that holds a picture, illustration, etc.
Object name is nihms396431ig7.jpg and the spatial An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg manifolds are locally furnished with the covariant reference G and current g metric tensors in the neighborhoods An external file that holds a picture, illustration, etc.
Object name is nihms396431ig8.jpg of X and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig9.jpg of x, respectively. These metric tensors are required for the mapping between the co- and contravariant objects in the Lagrangean and Eulerian manifolds [31]. In order to impose the quasi-incompressible nature of the biological tissues, the deformation gradient F is decomposed into volumetric Fvol:= J1/31 and unimodular [F with macron]: = J−1/3 F parts

equation M5
(4)

The ventricular myocardium is represented as a continuum with a hierarchical architecture consisting of discrete interconnected sheets of unidirectionally aligned muscle fibers [4, 11]. We consider the orthotropic nature of cardiac tissue with two spatially varying preferred directions which characterize the local orientation of myofibers and sheets forming the myocardium [22]. We introduce the Lagrangean unit vectors f0 and s0 such that |f0|G = 1 and |s0|G = 1. Under the action of ϕt, the Eulerian counterparts are obtained through the isochoric tangent maps as

equation M6
(5)

respectively.

2.2. Field equations of cardiac electromechanics

The electromechanical state of the heart is governed by three primary field (state) variables

equation M7
(6)

namely the deformation map ϕt, the extracellular potential Φe and the transmembrane potential or the so-called action potential Φ:= Φi − Φe, which is defined as the difference between the intracellular potential Φi and the extracellular potential Φe. The spatio-temporal evolution of the state variables is described by three field equations. The electrophysiology of cardiac tissue is governed by the bidomain equations, which are based on the homogenized geometrical description of the diffusion-reaction equations of the extra- and intracellular media interacting through voltage-gated ion channels and ion exchangers [33, 59]. Contrary to the monodomain model, where the transmembrane potential is the only field variable related to the electrophysiology, the intra- and extracellular potentials are treated as independent field variables in the bidomain model.

2.3. Bidomain equations of cardiac electrophysiology

The electrophysiology of a material point X of cardiac tissue at time t is defined as State(X, t):= {Φi(X, t), Φe(X, t)}. In the absence of externally applied current, the total current flux An external file that holds a picture, illustration, etc.
Object name is nihms396431ig26.jpg of the body is preserved, leading to the following conservation law on the boundary of the subdomain [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig10.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg

equation M8
(7)

where An external file that holds a picture, illustration, etc.
Object name is nihms396431ig11.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms396431ig12.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig13.jpg are the intra- and extracellular current flux vectors and the surface normal vector, respectively. Transformation of the surface integral to a volume integral via the Gauss theorem and the application of the localization theorem yield the local or strong form of (7)

equation M9
(8)

The constitutive relations between the current fluxes and the potentials of the intra- and extracellular media follow the Ohm’s law

equation M10
(9)

Here, An external file that holds a picture, illustration, etc.
Object name is nihms396431ig14.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig15.jpg are the intra- and extracellular conductance tensors, respectively. Moreover, we define the transmembrane current iT as the sum of the ionic and capacitive currents from the intracellular space into the extracellular space such that

equation M11
(10)

Here, ζ denotes the membrane surface to volume ratio in a periodic representative unit of cardiac tissue. An external file that holds a picture, illustration, etc.
Object name is nihms396431ig16.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig17.jpg are the membrane capacitance per unit area and the current density flowing through the membrane ion channels, see Figure 1. The notation [•]:= D[•]/Dt is employed for the material time derivative throughout the manuscript. The transmembrane current iT acts as a source term from the intracellular medium to the extracellular medium, which can be accounted for in the current conservation equations of each subspace as follows

Figure 1
One-dimensional representaion of myocardium and approximation of the cell membrane by resistor-capacitor circuit. The circuit consists of a nonlinear resistor An external file that holds a picture, illustration, etc.
Object name is nihms396431ig23.jpg dependent on the ionic concentrations (sodium [Na]+, potassium [K]+ and leakage currents) ...
equation M12
(11)

Equations (8) and (11) constitute the strong form of the coupled bidomain model of cardiac electrophysiology. The summation of (11)1 and (11)2 defines the conservation of total current, rendering one of two equations redundant. Due to electrophysiological differences between the heart tissue and the surrounding medium, we can distinguish between the two domains, the cardiac domain An external file that holds a picture, illustration, etc.
Object name is nihms396431ig18.jpg and the extracardiac domain An external file that holds a picture, illustration, etc.
Object name is nihms396431ig19.jpg, such that the total domain of interest is An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg = An external file that holds a picture, illustration, etc.
Object name is nihms396431ig18.jpg [union or logical sum] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig19.jpg. The boundary of the cardiac domain, [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig18.jpg = [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig20.jpg [union or logical sum] [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg, consists of the endocardium [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig20.jpg, the boundary of the inner walls, and the epicardium [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg, the boundary of the outer walls, such that [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig20.jpg[partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg = [empty], see Figure 2. Similarly, the boundary of the extracardiac domain, [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig19.jpg = [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg [union or logical sum] [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig22.jpg, consists of the epicardium [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg and the outer boundary [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig22.jpg, such that [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg[partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig22.jpg = [empty]. We can further decompose each boundary into Dirichlet and Neumann-type boundaries, such that equation M13 and equation M14, where a:= {c, e}. In the absence of an extracardiac domain, the epicardium is identical with the boundary of the extracardiac domain, [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg = [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig19.jpg. We can then rewrite the Neumann boundary conditions at the interface [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg as follows

Figure 2
Illustration of the cardiac and extracardiac domains An external file that holds a picture, illustration, etc.
Object name is nihms396431ig18.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig19.jpg and the endocardium [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig20.jpg, the epicardium [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig21.jpg, and outer boundary [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig22.jpg.
equation M15
(12)

The extracardiac medium is represented by a single potential term Φ0 and the corresponding flux vector An external file that holds a picture, illustration, etc.
Object name is nihms396431ig24.jpg = − An external file that holds a picture, illustration, etc.
Object name is nihms396431ig25.jpg[nabla]Φ0, An external file that holds a picture, illustration, etc.
Object name is nihms396431ig25.jpg denoting the extracardiac conductance. It possesses no transmembrane currents, and the set of transient diffusion equation takes the simple form

equation M16
(13)

Equation (12) states that the total current An external file that holds a picture, illustration, etc.
Object name is nihms396431ig26.jpg can be discretized continuously in the intra- and extracardiac medium leading to the following boundary conditions:

equation M17
(14)

where iapp denotes the applied current at the boundary. The second set of equations stem from (11) and is bounded to the myocardium introducing the boundary conditions

equation M18
(15)

such that [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg = [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig27.jpg [union or logical sum] [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig28.jpg and [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig27.jpg[partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig28.jpg = [empty]. The bidomain equations (11)1 and (8) are recast into a more coherent form

equation M19
(16)

equation M20
(17)

in An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg along with the normalized conductances

equation M21
(18)

The bidomain equations are derived for a rigid unit volume of cardiac tissue where the gradient operators with respect to material and spatial setting are identical such that [nabla] = [nabla]x = [nabla]X. For the extension of the bidomain equations to the finite deformation setting, two assumptions are made: The electrical fluxes of Kirchhoff type are assumed to be linearly depending on the positive spatial gradient of the respective potentials

equation M22
(19)

Moreover, the Equations (16, 17) are assumed to be valid in the spatial setting, e.g. divergence operator is taken to be div(•) = [nabla]x · (•) and for unit volume of the deformed configuration.

Modification of the bidomain equations according to the assumptions made above and normalizing the flux terms with the volume change J, we end up with the modified set of equations describing the cardiac electrophysiology of a deforming medium

equation M23
(20)

equation M24
(21)

for a unit undeformed volume in An external file that holds a picture, illustration, etc.
Object name is nihms396431ig7.jpg. Equation (20) is a parabolic partial differential equation and equation (21) defines an elliptic partial differential equation. An external file that holds a picture, illustration, etc.
Object name is nihms396431ig29.jpg is the Fitzhugh-Nagumo type generic function defining the excitation induced ionic currents through the cell membrane from the intracellular to the extracellular domain.

2.4. Cardiac mechanics

The balance of linear momentum

equation M25
(22)

states the quasi-static mechanical equilibrium in terms of the Kirchhoff stress tensor [tau] and a given body force b. The balance of linear momentum (22) is subjected to essential boundary conditions ϕ = [phi with macron] on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig30.jpg and natural boundary conditions t = J−1τ · An external file that holds a picture, illustration, etc.
Object name is nihms396431ig13.jpg on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig31.jpg such that [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg = [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig30.jpg [union or logical sum] [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig31.jpg and [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig30.jpg[partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig31.jpg = [empty]. The momentum balance depends on the transmembrane potential Φ through the functional dependence of the Kirchhoff stress tensor.

2.5. Constitutive equations

The solution of the coupled initial boundary value problem requires the definition of the constitutive relations which characterize the Kirchhoff stress tensor [tau], the potential flux vectors and the current source An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg appearing in Equations (2022). The Kirchhoff stress tensor [tau] can be decomposed into the passive [tau]pas and active [tau]act parts

equation M26
(23)

see [38, 47]. The passive part [tau]pas is solely governed by mechanical deformation, while the active part [tau]act is generated by the excitation-induced contraction of the myocytes during the course of depolarization. Since the formulation is laid out in the Eulerian setting, we need to explicitly include the current metric g among the arguments of the constitutive functions. Since the parabolic and the elliptic parts of the bidomain equations are altered in order to replace the intracellular potential Φi with the transmembrane potential Φ, the following potential flux vectors are defined

equation M27
(24)

through the deformation-dependent anisotropic spatial conduction tensors An external file that holds a picture, illustration, etc.
Object name is nihms396431ig33.jpg(g; F, f0) and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig34.jpg (g; F, f0) governing the conduction speed of the non-planar depolarization fronts of the intra- and extracellular media in three-dimensional anisotropic cardiac tissue. In the monodomain model of cardiac electro-physiology, based on the assumption of coaxiality of the intra-and extracellular conductances, a unique effective monodomain conductance reads An external file that holds a picture, illustration, etc.
Object name is nihms396431ig35.jpg. The elliptic part of the bidomain equations can also be omitted by the coaxiality assumption [9]. The last constitutive relation characterizes the electrical source term of the Fitzhugh-Nagumo-type excitation equation (20)

equation M28
(25)

It is additively decomposed into the excitation-induced purely electrical part equation M29 and the stretch-induced mechano-electrical part equation M30. The former describes the effective current generation due to the inward and outward flow of ions across the cell membrane. This ionic flow is triggered by a perturbation of the resting potential of a cardiac cell beyond some physical threshold upon the arrival of the depolarization front. The latter accounts for the opening of ion channels under the action of deformation [28, 38]. Apart from the primary field variables, the recovery variable r appears among the arguments of equation M31 in (25). It characterizes the repolarization response of the action potential. The evolution of the recovery variable r determines the shape and the duration of the action potential inherent to each cardiac cell and may change throughout the heart. Its evolution is commonly modeled by a local ordinary differential equation

equation M32
(26)

From an algorithmic point of view, the local nature of the evolution equation (26) allows us to treat the recovery variable as an internal variable. This is one of the key features of the proposed formulation that preserves the modular global structure of the field equations as set out in our recent works [17, 18]. Furthermore, as mentioned in Section 1, cardiac tissue possesses an anisotropic and inhomogeneous micro-structure. This undoubtedly necessitates the explicit incorporation of a position-dependent orientation of myocytes, possibly in terms of structural tensors, in the argument list of the constitutive functions for [tau], D and equation M33. At this stage, however, we have suppressed this dependency for the sake of conciseness by leaving details out until Section 4 where we introduce the active and passive response of the cardiac tissue.

Having the field equations and the functional forms of the constitutive equations at hand, we are now in a position to construct a unified finite element framework for the monolithic numerical solution of the strongly coupled problem of three-field cardiac electromechanics.

3. The finite element formulation

In this section, we present a monolithic set of equations derived from a Galerkin type finite element formulation. To this end, we construct the weak forms of the Equations (2022) and consistently linearize them along the three independent field variables, namely the deformation map ϕ (X, t), the transmembrane potential Φ(X, t) and the extracellular potential Φe(X, t). An identical temporal as well as spatial discretization scheme will be employed for the mechanical and electrical fields. All field variables are discretized with isoparametric shape functions to transform the continuous integral equations of the nonlinear weighted residuals and their linearizations to a set of coupled, discrete algebraic equations. This set of algebraic equations is then solved monolithically via a Newton-type iterative solution scheme for the nodal degrees of freedom.

3.1. Weak formulation

A conventional Galerkin procedure is applied to construct the weak forms of the governing field equations (2022). For this purpose, the residual equations are multiplied by the weight functions δϕ, δΦ, and δΦe, which satisfy the essential boundary conditions, δϕ = 0 on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig36.jpg, δΦ = 0 on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig37.jpg, and δΦe = 0 on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig38.jpg and belong to the Hilbert space equation M34 consisting of functions which, along with their weak partial derivatives, are square integrable in the domain An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg. The weighted residual equations are then integrated over the solid volume. The conventional procedures, integration by parts along with the Gauss theorem, lead to the following weighted residual expressions for the balance of linear momentum

equation M35
(27)

the parabolic part of the bidomain equations

equation M36
(28)

and the elliptic part of the bidomain equations

equation M37
(29)

The explicit forms of the individual terms in (27) are defined as

equation M38
(30)

where the body force b and the surface traction t are assumed to be given. The weak residual expressions equation M39 and equation M40 in (28) stem from the parabolic part of the bidomain equations (20) and read

equation M41
(31)

along with the definitions

equation M42
(32)

Likewise, the residual expressions equation M43 and equation M44 in (29) of the elliptic part of the bidomain equations (21) are given by

equation M45
(33)

The surface flux terms equation M46 and equation M47 are the prescribed natural boundary conditions. In contrast to the mechanical external term in (30)2, equation M48 depends explicitly upon the field variables through the non-linear source term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg.

3.2. Consistent linearization

The weighted residual equations (30, 31, 33) are nonlinear functionals of the primary variables due to the geometric and constitutive terms. The solution of these equations thus requires a consistent linearization with respect to the field variables about ϕ = [phi with tilde], Φ = [Phi w/ tilde] and Φe = [Phi w/ tilde]e

equation M49
(34)

The incremental terms ΔGϕ and ΔG[var phi], which can be obtained through the Gâteaux derivative, may be expressed in the following decomposed forms

equation M50
(35)

based on the definitions in (2022). We then start with the elaboration of the increment equation M51 according to (30)1

equation M52
(36)

Linearization of the non-linear terms in (36) yields

equation M53
(37)

where £Δϕ [tau] denotes the objective Lie derivative along the increment Δϕ and can be expressed as

equation M54
(38)

The fourth-order spatial tangent moduli An external file that holds a picture, illustration, etc.
Object name is nihms396431ig39.jpg in (38) and the sensitivity of the Kirchhoff stresses to the action potential Cϕ[var phi] introduced in (37)2 are defined as

equation M55
(39)

respectively. Incorporation of (37, 38) into (36) results in the following expression

equation M56
(40)

The terms in (40) demonstrate the inherent nonlinearities arising from the mechanical material response, the geometry, and the coupled electromechanical stress response. Since the body force b and the traction boundary conditions t in (30)2 are prescribed, we have equation M57 yielding the identity equation M58. Recalling the explicit algorithmic form of equation M59 from (31)1, the increment equation M60 can be expressed as

equation M61
(41)

Analogous to (37)1, linearization of [nabla]x(δΦ) leads to

equation M62
(42)

Furthermore, based on the functional definition of the spatial potential fluxes, An external file that holds a picture, illustration, etc.
Object name is nihms396431ig40.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms396431ig41.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig42.jpg, we obtain

equation M63
(43)

where £Δϕ An external file that holds a picture, illustration, etc.
Object name is nihms396431ig43.jpg denotes the objective Lie derivative of the potential flux An external file that holds a picture, illustration, etc.
Object name is nihms396431ig43.jpg along the increment Δϕ

equation M64
(44)

for a ={i, ie, e}. In Equations (4344), the second-order deformation-dependent conduction tensors

equation M65
(45)

and the third-order mixed moduli

equation M66
(46)

are introduced, respectively. Substituting the results (4243) and (44) into (41), we obtain the following explicit expression for the increment

equation M67
(47)

where C[var phi]ϕ = Ci + Cie. The external term equation M68 in (28) depends non-linearly on the field variables through the source term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg defined in (25). For a given q on [partial differential] An external file that holds a picture, illustration, etc.
Object name is nihms396431ig44.jpg, we then obtain the following incremental form

equation M69
(48)

where

equation M70
(49)

is the linearization of the scalar-valued function An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg. The tangent terms An external file that holds a picture, illustration, etc.
Object name is nihms396431ig45.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig46.jpg in (48)2 are defined as

equation M71
(50)

respectively. Based on the decomposed form introduced before, the scalar tangent term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig46.jpg can be expressed as

equation M72
(51)

Inserting the result (49) into (48), we obtain the linearized external term

equation M73
(52)

Considering the weak form of the elliptic part of the bidomain equations equation M74 in (33), the increment equation M75 can be written as

equation M76
(53)

Substituting the results (42, 43) and (44) into (53), we finally obtain

equation M77
(54)

where C[var phi]eϕ = Ci + Ce. The increment operator Δ [•] in Equation (42) applies to Φe in the same way as Φ. This completes the linearization within the continuous spatial setting. The increment equation M78 vanishes, since the external term equation M79 in (29, 33) depends solely on the prescribed boundary flux term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig47.jpg.

3.3. Spatial discretization

In a final step, we perform the spatial discretization of the field variables to obtain algebraic counterparts of the residual expressions (30, 31, 33) and construct the element matrices from the linearized terms (40, 47, 54). To this end, we discretize the reference domain An external file that holds a picture, illustration, etc.
Object name is nihms396431ig6.jpg into element subdomains An external file that holds a picture, illustration, etc.
Object name is nihms396431ig48.jpg such that equation M80 with nel denoting the number of the subdomains An external file that holds a picture, illustration, etc.
Object name is nihms396431ig48.jpg in the body. We then interpolate the field variables and the associated weight functions on each element domain by introducing the discrete nodal values and An external file that holds a picture, illustration, etc.
Object name is nihms396431ig49.jpg-continuous shape functions

equation M81
(55)

where nen refers to the number of nodes per element. Based on the discretization (55), the spatial gradient of the weight functions and of the incremental fields read

equation M82
(56)

Incorporating the discretized representations (55, 56) in (27, 28, 29) along with (30, 31, 33), we obtain the discrete residual vectors

equation M83
(57)

where the operator A denotes the standard assembly of element contributions at the local element nodes i = 1,…,nen over nel subdomains. Following analogous steps, the discrete form of the linearized residual terms (34) can readily be obtained by substituting the discretized representations (5556) in (40, 47, 52, 54). Alternatively, the linearization of the discrete residual expressions (57)1–3 can be shown as follows

equation M84
(58)

The coupled element matrix is derived via incorporation of (55, 56) into (40, 47, 52, 54)

equation M85
(59)

where the components of element matrix An external file that holds a picture, illustration, etc.
Object name is nihms396431ig50.jpg are defined as

equation M86
(60)

4. Constitutive equations

In this section, we specify the constitutive equations that are utilized in the representative numerical examples in Section 5. However, the generic coupled formulation described in Section 3 is independent of the particular choice of constitutive functions characterizing the physiology of cardiac tissue. We identify the explicit expressions for the Kirchhoff stress [tau], the potential fluxes qi, qe, qie, and the current source An external file that holds a picture, illustration, etc.
Object name is nihms396431ig29.jpg, whose functional dependencies have already been briefly outlined in Section 2.5. These equations include not only the explicit functional evaluations, but also the accompanying ordinary differential equations governing the temporal evolution of the related internal variables. This, in turn, requires the construction of algorithmic procedures for the local update of these internal variables at the quadrature point level. We also illustrate the computation of the consistent tangent moduli, to ensure quadratic convergence of the underlying Newton iteration.

4.1. Active and passive stress response

Cardiac tissue possesses a highly anisotropic micro-structure that is chiefly made up of unevenly distributed myofibers. We characterize the orthotropic tissue microstructure in terms of two characteristic directions, the local orientation of the myofibers and the normal to the sheets forming the myocardium [16, 22]. We model cardiac tissue as hyperelastic material characterized in terms of a free energy function, such that the stresses are derived from the potential as conjugate pairs of the strain measures. To achieve an objective representation of the free energy function independent of the choice of reference frame and superimposed rigid body transformations, we introduce a set of invariants associated with the unimodular part of the right Cauchy-Green tensor and the two unit vectors f0 and s0 [55]

equation M87
(61)

In line with [13, 32, 53], the free energy function is additively decomposed into volumetric and isochoric parts

equation M88
(62)

The volumetric part is used to enforce quasi-incompressiblity and the deviatoric part characterizes the orthotropic isochoric response in terms of Fung-type free energy functions according to [22]

equation M89
(63)

The additive decomposition of the free energy (62) leads to the decomposed stress response through the Doyle-Ericksen formula

equation M90
(64)

in terms of the pressure [p with hat]:= JU′(J), the Kirchhoff stress [tau]:= 2 [partial differential]g [Psi] and the contravariant metric tensor in the Eulerian manifold g[sharp] = g−1. The isochoric part of the Kirchhoff stress tensor τiso is obtained through the contraction of [tau] with the isochoric projection tensor equation M91 where An external file that holds a picture, illustration, etc.
Object name is nihms396431ig51.jpg:= −g g[sharp]is the fourth-order identity tensor. The explicit form of [tau] reads

equation M92
(65)

where the deformation-dependent scalar coefficients [Psi]1, [Psi]4f, [Psi]4s, and [Psi]8fs are defined as follows

equation M93
(66)

The active Kirchhoff stress [tau]act is generated by the excitation-induced contraction of the cardiac muscle cells, this part of the stress tensor is approximated to be of purely anisotropic form

equation M94
(67)

The above equation implies that the direction of the active stress tensor is dictated by the deformed structural tensor f [multiply sign in circle] f, while its magnitude is chiefly determined by the transmembrane potential-dependent active fiber tension σ(Φ). To model the twitch-like response of the fiber tension σ (Φ), we adopt an evolution equation [38]

equation M95
(68)

in which the parameter kσ controls the saturated value of σ for a given potential Φ and a given resting potential Φr, which is about −80 mV for cardiac cells. That is, [sigma with dot above] vanishes identically when σ admits the value σ = kσ (Φ − Φr) for ε(Φ) ≠ 0. Contrary to its Heaviside form proposed in [38], we use the following smoothly varying form for the rate switch function [18]

equation M96
(69)

in terms of the parameters ε0 and ε that characterize the two limiting values of the function for Φ < [Phi w/ macron] and Φ > [Phi w/ macron] about the phase shift [Phi w/ macron], respectively. The transition rate of ε from ε0 to ε about [Phi w/ macron] is determined by the parameter ξ. To compute the current value of σ, we discretize the above equations in time using an implicit backward Euler scheme. The closed-form algorithmic expression for the active fiber tension reads

equation M97
(70)

in terms of the current action potential Φ and the fiber tension σn of the previous time step tn, which is stored as a history variable at each quadrature point of the finite element discretization. Next, we determine the algorithmic tangent moduli based on the definition (38)

equation M98
(71)

The volumetric part is given by the common expression

equation M99
(72)

where [kappa macron]:= J2U″(J) is the bulk modulus. The symmetric fourth identity tensor has the indicial representation equation M100 in terms of the components of the contravariant metric g[sharp]. The isochoric part of the moduli has the following form

equation M101
(73)

where the unimodular tangent term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig52.jpg is derived from the stress expression (65)

equation M102
(74)

in terms of the fourth order structural tensors

equation M103
(75)

Similarly, the sensitivity of the Kirchhoff stress tensor to the transmembrane potential then follows from (39)2

equation M104
(76)

4.2. Spatial potential flux

We assume that the second-order conduction tensors An external file that holds a picture, illustration, etc.
Object name is nihms396431ig53.jpg in the intra- and extracellular domains a:= {i,e} can be additively decomposed

equation M105
(77)

into the conductance along the fiber direction equation M106 and the transverse conductance across the fiber direction equation M107, which can be determined experimentally using microelectrode array recordings [7]. The third-order mixed moduli Ca which is defined in (46) and enter the formulations in (47) and (54) follows from inseting (46) into (32)

equation M108
(78)

4.3. Current source

Last, we specify the constitutive equation for the electrical source term An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg. To set up the model equations and parameters in the non-dimensional space, we introduce the non-dimensional transmembrane potential [var phi] and the non-dimensional time τ through the following conversion formulae

equation M109
(79)

The non-dimensional potential [var phi] is related to the physical transmembrane potential Φ through the factor β[var phi] and the potential difference δ[var phi], both in millivolt (mV). Likewise, the dimension-less time τ is converted to the physical time t by multiplying it with the factor βt in millisecond (ms). We then obtain the following conversion expressions

equation M110
(80)

for the normalized source term f[var phi], and the non-dimensional counterparts h:= [partial differential][var phi] f[var phi] and h:= 2 [partial differential]g f[var phi] of the tangent terms defined in (50). The additive split of An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg, introduced in (25) Section 2.5, along with (80)1 implies the equivalent decomposition of equation M111 into the purely electrical part equation M112 and the stretch-induced mechano-electrical part equation M113. This leads to the dimensionless counterpart of (51)

equation M114
(81)

Throughout the subsequent treatments, the Aliev-Panfilov model [2] is used for the description of the source term equation M115 and the evolution of the recovery variable r, following the algorithmic treatment introduced in [17, 18].

5. Representative numerical examples

In this section, we illustrate the performance of the fully coupled algorithm established in the previous sections. First, we examine an orthotropic block of tissue subjected to an initially smooth excitation wave, which turns into a reentrant spiral wave after a slight perturbation. Second, we analyze the excitation of a generic biventricular heart to illustrate the overall electro-physiological and electro-mechanical features of the excitation-contraction coupling in the heart. The convergence behaviour of the proposed algorithm will be discussed for each example. The stretch induced mechano-electric feedback is omitted throughout the subsequent treatments.

5.1. Excitation-induced contraction in spiral re-entry problem

A three-dimensional 100 mm × 100 mm × 12 mm slice of an orthotropic cardiac tissue is discretized into 21 × 21 × 2 eight-node coupled brick elements as described in [17, 19]. The myofibers are aligned in the horizontal x-direction and the sheets are aligned in the vertical y-direction. The material properties of the myocardium are summarized in Table 1. The initial value of the transmembrane potential is set to its resting value of Φ(X, t0) = −80 mV in the entire domain, except in the x = 0 plane where the transmembrane potential is set to Φ (X, t0) = −40 mV to trigger an excitation wave front. The initial value of the extracellular potential is set to Φe(X, t0) = 0 in the entire domain. The intra- and extracellular conductances are kept proportional to recover the special case of monodomain electromechanics. The nodes in the z = 0 plane are supported by uncoupled linear springs with stiffnesses of kx = ky = 10−3 N/mm and kz = 10−1 N/mm to provide a fairly unconstrained representation. Flux-free Neumann boundary conditions are applied to the parabolic and elliptic parts of the initial boundary value problem. We apply 150 fixed time increments of Δt = 3 ms followed by 1000 fixed time increments of Δt = 1 ms. Once the wave front has formed, it travels in the horizontal x–direction, depolarizing the tissue sample and causing contraction along the horizontal myocyte direction, see frame t = 75 ms in Figure 3. Then, the myocytes start to restitute their undeformed state as the block repolarizes, see frame t = 420 ms. To initiate a spiral wave, we follow the procedure suggested in [17]. We apply an additional current source of An external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms396431ig32.jpg + 5β[var phi]/βt between the interval t [set membership] [440, 460] ms in the rectangular region bounded by the coordinates x [set membership] [40,50] mm, y [set membership] [0, 55] mm and z [set membership] [0, 12]. The snapshots in Figure 3 for t > 460 ms demonstrate the initition, development and stable rotation of the reentrant spiral wave. In order to assess the computational efficiency of the proposed algorithm and its computational cost relative to the monodomain formulation in [18], we perform a successive mesh refinement with 11 ×11 ×2, 15 × 15 × 2 and 25 × 25 × 2 elements, respectively. Table 2 summarizes the norms of the global residual of the first five iterations for the 21 ×21 ×2 element discretization. Convergence is linear for the first few steps and then turns quadratic until the residual tolerance is reached. For the rest of the simulation the convergence remains quadratic throughout.

Figure 3
Initiation and rotation of a spiral re-entry in excitable and deformable cardiac tissue. The first, second and third row denote the transmembrane potential Φ, extracellular potential Φe and active fiber tension σ, respectively. ...
Table 1
Material parameters for spiral re-entry problem.
Table 2
Quadratic convergence of global residual for spiral re-entry problem.

Table 3 illustrates the total computation times until t = 1450 ms for both the monodomain and the bidomain model for the different discretizations. All simulations are carried out on a single 3.40 GHz Intel Pentium(R) processor on 32bit Windows XP operating system with 3.24GB RAM. The total computation time includes the writing of the whole mesh database and nodal data information of the state variables into.VTK type files for the visualization toolkit after each time step. As the mesh size decreases or the number of elements increases, the computation times of the bidomain formulation converge to those of the monodomain formulation. This result indicates that the bidomain formulation generates almost no extra cost, especially at large scale simulations. This is mainly due to the weak coupling between the elliptic and and parabolic parts of the bidomain equations in the case of proportional intra- and extracellular conductances.

Table 3
Computational times of mono- and bidomain-based spiral re-entry problem.

5.2. Excitation-contraction in generic biventricular heart

Next, we consider a three-dimensional generic biventricular heart, constructed from two truncated ellipsoids, see [52]. The dimensions and spatial discretization of the geometry are adopted from our previous contribution [18]. The heart is discretized with 13348 four-node coupled tetrahedral elements connected through 3059 nodes. The myofiber orientation is varied from +70° in the endocardium, the inner wall, to −70° in the epicardium, the outer wall, as described in [18]. The material parameters used for this model problem are given in Table 4. The displacement degrees of freedoms on the top base surface of the ventricles are fixed and flux-free boundary conditions are imposed on all surfaces. Excitation is initiated by assigning an elevated transmembrane potential of Φ0 = −10 mV at the nodes located at the upper part of septum. The remainder of the heart is assigned the resting potential of Φ(X, t0) = −80 mV. The initial value of the extracellular potential is set to Φe(X, t0) = 0 for the entire domain. An adaptive time stepping algorithm is used [57, 64], with a starting time increment value Δt(t0) = 5 ms bounded in Δt [set membership] [0.01, 10] ms.

Table 4
Material parameters for generic biventricular heart problem.

The deformed geometry of the generic heart and the color plots corresponding to the transmembrane potential, the extracellular potential and the active fiber tension in depolarization and repolarization phases are depicted in Figure 4 and Figure 5, respectively. The algorithm conceptually captures both the depolarization and the repolarization phases. Unfortunately, our relatively corse benchmark mesh cannot resolve the discontinuity at the apex accurately. The steep curvature at the apex induces a discontinuity in the conductances along the fiber direction, which artificially elongates the action potential duration and leads to an artificially slow reduction of the active fiber tension in the apex region. We have recently established a more accurate feature based Poisson interpolation, which resolves the apex discontinuity more accurately and avoids these spurious contractions [65]. The apex singularity will inherently disappear when we solve the underlying problem on real patient-specific geometries [64].

Figure 4
Coupled excitation induced contraction of generic heart model at various stages of depolarization. The first, second and third row denote the transmembrane potential Φ, extracellular potential Φe and active fiber tension σ, respectively. ...
Figure 5
Coupled excitation induced contraction of generic heart model at various stages of repolarization. The first, second and third row denote the transmembrane potential Φ, extracellular potential Φe and active fiber tension σ, respectively. ...

The convergence behaviour of the global Newton-Raphson iteration is depicted in Table 5 for various stages of depolarization and repolarization. The depolarization wave-front renders the problem stiff and smaller time increments are required. However, the repolarization phase of the problem is numerically easier to tackle because of the lack of steep potential gradients. Since the the intra- and extracellular conductances are not proportional, an equivalent monodomain conductance cannot be obtained for this particular problem. Assuming an effective monodomain conductance by means of geometric average of the conductances along the fiber and transverse directions, we have compared the computational costs of the bidomain and monodomain based electromechanics formulations. Since mono- and bidomain based simulations are not exactly equivalent though, we compare the computation time required for the depolarization wave to reach the bottom of the left ventricle. Table 6 shows the corresponding computation times for both simulations. The anisotropic distribution of the fiber directions and non-proportional conductances add additional complexity to the problem, which particularly affects the elliptic part of the bidomain equations. This is reflected by the additional 75% in computational cost in the excitation-contraction of generic heart problem.

Table 5
Quadratic convergence of global residual for generic biventricular heart problem.
Table 6
Computational times of mono- and bidomain-based generic biventricular heart problem.

6. Discussion and outlook

We have proposed a fully coupled finite element method for cardiac electromechanics, which is inherently modular, unconditionally stable, computationally efficient, and robust. Modularity is obtained by defining constitutive equations independently for source and flux terms of all three fields. Unconditional stability is ensured by using the backward Euler method for the temporal discretization of the local recovery variable r and the transmembrane potential Φ. Efficiency is guaranteed by the Newton Raphson iteration in combination with a consistent algorithmic linearization of the underlying weak forms of the governing equations. Robustness is enhanced by the fully coupled, monolithic solution strategy.

In weakly coupled multifield problems, decoupling based on operator-splitting methods can improve computational efficiency by dividing a big non-symmetric problem into smaller symmetric subproblems. Depending on the degree of coupling, operator splitting might induce instabilities, which can be overcome by refining the time step size. In general, there is an inherent competition between saving computer time by solving smaller symmetric subsystems and losing computer time by smaller time step sizes relative to monolithic schemes.

For the particular problem of cardiac electromechanics, steep potential gradients and strong coupling between the electrical and mechanical fields could render operator splitting inefficient. As we have demonstrated here, the use of the classical Newton Raphson solution strategy, along with consistently linearized residuals, ensures the quadratic convergence of the iterative solver. Depending on the number of overall degrees of freedom, an consistently linearized implicit algorithm might be superior to decoupled solution strategies based on operator-splitting schemes. Specifically, in our setup, the proposed implicit approach allows us to use straightforward time-adaptive schemes, in which the time step size is adjusted based on the number of iteration steps. In addition, the proposed approach is highly modular and generic. It is independent from the particular choice of electrical and mechanical constitutive models, and can be implemented into basically any commercial nonlinear finite element environment.

We believe that, despite intense research, the use of finite element methods in cardiac electromechanics is still under-appreciated. With this manuscript, we hope to encourage future research along this direction. Discretizing all relevant fields with unique discretization schemes both in time and space, like we have proposed here, enables the use of generic finite element packages and commercially available solvers. These typically come with well-established tools like error estimators or readily available temporally and spatially adaptive schemes. Overall, robust and efficient solution algorithms in electromechanics have the potential to enable the design and refinement of treatment strategies for heart disease, the number one cause of death in the industrialized nations.

  • Cardiac contraction is initiated by smoothly propagating excitation waves.
  • Existing models solve this excitation-contraction coupling in a staggered sense.
  • Here we propose a monolithic, fully implicit excitation-contraction algorithm.
  • We represent excitation using the bidomain model of electrophysiology.
  • Our model accurately predicts excitation and contraction in an idealized heart.

Acknowledgments

SG has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant PCIG09-GA-2011-294161. EK has received funding from the National Science Foundation CAREER award CMMI-0952021 and from the National Institutes of Health Grant U54 GM072970. MK has received funding from the Deutsche Forschungsgemein-schaft (DFG) under grant Ka-1163/18-1.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Abilez OJ, Wong J, Prakash R, Deisseroth K, Zarins CK, Kuhl E. Multiscale computational models for optogenetic control of cardiac function. Biophysical Journal. 2011;101:1326–1334. [PubMed]
2. Aliev RR, Panfilov AV. A simple two-variable model of cardiac excitation. Chaos, Solitons and Fractals. 1996;7:293–301.
3. Bacharova L, Mateasik A, Krause R, Prinzen FW, Auricchio A, Potse M. The effect of reduced intracellular coupling on electrocardiographic signs of left ventricular hypertrophy. Journal of Electrocardiology. 2011;44:571–576. [PubMed]
4. Borg TK, CJB The collagen matrix of the heart. Federation Proceedings. 1981;40:2037–2041. [PubMed]
5. Bothe W, Kuhl E, Kvitting JP, Rausch MK, Göktepe S, Swanson JC, Farahmandnia S, Ingels NB, Miller DC. Rigid, complete annuloplasty rings increase anterior mitral leaflet strains in the normal beating ovine heart. Circulation. 2011;124:S81–S96. [PMC free article] [PubMed]
6. Brera M, Jerome JW, Mori Y, Sacco R. A conservative and monotone mixed-hybridized finite element approximation of transport problems in heterogeneous domains. Computer Methods in Applied Mechanics and Engineering. 2010;199:2709–2720.
7. Chen MQ, Wong J, Kuhl E, Giovangrandi LB, Kovacs GTA. Characterization of electrophysiological conduction in cardiomyocyte co-cultures using co-occurrence analysis. Computer Methods in Biomechanics and Biomedical Engineering. 2012 doi: 10.1080/10255842.2011.615310. [PubMed] [Cross Ref]
8. Costa KD, Holmes JW, McCulloch AD. Modelling cardiac mechanical properties in three dimensions. Philosophical Transactions of the Royal Society London A. 2001;359:1233–1250.
9. Dal H, Göktepe S, Kaliske M, Kuhl E. A fully implicit finite element method for bidomain models of cardiac electrophysiology. Computer Methods in Biomechanics and Biomedical Engineering. 2011 doi: 10.1080/10255842.2011.554410. [PubMed] [Cross Ref]
10. Dal H, Göktepe S, Kaliske M, Kuhl E. A three-field, bi-domain based approach to the strongly coupled electromechanics of the heart. PAMM. 2011;11:931–934.
11. Ennis DB, Nguyen TC, Riboh JC, Wigström L, Harrington KB, Daughters GT, Ingels NB, Miller DC. Myofiber angle distributions in the ovine left ventricle do not conform to computationally optimized predictions. Journal of Biomechanics. 2008;41:3219–3224. [PMC free article] [PubMed]
12. Fitzhugh R. Impulses and physiological states in theoretical models of nerve induction. Biophysical Journal. 1961;1:455–466. [PubMed]
13. Flory PJ. Thermodynamic relations for high elastic materials. Transactions of the Faraday Society. 1961;57:829–838.
14. Franzone PC, Pavarino LF. A parallel solver for reaction-diffusion systems in computational electrophysiology. Mathematical Models and Methods in Applied Sciences. 2004;14:883–911.
15. Gerardo-Giorda L, Mirabella L, Nobile F, Perego M, Veneziani A. A model-based block-triangular preconditioner for the bidomain system in electrocardiology. Journal of Computational Physics. 2009;228:3625–3639.
16. Göktepe S, Acharya SNS, Wong J, Kuhl E. Computational modeling of passive myocardium. International Journal for Numerical Methods in Biomedical Engineering. 2011;27:1–12.
17. Göktepe S, Kuhl E. Computational modeling of cardiac electrophysiology: A novel finite element approach. International Journal for Numerical Methods in Engineering. 2009;79:156–178.
18. Göktepe S, Kuhl E. Electromechanics of the heart: a unified approach to the strongly coupled excitation–contraction problem. Computational Mechanics. 2010;45:227–243.
19. Göktepe S, Wong J, Kuhl E. Atrial and ventricular fibrillation - computational simulation of spiral waves in cardiac tissue. Archive of Applied Mechanics. 2010;80:569–580.
20. Guccione JM, McCulloch AD. Mechanics of active contraction in cardiac muscle: Part I - constitutive relations for fiber stress that describe deactivation. Journal of Biomechanical Engineering. 1993;115:72–81. [PubMed]
21. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conductance and excitation in nerve. Journal of Physiology. 1952;117:500–544. [PubMed]
22. Holzapfel GA, Ogden RW. Constitutive modelling of passive myocardium: A structurally based framework for material characterization. Philosophical Transactions Series A, Mathematical, Physical, and Engineering Sciences. 2009;367:3445–3475. [PubMed]
23. Hunter PJ, Lie WW, McCulloch AD, Noble D. Multiscale modeling: Physiome project standards, tools, and databases. IEEE Computer. 2006;39:48–54.
24. Hunter PJ, Pullan AJ, Smaill BH. Modeling total heart function. Annular Review of Biomedical Engineering. 2003;5:147–177. [PubMed]
25. Itoh A, Krishnamurthy G, Swanson J, Ennis D, Bothe W, Kuhl E, Karlsson M, Davis L, Miller DC, Ingels NB. Active stiffening of mitral valve leaflets in the beating heart. American Journal of Physiology Heart and Circulation Physiology. 2009;296:1766–1773. [PubMed]
26. Keener JP, Bogar K. A numerical method for the solution of the bidomain equations in cardiac tissue. Chaos: An Interdisciplinary Journal of Nonlinear Science. 1998;8:234–241. [PubMed]
27. Kerckhoffs R, Healy S, Usyk T, McCulloch A. Computational methods for cardiac electromechanics. Proceedings of the IEEE. 2006;94:769–783.
28. Kohl P, Hunter P, Noble D. Stretch-induced changes in heart rate and rhythm: Clinical observations, experiments and mathematical models. Progress in Biophysics and Molecular Biology. 1999;71:91–138. [PubMed]
29. Kotikanyadanam M, Göktepe S, Kuhl E. Computational modeling of electrocardiograms - a finite element approach towards cardiac excitation. International Journal for Numerical Methods in Biomedical Engineering. 2010;26:524–533.
30. Krishnamurthy G, Ennis DB, Itoh A, Bothe W, Swanson-Birchill JC, Karlsson M, Kuhl E, Miller DC, Ingels NB. Material properties of the ovine mitral valve anterior leaflet in vivo from inverse finite element analysis. American Journal of Physiology Heart and Circulation Physiology. 2008;295:H1141–H1149. [PubMed]
31. Marsden JE, Hughes TJR. Mathematical Foundations of Elasticity. Prentice-Hall; Englewood Cliffs, New Jersey: 1983.
32. Miehe C. Aspects of the formulation and finite element implementation of large strain isotropic elasticity. International Journal for Numerical Methods in Engineering. 1994;37:1981–2004.
33. Miller WT, Geselowitz DB. Simulation studies of the electrocardiogram i: the normal heart. Circulation Research. 1978;43:301–315. [PubMed]
34. Mirabella L, Nobile F, Veneziani A. An a posteriori error estimator for model adaptivity in electrocardiology. Computer Methods in Applied Mechanics and Engineering. 2011;200:2727–2737.
35. Munteanu M, Pavarino LF, Scacchi S. A scalable newton–krylov–schwarz method for the bidomain reaction-diffusion system. SIAM Journal on Scientific Computing. 2009;31:3861–3883.
36. Nagumo J, Arimoto S, Yoshizawa S. Active pulse transmission line simulating nerve axon. Proceedings of the Instiute of Radio Engineers. 1962;50:2061–2070.
37. Nash MP, Hunter PJ. Computational mechanics of the heart. from tissue structure to ventricular function. Journal of Elasticity. 2000;61:113–141.
38. Nash MP, Panfilov AV. Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biology. 2004;85:501–522. [PubMed]
39. Nickerson D, Nash M, Nielsen P, Smith N, Hunter P. Computational multiscale modeling in the IUPS physiome project: Modeling cardiac electromechanics. Systems Biology. 2006;50:617–630.
40. Pennacchio M, Simoncini V. Efficient algebraic solution of reaction–diffusion systems for the cardiac excitation process. Journal of Computational and Applied Mathematics. 2002;145:49–70.
41. Pennacchio M, Simoncini V. Algebraic multigrid preconditioners for the bidomain reaction–diffusion system. Applied Numerical Mathematics. 2009;59:3033–3050.
42. Plank G, Burton RAB, Hales P, Bishop M, Mansoori T, Bernabeau MO, Garny A, Prassl AJ, Bollensdorff C, Mason F, Mahmood F, Rodriguez B, Grau V, Schneider JE, Gavaghan G, Kohl P. Generation of histo-anatomically representative models of the individual heart: tools and application. Philosophical Transactions of the Royal Society London A. 2008;197:4051–4061. [PMC free article] [PubMed]
43. Plank G, Liebmann M, dos Santos RW, Vigmond EJ, Haase G. IEEE Transactions on Biomedical Engineering. 2007;54:585–596. [PubMed]
44. Potse M, Dubé B, Richer J, Vinet A, Gulrajani RM. A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart. IEEE Transactions on Biomedical Engineering. 2006;53:2425–2435. [PubMed]
45. Prassl A, Kickinger F, Ahammer H, Grau V, Schneider J, Hofer E, Vigmond E, Trayanova N, Plank G. Automatically generated, anatomically accurate meshes for cardiac electrophysiology problems. Biomedical Engineering, IEEE Transactions on. 2009;56:1318–1330. [PMC free article] [PubMed]
46. Rohmer D, Sitek A, Gullberg GT. Reconstruction and visualization of fiber and laminar structure in the normal human heart from ex vivo diffusion tensor magnetic resonance imaging (DTMRI) data. Investigative Radiology. 2007;42:777–789. [PubMed]
47. Sainte-Marie J, Chapelle D, Cimrman R, Sorine M. Modeling and estimation of the cardiac electromechanical activity. Computers & Structures. 2006;84:1743–1759.
48. dos Santos R, Plank G, Bauer S, Vigmond E. Parallel multigrid preconditioner for the cardiac bidomain model, Biomedical Engineering. IEEE Transactions on. 2004;51:1960–1968. [PubMed]
49. Scacchi S. A hybrid multilevel schwarz method for the bidomain model. Computer Methods in Applied Mechanics and Engineering. 2008;197:4051–4061.
50. Scacchi S. A multilevel hybrid newton-krylov-schwarz method for the biodomain model of electrocardiology. Computer Methods in Applied Mechanics and Engineering. 2011;200:717–725.
51. Scacchi S, Franzone PC, Pavarino L, Taccardi B. A reliability analysis of cardiac repolarization time markers. Mathematical Biosciences. 2009;219:113–128. [PubMed]
52. Sermesant M, Rhode K, Sanchez-Ortiz G, Camara O, Andriantsimiavona R, Hegde S, Rueckert D, Lambiase P, Bucknall C, Rosenthal E, Delingette H, Hill D, Ayache N, Razavi R. Simulation of cardiac pathologies using an electromechanical biventricular model and XMR interventional imaging. Medical Image Analysis. 2005;9:467–480. [PubMed]
53. Simo JC, Taylor RL. Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering. 1991;85:273–310.
54. Skouibine K, Trayanova N, Moore P. A numerically efficient model for simulation of defibrillation in an active bidomain sheet of myocardium. Mathematical Biosciences. 2000;166:85–100. [PubMed]
55. Spencer AJM. Theory of invariants. In: Eringen A, editor. Continuum Physics. Vol. 1. Academic Press; New York: 1971.
56. Sundnes J, Lines GT, Tveito A. An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso. Mathematical Biosciences. 2005;194:233–248. [PubMed]
57. Taylor RL. User Manual. University of California; Berkeley: 2011. FEAP - A Finite Element Analysis Program, Version 8.3.
58. Tsamis A, Bothe W, Kvitting JP, Swanson JC, Miller DC, Kuhl E. Active contraction of cardiac muscle: In vivo characterization of mechanical activation sequences in the beating heart. Journal of the Mechanical Behavior of Biomedical Materials. 2011;4:1167–1176. [PMC free article] [PubMed]
59. Tung L. PhD thesis. MIT; 1978. A Bidomain model for describing ischaemic myocardial DC potentials.
60. Usyk TP, LeGrice IJ, McCulloch AD. Computational model of three-dimensional cardiac electromechanics. Computing and Visualization in Science. 2002;4:249–257.
61. Vigmond E, Aguel F, Trayanova N. Computational techniques for solving the bidomain equations in three dimensions, Biomedical Engineering. IEEE Transactions on. 2002;49:1260–1269. [PubMed]
62. Vigmond E, dos Santos RW, Prassl A, Deo M, Plank G. Solvers for the cardiac bidomain equations. Progress in Biophysics and Molecular Biology. 2008;96:3–18. [PMC free article] [PubMed]
63. Vigmond EJ, Hughes M, Plank G, Leon LJ. Computational tools for modeling electrical activity in cardiac tissue. Journal of Electrocardiology. 2003;36:69–74. [PubMed]
64. Wong J, Göktepe S, Kuhl E. Computational modeling of electrochemical coupling: A novel finite element approach towards ionic models for cardiac electrophysiology. Computer Methods in Applied Mechanics and Engineering. 2011;200:3139–3158.
65. Wong J, Kuhl E. Generating fiber orientation maps in human heart models using Poisson interpolation. 2012. submitted for publication. [PubMed]
66. Zhang Y, Bazilevs Y, Goswami S, Bajaj CL, Hughes TJR. Patient-specific vascular nurbs modeling for isogeometric analysis of blood flow. Computer Methods in Applied Mechanics and Engineering. 2007;196:2943–2959. [PMC free article] [PubMed]