Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3501134

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Governing equations of the cardiac electromechanics
- 3. The finite element formulation
- 4. Constitutive equations
- 5. Representative numerical examples
- 6. Discussion and outlook
- References

Authors

Related links

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.004PMCID: PMC3501134

NIHMSID: NIHMS396431

See other articles in PMC that cite the published article.

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.

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 . In the second step, the parabolic part is solved for a constant external potential * _{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

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.

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

A *body B* is a three-dimensional manifold consisting of material points
*B*. The motion of the body is defined by a one-parameter function of time via bijective mappings

$$\mathit{\chi}(\mathcal{P},t)=\{\begin{array}{l}B\to \mathcal{B}(\mathcal{P},t)\in {\mathbb{R}}^{3}\times {\mathbb{R}}_{+}\hfill \\ \mathcal{P}\mapsto \mathit{x}={\mathit{\chi}}_{t}(\mathcal{P})=\mathit{\chi}(\mathcal{P},t).\hfill \end{array}$$

(1)

The point ** x** =
(
,

$${\mathit{\varphi}}_{t}(\mathit{X})=\{\begin{array}{l}{\mathcal{B}}_{0}\to \mathcal{B}\in {\mathbb{R}}^{3}\hfill \\ \mathit{X}\mapsto \mathit{x}=\mathit{\varphi}(\mathit{X},t)\hfill \end{array}$$

(2)

maps the reference configuration ** X**
of a material point onto the spatial counterpart

$$\mathit{F}:{T}_{X}{\mathcal{B}}_{0}\to {T}_{x}\mathcal{B};\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathit{F}:={\nabla}_{X}{\mathit{\varphi}}_{t}(\mathit{X})$$

(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 * _{X}* [•] and

$$\mathit{F}={\mathit{F}}_{\mathit{vol}}\overline{\mathit{F}}.$$

(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 *f*_{0} and *s*_{0} such that |*f*_{0}|**_{G}** = 1 and |

$$\mathit{f}=\overline{\mathit{F}}{\mathit{f}}_{0}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathit{s}=\overline{\mathit{F}}{\mathit{s}}_{0},$$

(5)

respectively.

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

$$\text{State}\phantom{\rule{0.16667em}{0ex}}:=\{\mathit{\varphi}(\mathit{X},t),\mathrm{\Phi}(\mathit{X},t),{\mathrm{\Phi}}_{e}(\mathit{X},t))\},$$

(6)

namely the deformation map *ϕ _{t}*, the extracellular potential Φ

The electrophysiology of a material point **X** of cardiac tissue at time *t* is defined as State(**X**, *t*):= {Φ* _{i}*(

$${\int}_{\partial \mathcal{S}}{\mathfrak{q}}_{t}\xb7\mathfrak{n}\phantom{\rule{0.16667em}{0ex}}dA=0\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}};\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathfrak{q}}_{t}={\mathfrak{q}}_{i}+{\mathfrak{q}}_{e},$$

(7)

where , and 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)

$$\text{div}\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{t}=0.$$

(8)

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

$${\mathfrak{q}}_{i}=-{\widehat{\mathfrak{D}}}_{i}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{i}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}};\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathfrak{q}}_{e}=-{\widehat{\mathfrak{D}}}_{e}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{e}.$$

(9)

Here,
and
are the intra- and extracellular conductance tensors, respectively. Moreover, we define the transmembrane current *i _{T}* as the sum of the ionic and capacitive currents from the intracellular space into the extracellular space such that

$${i}_{T}=\zeta \left({\mathcal{C}}_{m}\stackrel{.}{\mathrm{\Phi}}+{\mathcal{J}}_{\mathit{ion}}(\mathrm{\Phi},r)\right)\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathrm{\Phi}={\mathrm{\Phi}}_{i}-{\mathrm{\Phi}}_{e}.$$

(10)

Here, *ζ* denotes the membrane surface to volume ratio in a periodic representative unit of cardiac tissue.
and
are the membrane capacitance per unit area and the current density flowing through the membrane ion channels, see Figure 1. The notation [•]:= D[•]/D*t* is employed for the material time derivative throughout the manuscript. The transmembrane current *i _{T}* 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

One-dimensional representaion of myocardium and approximation of the cell membrane by resistor-capacitor circuit. The circuit consists of a nonlinear resistor
dependent on the ionic concentrations (sodium [Na]^{+}, potassium [K]^{+} and leakage currents) **...**

$$\text{div}\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{i}=-{i}_{T}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}};\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{div}\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{e}={i}_{T}.$$

(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
and the extracardiac domain
, such that the total domain of interest is
=
. The boundary of the cardiac domain,
=
, consists of the endocardium
, the boundary of the inner walls, and the epicardium
, the boundary of the outer walls, such that
∩
= , see Figure 2. Similarly, the boundary of the extracardiac domain,
=
, consists of the epicardium
and the outer boundary
, such that
∩
= . We can further decompose each boundary into Dirichlet and Neumann-type boundaries, such that
$\partial {\mathcal{B}}^{a}=\partial {\mathcal{B}}_{\phi}^{a}\cup \partial {\mathcal{B}}_{q}^{a}$ and
$\partial {\mathcal{B}}_{\phi}^{a}\cap \partial {\mathcal{B}}_{q}^{a}=\varnothing $, where *a*:= {*c*, *e*}. In the absence of an extracardiac domain, the epicardium is identical with the boundary of the extracardiac domain,
=
. We can then rewrite the Neumann boundary conditions at the interface
as follows

Illustration of the cardiac and extracardiac domains
and
and the endocardium
, the epicardium
, and outer boundary
.

$${\mathfrak{q}}_{t}\xb7\mathfrak{n}={\mathfrak{q}}_{0}\xb7\mathfrak{n}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{on}\phantom{\rule{0.38889em}{0ex}}\partial {\mathcal{B}}_{q}^{s}.$$

(12)

The extracardiac medium is represented by a single potential term Φ_{0} and the corresponding flux vector
= −
Φ_{0},
denoting the extracardiac conductance. It possesses no transmembrane currents, and the set of transient diffusion equation takes the simple form

$$\text{div}\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{0}=0\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{in}\phantom{\rule{0.38889em}{0ex}}{\mathcal{B}}^{e}.$$

(13)

Equation (12) states that the total current can be discretized continuously in the intra- and extracardiac medium leading to the following boundary conditions:

$${\mathrm{\Phi}}_{0}={\overline{\mathrm{\Phi}}}_{0}\phantom{\rule{0.38889em}{0ex}}\text{on}\phantom{\rule{0.38889em}{0ex}}\partial {\mathcal{B}}_{\phi}^{e}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}{\mathfrak{q}}_{0}\xb7\mathfrak{n}={i}_{\mathit{app}}\phantom{\rule{0.38889em}{0ex}}\text{on}\phantom{\rule{0.38889em}{0ex}}\partial {\mathcal{B}}_{q}^{e},$$

(14)

where *i _{app}* 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

$${\mathrm{\Phi}}_{i,e}={\overline{\mathrm{\Phi}}}_{i,e}\phantom{\rule{0.38889em}{0ex}}\text{on}\phantom{\rule{0.38889em}{0ex}}\partial {\mathcal{B}}_{\phi}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}{\mathfrak{q}}_{i,e}\xb7\mathfrak{n}={\overline{\mathfrak{q}}}_{i,e}\phantom{\rule{0.38889em}{0ex}}\text{on}\phantom{\rule{0.38889em}{0ex}}\partial {\mathcal{B}}_{q},$$

(15)

such that
=
and
∩
= . The bidomain equations (11)_{1} and (8) are recast into a more coherent form

$$\stackrel{.}{\mathrm{\Phi}}=\text{div}({\mathfrak{D}}_{i}\xb7\nabla \mathrm{\Phi})+\text{div}({\mathfrak{D}}_{i}\xb7\nabla {\mathrm{\Phi}}_{e})+{f}^{\mathrm{\Phi}}(\mathrm{\Phi},r),$$

(16)

$$0=-\text{div}\phantom{\rule{0.16667em}{0ex}}({\mathfrak{D}}_{i}\xb7\nabla \mathrm{\Phi})-\text{div}\phantom{\rule{0.16667em}{0ex}}(\mathfrak{D}\xb7\nabla {\mathrm{\Phi}}_{e})$$

(17)

in along with the normalized conductances

$$\mathfrak{D}\phantom{\rule{0.16667em}{0ex}}:=\frac{1}{\zeta {\mathcal{C}}_{m}}({\widehat{\mathfrak{D}}}_{i}+{\widehat{\mathfrak{D}}}_{e})\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}{\mathfrak{D}}_{i}\phantom{\rule{0.16667em}{0ex}}:=\frac{1}{\zeta {\mathcal{C}}_{m}}{\widehat{\mathfrak{D}}}_{i}.$$

(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 = * _{x}* =

$${\mathfrak{q}}_{a}\phantom{\rule{0.16667em}{0ex}}:={\mathfrak{D}}_{a}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{a}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{where}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}a=\{i,e\}.$$

(19)

Moreover, the Equations (16, 17) are assumed to be valid in the spatial setting, e.g. divergence operator is taken to be div(•) = * _{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

$$\stackrel{.}{\mathrm{\Phi}}=J\phantom{\rule{0.16667em}{0ex}}\text{div}({J}^{-1}{\mathcal{D}}_{i}\xb7{\nabla}_{x}\mathrm{\Phi})+{J}^{-1}\phantom{\rule{0.16667em}{0ex}}\text{div}({J}^{-1}{\mathfrak{D}}_{i}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{e})+{\mathcal{F}}^{\phi},$$

(20)

$$0=J\phantom{\rule{0.16667em}{0ex}}\text{div}\phantom{\rule{0.16667em}{0ex}}\left({J}^{-1}{\mathfrak{D}}_{i}\xb7{\nabla}_{x}\mathrm{\Phi}\right)+J\phantom{\rule{0.16667em}{0ex}}\text{div}\phantom{\rule{0.16667em}{0ex}}\left({J}^{-1}\mathfrak{D}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{e}\right)$$

(21)

for a unit undeformed volume in . Equation (20) is a parabolic partial differential equation and equation (21) defines an elliptic partial differential equation. is the Fitzhugh-Nagumo type generic function defining the excitation induced ionic currents through the cell membrane from the intracellular to the extracellular domain.

The balance of linear momentum

$$J\phantom{\rule{0.16667em}{0ex}}\text{div}[{J}^{-1}\widehat{\mathit{\tau}}]+\mathit{b}=\mathbf{0}\phantom{\rule{0.38889em}{0ex}}\text{in}\phantom{\rule{0.38889em}{0ex}}\mathcal{B}$$

(22)

states the quasi-static mechanical equilibrium in terms of the Kirchhoff stress tensor ** and a given body force **** b**. The balance of linear momentum (22) is subjected to essential boundary conditions

The solution of the coupled initial boundary value problem requires the definition of the constitutive relations which characterize the Kirchhoff stress tensor **, the potential flux vectors and the current source
appearing in Equations (20–22). The Kirchhoff stress tensor **** can be decomposed into the passive **^{pas} and active ^{act} parts

$$\widehat{\mathit{\tau}}={\widehat{\mathit{\tau}}}^{\text{pas}}(\mathit{g};\mathit{F})+{\widehat{\mathit{\tau}}}^{\text{act}}(\mathit{g};\mathit{F},\mathrm{\Phi}),$$

(23)

see [38, 47]. The passive part ^{pas} is solely governed by mechanical deformation, while the active part ^{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 Φ

$${\mathfrak{q}}_{i}={\mathfrak{D}}_{i}\xb7{\nabla}_{x}(\mathrm{\Phi})\phantom{\rule{0.16667em}{0ex}};\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{ie}={\mathfrak{D}}_{i}\xb7{\nabla}_{x}({\mathrm{\Phi}}_{e})\phantom{\rule{0.16667em}{0ex}};\phantom{\rule{0.16667em}{0ex}}{\mathfrak{q}}_{e}=\mathfrak{D}\xb7{\nabla}_{x}({\mathrm{\Phi}}_{e})$$

(24)

through the deformation-dependent anisotropic spatial conduction tensors
(** g; F**,

$${\widehat{\mathcal{F}}}^{\phi}={\widehat{\mathcal{F}}}_{\text{e}}^{\phi}(\mathrm{\Phi},r)+{\widehat{\mathcal{F}}}_{\text{m}}^{\phi}(\mathit{g};\mathit{F},{\mathit{f}}_{0},\mathrm{\Phi}).$$

(25)

It is additively decomposed into the excitation-induced purely electrical part
${\widehat{\mathcal{F}}}_{\text{e}}^{\phi}(\mathrm{\Phi},r)$ and the stretch-induced mechano-electrical part
${\widehat{\mathcal{F}}}_{\text{m}}^{\phi}(\mathit{g};\mathit{F},{\mathit{f}}_{0},\mathrm{\Phi})$. 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
${\widehat{\mathcal{F}}}_{\text{e}}^{\phi}$ 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

$$\stackrel{.}{r}={\widehat{f}}^{r}(\mathrm{\Phi},r).$$

(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 **, **** and
${\widehat{\mathcal{F}}}_{m}^{\phi}$. 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.

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 (20–22) and consistently linearize them along the three independent field variables, namely the deformation map ** ϕ** (

A conventional Galerkin procedure is applied to construct the weak forms of the governing field equations (20–22). For this purpose, the residual equations are multiplied by the weight functions *δ*** ϕ**,

$${G}^{\varphi}={G}_{int}^{\varphi}(\delta \mathit{\varphi},\mathit{\varphi},\mathrm{\Phi})-{G}_{\text{ext}}^{\varphi}(\delta \mathit{\varphi})=0,$$

(27)

the parabolic part of the bidomain equations

$${G}^{\phi}={G}_{int}^{\phi}(\delta \mathrm{\Phi},\mathit{\varphi},\mathrm{\Phi},{\mathrm{\Phi}}_{e})-{G}_{\text{ext}}^{\phi}(\delta \mathrm{\Phi},\mathit{\varphi},\mathrm{\Phi},{\mathrm{\Phi}}_{e})=0,$$

(28)

and the elliptic part of the bidomain equations

$${G}^{{\phi}_{e}}={G}_{int}^{{\phi}_{e}}(\delta {\mathrm{\Phi}}_{e},\mathit{\varphi},\mathrm{\Phi},{\mathrm{\Phi}}_{e})-{G}_{\text{ext}}^{{\phi}_{e}}(\delta {\mathrm{\Phi}}_{e},\mathit{\varphi},\mathrm{\Phi},{\mathrm{\Phi}}_{e})=0.$$

(29)

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

$$\begin{array}{l}{G}_{int}^{\varphi}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathit{\varphi})\phantom{\rule{0.16667em}{0ex}}:\phantom{\rule{0.16667em}{0ex}}\widehat{\mathit{\tau}}\phantom{\rule{0.16667em}{0ex}}dV,\\ {G}_{\text{ext}}^{\varphi}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\mathcal{B}}\delta \mathit{\varphi}\xb7\mathit{b}\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\partial {\mathcal{S}}_{t}}\delta \mathit{\varphi}\xb7\overline{\mathit{t}}\phantom{\rule{0.16667em}{0ex}}dA,\end{array}$$

(30)

where the body force ** b** and the surface traction

$$\begin{array}{l}{G}_{int}^{\phi}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\mathcal{B}}\left(\delta \mathrm{\Phi}\phantom{\rule{0.16667em}{0ex}}\stackrel{.}{\mathrm{\Phi}}+{\nabla}_{x}(\delta \mathrm{\Phi})\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{ie})\right)dV,\\ {G}_{\text{ext}}^{\phi}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\mathcal{B}}\delta \mathrm{\Phi}{\widehat{\mathcal{F}}}^{\phi}dV+{\int}_{\partial {\mathcal{S}}_{q}}\delta \mathrm{\Phi}\phantom{\rule{0.16667em}{0ex}}\overline{\mathfrak{q}}\phantom{\rule{0.16667em}{0ex}}dA\end{array}$$

(31)

along with the definitions

$${\widehat{\mathfrak{q}}}_{i}={\mathfrak{D}}_{i}\xb7{\nabla}_{x}\mathrm{\Phi},\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathfrak{q}}}_{ie}={\mathfrak{D}}_{i}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{e},\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathfrak{q}}}_{e}=\mathfrak{D}\xb7{\nabla}_{x}{\mathrm{\Phi}}_{e}.$$

(32)

Likewise, the residual expressions ${G}_{int}^{{\phi}_{e}}$ and ${G}_{\text{ext}}^{{\phi}_{e}}$ in (29) of the elliptic part of the bidomain equations (21) are given by

$$\begin{array}{l}{G}_{\text{int}}^{{\phi}_{e}}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\mathcal{B}}{\nabla}_{x}(\delta {\mathrm{\Phi}}_{e})\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{e})\phantom{\rule{0.16667em}{0ex}}dV,\\ {G}_{\text{ext}}^{{\phi}_{e}}\phantom{\rule{0.16667em}{0ex}}:={\int}_{\partial \mathcal{B}}\delta {\mathrm{\Phi}}_{e}\phantom{\rule{0.16667em}{0ex}}{\overline{\mathfrak{q}}}_{t}\phantom{\rule{0.16667em}{0ex}}dA.\end{array}$$

(33)

The surface flux terms
$\overline{\mathfrak{q}}=({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{ie})\xb7\mathfrak{n}$ and
${\overline{\mathfrak{q}}}_{t}=({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{e})\xb7\mathfrak{n}$ are the prescribed natural boundary conditions. In contrast to the mechanical external term in (30)_{2},
${G}_{\text{ext}}^{\phi}$ depends explicitly upon the field variables through the non-linear source term
.

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 ** ϕ** =

$$\begin{array}{l}\text{Lin}\phantom{\rule{0.16667em}{0ex}}{G}^{\varphi}{\mid}_{\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}}}\phantom{\rule{0.16667em}{0ex}}:=\phantom{\rule{0.16667em}{0ex}}{G}^{\varphi}(\delta \mathit{\varphi},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}})+\mathrm{\Delta}{G}^{\varphi}(\delta \mathit{\varphi},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}};\mathrm{\Delta}\mathit{\varphi},\mathrm{\Delta}\mathrm{\Phi})=0,\\ \text{Lin}\phantom{\rule{0.16667em}{0ex}}{G}^{\phi}{\mid}_{\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\stackrel{\sim}{\mathrm{\Phi}}}_{e}}\phantom{\rule{0.16667em}{0ex}}:=\phantom{\rule{0.16667em}{0ex}}{G}^{\phi}(\delta \mathrm{\Phi},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\stackrel{\sim}{\mathrm{\Phi}}}_{e})+\mathrm{\Delta}{G}^{\phi}(\delta \mathrm{\Phi},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\stackrel{\sim}{\mathrm{\Phi}}}_{e};\mathrm{\Delta}\mathit{\varphi},\mathrm{\Delta}\mathrm{\Phi},\mathrm{\Delta}{\mathrm{\Phi}}_{e})=0,\\ \text{Lin}\phantom{\rule{0.16667em}{0ex}}{G}^{{\phi}_{e}}{\mid}_{\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\mathrm{\Phi}}_{e}}\phantom{\rule{0.16667em}{0ex}}:=\phantom{\rule{0.16667em}{0ex}}{G}^{{\phi}_{e}}(\delta {\mathrm{\Phi}}_{e},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\stackrel{\sim}{\mathrm{\Phi}}}_{e})+\mathrm{\Delta}{G}^{{\phi}_{e}}(\delta {\mathrm{\Phi}}_{e},\stackrel{\sim}{\mathit{\varphi}},\stackrel{\sim}{\mathrm{\Phi}},{\stackrel{\sim}{\mathrm{\Phi}}}_{e};\mathrm{\Delta}\mathit{\varphi},\mathrm{\Delta}\mathrm{\Phi},\mathrm{\Delta}{\mathrm{\Phi}}_{e})=0.\end{array}$$

(34)

The incremental terms Δ*G ^{ϕ}* and Δ

$$\mathrm{\Delta}{G}^{\gamma}=\mathrm{\Delta}{G}_{int}^{\gamma}-\mathrm{\Delta}{G}_{\text{ext}}^{\gamma}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\gamma =\varphi ,\phi ,{\phi}_{e},$$

(35)

based on the definitions in (20–22). We then start with the elaboration of the increment
$\mathrm{\Delta}{G}_{int}^{\varphi}$ according to (30)_{1}

$$\mathrm{\Delta}{G}_{int}^{\varphi}={\int}_{\mathcal{B}}\mathrm{\Delta}({\nabla}_{x}(\delta \mathit{\varphi}))\phantom{\rule{0.16667em}{0ex}}:\phantom{\rule{0.16667em}{0ex}}\widehat{\mathit{\tau}}+{\nabla}_{x}(\delta \mathit{\varphi})\phantom{\rule{0.16667em}{0ex}}:\phantom{\rule{0.16667em}{0ex}}\mathrm{\Delta}\widehat{\mathit{\tau}}\phantom{\rule{0.16667em}{0ex}}dV.$$

(36)

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

$$\begin{array}{l}\mathrm{\Delta}({\nabla}_{x}(\delta \mathit{\varphi}))=-{\nabla}_{x}(\delta \mathit{\varphi}){\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi});\\ \mathrm{\Delta}\widehat{\mathit{\tau}}={\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}\widehat{\mathit{\tau}}+{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\widehat{\mathit{\tau}}+\widehat{\mathit{\tau}}{\nabla}_{x}^{T}(\mathrm{\Delta}\mathit{\varphi})+{\mathfrak{C}}^{\varphi \phi}\mathrm{\Delta}\mathrm{\Phi},\end{array}$$

(37)

where **£**_{Δ}_{ϕ}** denotes the objective Lie derivative along the increment Δ****_{ϕ}** and can be expressed as

$$\begin{array}{l}{\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}\widehat{\mathit{\tau}}={\mathbb{C}}^{\varphi \varphi}:{\scriptstyle \frac{1}{2}}{\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}\mathit{g}={\mathbb{C}}^{\varphi \varphi}:(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}));\\ {\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}\mathit{g}=\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})+{\nabla}_{x}^{T}(\mathrm{\Delta}\mathit{\varphi})\mathit{g}.\end{array}$$

(38)

The fourth-order spatial tangent moduli
in (38) and the sensitivity of the Kirchhoff stresses to the action potential * ^{ϕ}* introduced in (37)

$${\mathbb{C}}^{\varphi \varphi}:=2{\partial}_{\mathit{g}}\widehat{\mathit{\tau}}(\mathit{g};\mathit{F},\mathrm{\Phi})\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathfrak{C}}^{\varphi \phi}:={\partial}_{\mathrm{\Phi}}\widehat{\mathit{\tau}}(\mathit{g};\mathit{F},\mathrm{\Phi}),$$

(39)

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

$$\mathrm{\Delta}{G}_{int}^{\varphi}={\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathit{\varphi}):{\mathbb{C}}^{\varphi \varphi}:(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}))\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathit{\varphi}):({\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\widehat{\mathit{\tau}})\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathit{\varphi}):{\mathfrak{C}}^{\varphi \phi}\mathrm{\Delta}\mathrm{\Phi}\phantom{\rule{0.16667em}{0ex}}dV.$$

(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

$$\mathrm{\Delta}{G}_{\text{int}}^{\phi}={\int}_{\mathcal{B}}\delta \mathrm{\Phi}\frac{\mathrm{\Delta}\mathrm{\Phi}}{\mathrm{\Delta}t}+\mathrm{\Delta}({\nabla}_{x}\delta \mathrm{\Phi})\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{ie})\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathrm{\Phi})\xb7\mathrm{\Delta}({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{ie})\phantom{\rule{0.16667em}{0ex}}dV.$$

(41)

Analogous to (37)_{1}, linearization of * _{x}*(

$$\mathrm{\Delta}({\nabla}_{x}(\delta \mathrm{\Phi}))=-{\nabla}_{x}(\delta \mathrm{\Phi}){\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}).$$

(42)

Furthermore, based on the functional definition of the spatial potential fluxes, , and , we obtain

$$\begin{array}{l}\mathrm{\Delta}{\widehat{\mathfrak{q}}}_{i}={\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}{\widehat{\mathfrak{q}}}_{i}+{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\xb7{\widehat{\mathfrak{q}}}_{i}+{\mathfrak{D}}_{i}\xb7{\nabla}_{x}(\mathrm{\Delta}\mathrm{\Phi}),\\ \mathrm{\Delta}{\widehat{\mathfrak{q}}}_{ie}={\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}{\widehat{\mathfrak{q}}}_{ie}+{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\xb7{\widehat{\mathfrak{q}}}_{ie}+{\mathfrak{D}}_{i}\xb7{\nabla}_{x}(\mathrm{\Delta}{\mathrm{\Phi}}_{e}),\\ \mathrm{\Delta}{\widehat{\mathfrak{q}}}_{e}={\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}{\widehat{\mathfrak{q}}}_{e}+{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\xb7{\widehat{\mathfrak{q}}}_{e}+\mathfrak{D}\xb7{\nabla}_{x}(\mathrm{\Delta}{\mathrm{\Phi}}_{e}),\end{array}$$

(43)

where **£**_{Δ}* _{ϕ}*
denotes the objective Lie derivative of the potential flux
along the increment Δ

$${\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}{\widehat{\mathfrak{q}}}_{a}={\mathfrak{C}}_{a}\phantom{\rule{0.16667em}{0ex}}:\phantom{\rule{0.16667em}{0ex}}{\scriptstyle \frac{1}{2}}{\mathbf{\pounds}}_{\mathrm{\Delta}\mathit{\varphi}}\mathit{g}={\mathfrak{C}}_{a}\phantom{\rule{0.16667em}{0ex}}:\phantom{\rule{0.16667em}{0ex}}(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}))$$

(44)

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

$${\mathfrak{D}}_{i}:={\mathfrak{D}}_{i}(\mathit{g};\mathit{F},{\mathit{f}}_{0}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathfrak{D}:=\mathfrak{D}(\mathit{g};\mathit{F},{\mathit{f}}_{0})$$

(45)

and the third-order mixed moduli

$${\mathfrak{C}}_{a}:=2{\partial}_{\mathit{g}}{\widehat{\mathfrak{q}}}_{a}(\mathit{g};\mathit{F},\mathrm{\Phi})\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}a=\{i,ie,e\}$$

(46)

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

$$\mathrm{\Delta}{G}_{int}^{\phi}{\int}_{\mathcal{B}}\delta \mathrm{\Phi}\frac{\mathrm{\Delta}\mathrm{\Phi}}{\mathrm{\Delta}t}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathrm{\Phi})\xb7{\mathfrak{D}}_{i}\xb7({\nabla}_{x}(\mathrm{\Delta}\mathrm{\Phi})+{\nabla}_{x}(\mathrm{\Delta}{\mathrm{\Phi}}_{e}))\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta \mathrm{\Phi})\xb7{\mathfrak{C}}^{\phi \varphi}:(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}))\phantom{\rule{0.16667em}{0ex}}dV,$$

(47)

where * ^{ϕ}* =

$$\mathrm{\Delta}{G}_{\text{ext}}^{\phi}:={\int}_{\mathcal{B}}\delta \mathrm{\Phi}\mathrm{\Delta}\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathcal{F}}}^{\phi}\phantom{\rule{0.16667em}{0ex}}dV$$

(48)

where

$$\mathrm{\Delta}{\widehat{\mathcal{F}}}^{\phi}=\mathcal{H}:(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}))+\mathcal{H}\mathrm{\Delta}\mathrm{\Phi}$$

(49)

is the linearization of the scalar-valued function
. The tangent terms
and
in (48)_{2} are defined as

$$\mathcal{H}:=2{\partial}_{\mathit{g}}{\widehat{\mathcal{F}}}^{\phi}(\mathit{g};\mathit{F},\mathrm{\Phi})\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\mathcal{H}:={\partial}_{\mathrm{\Phi}}{\widehat{\mathcal{F}}}^{\phi}(\mathit{g};\mathit{F},\mathrm{\Phi}),$$

(50)

respectively. Based on the decomposed form introduced before, the scalar tangent term can be expressed as

$$\mathcal{H}={\mathcal{H}}_{\text{e}}+{\mathcal{H}}_{\text{m}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathcal{H}}_{\text{e}}:={\partial}_{\mathrm{\Phi}}{\widehat{\mathcal{F}}}_{\text{e}}^{\phi},\phantom{\rule{0.16667em}{0ex}}{\mathcal{H}}_{\text{m}}:={\partial}_{\mathrm{\Phi}}{\widehat{\mathcal{F}}}_{\text{m}}^{\phi}.$$

(51)

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

$$\mathrm{\Delta}{G}_{\text{ext}}^{\phi}={\int}_{\mathcal{B}}\delta \mathrm{\Phi}(\mathcal{H}:(\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi}))+\mathcal{H}\mathrm{\Delta}\mathrm{\Phi})\phantom{\rule{0.16667em}{0ex}}dV.$$

(52)

Considering the weak form of the elliptic part of the bidomain equations ${G}_{\text{int}}^{{\phi}_{e}}$ in (33), the increment $\mathrm{\Delta}{G}_{\text{int}}^{{\phi}_{e}}$ can be written as

$$\mathrm{\Delta}{G}_{\text{int}}^{{\phi}_{e}}={\int}_{\mathcal{B}}\mathrm{\Delta}({\nabla}_{x}\delta {\mathrm{\Phi}}_{e})\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{e})+{\nabla}_{x}(\delta {\mathrm{\Phi}}_{e})\xb7\mathrm{\Delta}({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{e})\phantom{\rule{0.16667em}{0ex}}dV.$$

(53)

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

$$\mathrm{\Delta}{G}_{\text{int}}^{{\phi}_{e}}={\int}_{\mathcal{B}}{\nabla}_{x}(\delta {\mathrm{\Phi}}_{e})\xb7({\mathfrak{D}}_{i}\xb7{\nabla}_{x}(\mathrm{\Delta}\mathrm{\Phi})+\mathfrak{D}\xb7{\nabla}_{x}(\mathrm{\Delta}{\mathrm{\Phi}}_{e}))\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{\mathcal{B}}{\nabla}_{x}(\delta {\mathrm{\Phi}}_{e})\xb7{\mathfrak{C}}^{{\phi}_{e}\varphi}:\mathit{g}{\nabla}_{x}(\mathrm{\Delta}\mathit{\varphi})\phantom{\rule{0.16667em}{0ex}}dV,$$

(54)

where ^{eϕ} = * _{i}* +

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
into element subdomains
such that
$\mathcal{B}\approx {\cup}_{e=1}^{{n}_{el}}{\mathcal{B}}_{el}^{h}$ with *n _{el}* denoting the number of the subdomains
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
-continuous shape functions

$$\begin{array}{ll}{\mathit{\varphi}}^{h}={\sum}_{i=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{i}\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathit{x}}}_{i}& \delta {\mathit{\varphi}}^{h}={\sum}_{i=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{i}\phantom{\rule{0.16667em}{0ex}}\delta {\widehat{\mathit{x}}}_{i},\\ {\mathrm{\Phi}}^{h}={\sum}_{j=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathrm{\Phi}}}_{j}& \delta {\mathrm{\Phi}}^{h}={\sum}_{j=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}\delta {\widehat{\mathrm{\Phi}}}_{j},\\ {\mathrm{\Phi}}_{e}^{h}={\sum}_{k=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{k}\phantom{\rule{0.16667em}{0ex}}{\widehat{\mathrm{\Phi}}}_{ek}& \delta {\mathrm{\Phi}}_{e}^{h}={\sum}_{k=1}^{{\text{n}}_{\text{en}}}{\mathcal{N}}^{k}\phantom{\rule{0.16667em}{0ex}}\delta {\widehat{\mathrm{\Phi}}}_{ek},\end{array}$$

(55)

where *n _{en}* 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

$$\begin{array}{ll}{\nabla}_{x}(\delta {\mathit{\varphi}}^{h})={\sum}_{i=1}^{{\text{n}}_{\text{en}}}\delta {\widehat{\mathit{x}}}_{i}\otimes {\nabla}_{x}{\mathcal{N}}^{i}& {\nabla}_{x}\mathrm{\Delta}{\mathit{\varphi}}^{h}={\sum}_{i=1}^{{\text{n}}_{\text{en}}}\mathrm{\Delta}{\widehat{\mathit{x}}}_{i}\otimes {\nabla}_{x}{\mathcal{N}}^{i}\\ {\nabla}_{x}(\delta {\mathrm{\Phi}}^{h})={\sum}_{j=1}^{{\text{n}}_{\text{en}}}\delta {\widehat{\phi}}_{j}\otimes {\nabla}_{x}{\mathcal{N}}^{j}& {\nabla}_{x}\mathrm{\Delta}{\mathrm{\Phi}}^{h}={\sum}_{j=1}^{{\text{n}}_{\text{en}}}\mathrm{\Delta}{\widehat{\mathrm{\Phi}}}_{j}\otimes {\nabla}_{x}{\mathcal{N}}^{j}\\ {\nabla}_{x}(\delta {\mathrm{\Phi}}_{e}^{h})={\sum}_{k=1}^{{\text{n}}_{\text{en}}}\delta {\widehat{\mathrm{\Phi}}}_{ek}\otimes {\nabla}_{x}{\mathcal{N}}^{k}& {\nabla}_{x}\mathrm{\Delta}{\mathrm{\Phi}}_{e}^{h}={\sum}_{k=1}^{{\text{n}}_{\text{en}}}\mathrm{\Delta}{\widehat{\mathrm{\Phi}}}_{ek}\otimes {\nabla}_{x}{\mathcal{N}}^{k}.\end{array}$$

(56)

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

$$\begin{array}{l}{\mathbb{R}}^{\varphi}={\text{A}}_{e=1}^{{\text{n}}_{\text{el}}}{\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}{\mathcal{N}}^{i}\xb7\widehat{\tau}-{\mathcal{N}}^{i}\mathit{b}\phantom{\rule{0.16667em}{0ex}}dV-{\int}_{\partial {\mathcal{S}}_{t}^{e}}{\mathcal{N}}^{i}\xb7\overline{\mathit{t}}\phantom{\rule{0.16667em}{0ex}}dA=\mathbf{0},\\ {\mathbb{R}}^{\phi}={\text{A}}_{e=1}^{{\text{n}}_{\text{el}}}{\int}_{{\mathcal{B}}_{el}^{h}}{\mathcal{N}}^{i}\frac{\mathrm{\Phi}-{\mathrm{\Phi}}_{n}}{\mathrm{\Delta}t}+{\nabla}_{x}{\mathcal{N}}^{i}\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{ie})\phantom{\rule{0.16667em}{0ex}}dV-{\int}_{{\mathcal{B}}_{el}^{h}}{\mathcal{N}}^{i}{\widehat{\mathcal{F}}}^{\phi}\phantom{\rule{0.16667em}{0ex}}dV-{\int}_{\partial {\mathcal{B}}^{el}}{\mathcal{N}}^{i}\overline{\mathfrak{q}}\phantom{\rule{0.16667em}{0ex}}dA=0,\\ {\mathbb{R}}^{{\phi}_{e}}={\text{A}}_{e=1}^{{\text{n}}_{\text{el}}}{\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}{\mathcal{N}}^{i}\xb7({\widehat{\mathfrak{q}}}_{i}+{\widehat{\mathfrak{q}}}_{e})\phantom{\rule{0.16667em}{0ex}}dV-{\int}_{\partial {\mathcal{B}}_{el}^{h}}{\mathcal{N}}^{i}{\overline{\mathfrak{q}}}_{t}\phantom{\rule{0.16667em}{0ex}}dA=0\end{array}$$

(57)

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

$$\text{Lin}\phantom{\rule{0.16667em}{0ex}}\mathbb{R}=\mathbb{R}+\frac{\partial \mathbb{R}}{\partial \mathbb{U}}\mathrm{\Delta}\mathbb{U}\phantom{\rule{0.16667em}{0ex}};\phantom{\rule{0.16667em}{0ex}}\mathbb{R}=\left\{\begin{array}{l}{\mathbb{R}}^{\varphi}\hfill \\ {\mathbb{R}}^{\phi}\hfill \\ {\mathbb{R}}^{{\phi}_{e}}\hfill \end{array}\right\}\phantom{\rule{0.16667em}{0ex}};\phantom{\rule{0.16667em}{0ex}}\mathbb{U}=\left\{\begin{array}{l}{\widehat{\mathit{x}}}^{h}\hfill \\ {\widehat{\mathrm{\Phi}}}^{h}\hfill \\ {\widehat{\mathrm{\Phi}}}_{e}^{h}\hfill \end{array}\right\}.$$

(58)

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

$$\mathbb{K}=\frac{\partial \mathbb{R}}{\partial \mathbb{U}}=\underset{e=1}{\overset{{n}_{el}}{\text{A}}}{\mathbb{K}}_{el}\phantom{\rule{0.16667em}{0ex}};\phantom{\rule{0.16667em}{0ex}}{\mathbb{K}}_{el}:=\left[\begin{array}{lll}{\mathbb{K}}^{\varphi \varphi}\hfill & {\mathbb{K}}^{\varphi \phi}\hfill & {\mathbb{K}}^{\varphi {\phi}_{e}}\hfill \\ {\mathbb{K}}^{\phi \varphi}\hfill & {\mathbb{K}}^{\phi \phi}\hfill & {\mathbb{K}}^{\phi {\phi}_{e}}\hfill \\ {\mathbb{K}}^{{\phi}_{e}\varphi}\hfill & {\mathbb{K}}^{{\phi}_{e}\phi}\hfill & {\mathbb{K}}^{{\phi}_{e}{\phi}_{e}},\hfill \end{array}\right]$$

(59)

where the components of element matrix are defined as

$$\begin{array}{l}{\mathbb{K}}^{\varphi \varphi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}\xb7\mathbb{C}\xb7{\nabla}_{x}{\mathcal{N}}^{j}+{\nabla}_{x}{\mathcal{N}}^{i}\xb7\widehat{\tau}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{\varphi \phi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}\xb7{\mathfrak{C}}^{\varphi \phi}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{\varphi {\phi}_{e}}=\mathbf{0},\\ {\mathbb{K}}^{\phi \varphi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}{\mathcal{N}}^{i}\xb7{\mathfrak{C}}^{\phi \varphi}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV+{\int}_{{\mathcal{B}}_{el}^{h}}{\mathcal{N}}^{i}\mathcal{H}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV\\ {\mathbb{K}}^{\phi \phi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}\xb7{\mathfrak{D}}_{i}\xb7{\nabla}_{x}{\mathcal{N}}^{j}+{\mathcal{N}}^{i}\left(\frac{1}{\mathrm{\Delta}t}-\mathcal{H}\right){\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{\phi {\phi}_{e}}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}.{\mathfrak{D}}_{i}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{{\phi}_{e}\varphi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}{\mathcal{N}}^{i}.{\mathfrak{C}}^{{\phi}_{e}\varphi}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{{\phi}_{e}\phi}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}.{\mathfrak{D}}_{i}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV,\\ {\mathbb{K}}^{{\phi}_{e}{\phi}_{e}}={\int}_{{\mathcal{B}}_{el}^{h}}{\nabla}_{x}^{T}{\mathcal{N}}^{i}.\mathfrak{D}\xb7{\nabla}_{x}{\mathcal{N}}^{j}\phantom{\rule{0.16667em}{0ex}}dV.\end{array}$$

(60)

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 **, the potential fluxes *** _{i}*,

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 ** C̄** and the two unit vectors

$$\begin{array}{ll}{\overline{I}}_{1}:=\text{tr}(\overline{\mathit{C}}),\hfill & {\overline{I}}_{4\text{f}}:=\overline{\mathit{C}}:({\mathit{f}}_{0}\otimes {\mathit{f}}_{0}),\hfill \\ {\overline{I}}_{4\text{s}}:=\overline{\mathit{C}}:({\mathit{s}}_{0}\otimes {\mathit{s}}_{0}),\hfill & {\overline{I}}_{8\text{fs}}:=\overline{\mathit{C}}:\text{sym}({\mathit{f}}_{0}\otimes {\mathit{s}}_{0}).\hfill \end{array}$$

(61)

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

$$\mathrm{\Psi}=U(J)+\overline{\mathrm{\Psi}}({\overline{I}}_{1},{\overline{I}}_{4\text{f}},{\overline{I}}_{4\text{s}},{\overline{I}}_{8\text{fs}}).$$

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

$$\overline{\mathrm{\Psi}}=\frac{a}{2b}exp[b({\overline{I}}_{1}-3)]+\frac{{a}_{\text{fs}}}{2{b}_{\text{fs}}}\{exp[{b}_{\text{fs}}{\overline{I}}_{8\text{fs}}^{2}]-1\}+\sum _{i=\text{f},\text{s}}\frac{{a}_{i}}{2{b}_{i}}\{exp[{b}_{i}{({\overline{I}}_{4i}-1)}^{2}]-1\}.$$

(63)

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

$$\mathit{\tau}={\mathit{\tau}}_{\text{vol}}+{\mathit{\tau}}_{\text{iso}}=\widehat{p}\phantom{\rule{0.16667em}{0ex}}{\mathit{g}}^{\u266f}+\overline{\mathit{\tau}}:\mathbb{P}$$

(64)

in terms of the pressure := *JU*′(*J*), the Kirchhoff stress **:= 2 ****_{g}** and the contravariant metric tensor in the Eulerian manifold

$$\overline{\mathit{\tau}}=2{\overline{\mathrm{\Psi}}}_{1}\overline{\mathit{b}}+2{\overline{\mathrm{\Psi}}}_{4\text{f}}\mathit{f}\otimes \mathit{f}+2{\overline{\mathrm{\Psi}}}_{4\text{s}}\phantom{\rule{0.16667em}{0ex}}\mathit{s}\otimes \mathit{s}+2{\overline{\mathrm{\Psi}}}_{8\text{fs}}\phantom{\rule{0.16667em}{0ex}}\text{sym}(\mathit{f}\otimes \mathit{s}),$$

(65)

where the deformation-dependent scalar coefficients _{1}, _{4f}, _{4s}, and _{8fs} are defined as follows

$$\begin{array}{l}{\overline{\mathrm{\Psi}}}_{1}:={\partial}_{{\overline{I}}_{1}}\mathrm{\Psi}={\scriptstyle \frac{a}{2}}exp[b({\overline{I}}_{1}-3)],\\ {\overline{\mathrm{\Psi}}}_{4\text{f}}:={\partial}_{{\overline{I}}_{4\text{f}}}\mathrm{\Psi}={a}_{\text{f}}({\overline{I}}_{4\text{f}}-1)exp[{b}_{\text{f}}{({\overline{I}}_{4\text{f}}-1)}^{2}],\\ {\overline{\mathrm{\Psi}}}_{4\text{s}}:={\partial}_{{\overline{I}}_{4\text{s}}}\mathrm{\Psi}={a}_{\text{s}}({\overline{I}}_{4\text{s}}-1)exp[{b}_{\text{s}}{({\overline{I}}_{4\text{f}}-1)}^{2}],\\ {\overline{\mathrm{\Psi}}}_{8\text{fs}}:={\partial}_{{\overline{I}}_{8\text{fs}}}\mathrm{\Psi}={a}_{\text{fs}}{\overline{I}}_{8\text{fs}}exp[{b}_{\text{fs}}{\overline{I}}_{8\text{fs}}^{2}].\end{array}$$

(66)

The active Kirchhoff stress ^{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

$${\widehat{\mathit{\tau}}}^{\text{act}}(\mathit{g};\mathrm{\Phi},\mathit{f})=\sigma (\mathrm{\Phi})\phantom{\rule{0.16667em}{0ex}}\mathit{f}\otimes \mathit{f}.$$

(67)

The above equation implies that the direction of the active stress tensor is dictated by the deformed structural tensor *f*** f**, while its magnitude is chiefly determined by the transmembrane potential-dependent active fiber tension

$$\stackrel{.}{\sigma}=\epsilon (\mathrm{\Phi})[{k}_{\sigma}(\mathrm{\Phi}-{\mathrm{\Phi}}_{r})-\sigma ],$$

(68)

in which the parameter *k _{σ}* controls the saturated value of

$$\epsilon (\mathrm{\Phi})={\epsilon}_{0}+({\epsilon}_{\infty}-{\epsilon}_{0})exp[-exp(-\xi (\mathrm{\Phi}-\mathrm{\Phi}))]$$

(69)

in terms of the parameters *ε*_{0} and *ε*_{∞} that characterize the two limiting values of the function for Φ < and Φ > about the phase shift , respectively. The transition rate of *ε* from *ε*_{0} to *ε*_{∞} about 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

$$\sigma (\mathrm{\Phi})=\frac{1}{1+\mathrm{\Delta}t\phantom{\rule{0.16667em}{0ex}}\epsilon (\mathrm{\Phi})}[{\sigma}_{n}+\mathrm{\Delta}t\phantom{\rule{0.16667em}{0ex}}\epsilon (\mathrm{\Phi}){k}_{\sigma}(\mathrm{\Phi}-{\mathrm{\Phi}}_{r})]$$

(70)

in terms of the current action potential Φ and the fiber tension *σ _{n}* of the previous time step

$${\mathbb{C}}^{\varphi \varphi}={\mathbb{C}}^{\text{vol}}+{\mathbb{C}}^{\text{iso}}.$$

(71)

The volumetric part is given by the common expression

$${\mathbb{C}}^{\text{vol}}=(\widehat{p}+\widehat{\kappa}){\mathit{g}}^{\u266f}\otimes {\mathit{g}}^{\u266f}-2p\mathbb{I}$$

(72)

where := *J*^{2}*U*″(*J*) is the bulk modulus. The symmetric fourth identity tensor has the indicial representation
${\mathbb{I}}^{\mathit{ijkl}}:={\scriptstyle \frac{1}{2}}({g}^{ik}{g}^{jl}+{g}^{il}{g}^{jk})$ in terms of the components of the contravariant metric *g*^{}. The isochoric part of the moduli has the following form

$${\mathbb{C}}^{\text{iso}}=\mathbb{P}:[\overline{\mathbb{C}}+{\scriptstyle \frac{2}{3}}(\overline{\mathit{\tau}}:\mathit{g})\mathbb{I}-{\scriptstyle \frac{2}{3}}(\mathbb{P}:\overline{\mathit{\tau}}\otimes {\mathit{g}}^{\u266f}+{\mathit{g}}^{\u266f}\otimes \overline{\mathit{\tau}})]:\mathbb{P}$$

(73)

where the unimodular tangent term is derived from the stress expression (65)

$$\overline{\mathbb{C}}=4{\overline{\mathrm{\Psi}}}_{1}^{\prime}(\overline{\mathit{b}}\otimes \overline{\mathit{b}})+4{\overline{\mathrm{\Psi}}}_{4\text{f}}^{\prime}\mathbb{F}+4{\overline{\mathrm{\Psi}}}_{4\text{s}}^{\prime}\mathbb{S}+4{\overline{\mathrm{\Psi}}}_{8\text{fs}}^{\prime}\phantom{\rule{0.16667em}{0ex}}\text{sym}(\mathit{f}\otimes \mathit{s})\otimes \text{sym}(\mathit{f}\otimes \mathit{s})$$

(74)

in terms of the fourth order structural tensors

$$\mathbb{F}:=\mathit{f}\otimes \mathit{f}\otimes \mathit{f}\otimes \mathit{f}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathbb{S}:=\mathit{s}\otimes \mathit{s}\otimes \mathit{s}\otimes \mathit{s}.$$

(75)

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

$${\mathfrak{C}}^{\varphi \phi}={\sigma}^{\prime}(\mathrm{\Phi})\mathit{f}\otimes \mathit{f}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\sigma}^{\prime}(\mathrm{\Phi}):={\partial}_{\mathrm{\Phi}}\sigma (\mathrm{\Phi}).$$

(76)

We assume that the second-order conduction tensors
in the intra- and extracellular domains *a*:= {*i*,*e*} can be additively decomposed

$${\mathfrak{D}}^{a}={d}_{\left|\right|}^{a}\mathit{f}\otimes \mathit{f}+{d}_{\perp}^{a}({\mathit{g}}^{\u266f}-\mathit{f}\otimes \mathit{f})$$

(77)

into the conductance along the fiber direction
${d}_{\left|\right|}^{a}$ and the transverse conductance across the fiber direction
${d}_{\perp}^{a}$, which can be determined experimentally using microelectrode array recordings [7]. The third-order mixed moduli * ^{a}* which is defined in (46) and enter the formulations in (47) and (54) follows from inseting (46) into (32)

$$\begin{array}{l}{\mathfrak{C}}_{i}=-2\phantom{\rule{0.16667em}{0ex}}{d}_{\perp}^{i}{\nabla}_{x}\mathrm{\Phi}\xb7\mathbb{I},\\ {\mathfrak{C}}_{ie}=-2\phantom{\rule{0.16667em}{0ex}}{d}_{\perp}^{i}{\nabla}_{x}{\mathrm{\Phi}}_{e}\xb7\mathbb{I},\\ {\mathfrak{C}}_{e}=-2\phantom{\rule{0.16667em}{0ex}}({d}_{\perp}^{i}+{d}_{\perp}^{e}){\nabla}_{x}{\mathrm{\Phi}}_{e}\xb7\mathbb{I}.\end{array}$$

(78)

Last, we specify the constitutive equation for the electrical source term
. To set up the model equations and parameters in the non-dimensional space, we introduce the non-dimensional transmembrane potential and the non-dimensional time *τ* through the following conversion formulae

$$\mathrm{\Phi}={\beta}_{\phi}\phi -{\delta}_{\phi}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}t={\beta}_{t}\tau .$$

(79)

The non-dimensional potential is related to the physical transmembrane potential Φ through the factor *β _{}* and the potential difference

$${\widehat{\mathcal{F}}}^{\phi}=\frac{{\beta}_{\phi}}{{\beta}_{t}}{\widehat{f}}^{\phi},\mathcal{H}=\frac{{\beta}_{\phi}}{{\beta}_{t}}\mathit{h}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mathcal{H}=\frac{1}{{\beta}_{t}}h$$

(80)

for the normalized source term * ^{}*, and the non-dimensional counterparts

$$h={h}_{\text{e}}+{h}_{\text{m}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{with}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{h}_{\text{e}}:={\partial}_{\phi}{\widehat{f}}_{\text{e}}^{\phi},{h}_{\text{m}}:={\partial}_{\phi}{\widehat{f}}_{\text{m}}^{\phi}.$$

(81)

Throughout the subsequent treatments, the Aliev-Panfilov model [2] is used for the description of the source term
${\widehat{f}}_{\text{e}}^{\phi}$ and the evolution of the recovery variable *r*, following the algorithmic treatment introduced in [17, 18].

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.

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**,

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

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**,

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

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

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.

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.

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.

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

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. [PMC free article] [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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |