Semi-implicit Non-conforming Finite-Element Schemes for Cardiac Electrophysiology: A Framework for Mesh-Coarsening Heart Simulations

The field of computational cardiology has steadily progressed toward reliable and accurate simulations of the heart, showing great potential in clinical applications such as the optimization of cardiac interventions and the study of pro-arrhythmic effects of drugs in humans, among others. However, the computational effort demanded by in-silico studies of the heart remains challenging, highlighting the need of novel numerical methods that can improve the efficiency of simulations while targeting an acceptable accuracy. In this work, we propose a semi-implicit non-conforming finite-element scheme (SINCFES) suitable for cardiac electrophysiology simulations. The accuracy and efficiency of the proposed scheme are assessed by means of numerical simulations of the electrical excitation and propagation in regular and biventricular geometries. We show that the SINCFES allows for coarse-mesh simulations that reduce the computation time when compared to fine-mesh models while delivering wavefront shapes and conduction velocities that are more accurate than those predicted by traditional finite-element formulations based on the same coarse mesh, thus improving the accuracy-efficiency trade-off of cardiac simulations.

1. Introduction

Computer simulations of the electrical activity of the heart have increasingly gained attention in the medical community, as they have steadily shown potential in the study of cardiac diseases and in the design of novel cardiac therapies. Current models of the human heart are able to represent the complex three-dimensional anatomical structure of the heart chambers, incorporating key functional features such as the Purkinje network and the cardiomyocyte orientation (Vadakkumpadan et al., 2009). Such advanced representation of the heart has enabled novel in-silico studies of undesired pro-arrhythmic effects of drugs in patients (Sahli Costabal et al., 2018), potentially reducing the number of subjects needed in a clinical trial by aiding the experiment design. Computational models of the heart have also shown promise in assisting the design of effective therapies for terminating atrial fibrillation (Trayanova et al., 2018). While these examples can only confirm the tremendous relevance of computational models in advancing the field of cardiology, they share the fundamental challenge of being highly demanding in terms of wall-clock time needed in computer simulations.

Mathematical models of the heart require the computer implementation of spatio-temporal discretization techniques in order to obtain a sequence of numerical representations of the physiological fields under study. Two fundamental aspects directly responsible for the computation time (CT) in a heart simulation are the ionic model used to account for subcellular electrochemical mechanisms, and the level of spatio-temporal discretization in terms of time-step size and mesh size (Sundnes et al., 2006). The choice of the mesh size typically faces a well-known trade-off problem of accuracy vs. efficiency, as decreasing the mesh size in a simulation results in more accurate numerical approximations, at the cost of increasing the number of degrees of freedom (DOFs), which drives the CT. Indeed, current simulations of the heart typically employ mesh sizes in the range of tens to hundreds of micrometers for domains with lengths in the order of centimeters, which ultimately translates into large systems of equations with several millions of DOFs that need to be solved at each time step. Such high dimensionality renders the solution of heart simulations extremely challenging for personal computers, and calls for improving their implementation in high-performance computing (HPC) platforms (Niederer et al., 2011a; Vazquez et al., 2011).

In the particular case of cardiac electrophysiology simulations, a common criterion to select the mesh size is the ability of the numerical simulation to recover an accurate conduction velocity (CV) and wavefront shape (Pathmanathan et al., 2010; Krishnamoorthi et al., 2013; Dupraz et al., 2015). It has been shown that both the wavefront shape and the CV suffer from a strong dependence on the spatial discretization, which for the case of finite-element (FE) discretization using linear basis functions results in a significant loss of accuracy for the case of mesh sizes >0.1mm (Pezzuto et al., 2016). In order to achieve larger mesh sizes, higher-order FE formulations have been proposed, which show that FE Lagrange basis functions of order 2 and 3 result in accurate CV for coarser meshes (Arthurs et al., 2012; Pezzuto et al., 2016). It should be noted, however, that higher-order FE schemes based on Lagrange basis functions necessarily increase the total number of DOFs in simulations when compared to linear-element formulations, as well as they require an additional computational effort for quadrature procedures, as higher-order basis functions demand the use of more quadrature-point evaluations (Cantwell et al., 2014). Recently, Hurtado and Rojas (2018) introduced a non-conforming finite-element scheme (NCFES) for the spatial discretization of the monodomain equation of cardiac electrophysiology that allows for the use of coarse meshes without significant loss of accuracy measured in terms of CV and wavefront shape. More specifically, hexahedral trilinear elements (Q1) were enhanced with non-conforming basis functions of degree 2 to create a non-conforming element (Q1NC) that is capable of representing a second-order polynomial within the element domain, a concept widely employed in the context of solid mechanics FE simulations (Wilson et al., 1973; Taylor et al., 1976). Further, they showed that the DOFs associated to the non-conforming basis functions can be solved at the element level, and therefore the number of global DOFs of the Q1NC scheme equals that of a standard Q1 FE scheme. As a result, Q1NC simulations delivered a CV and wavefront shape similar to that of second-order Lagrange formulations (Q2) at the computational cost in the order of a Q1 formulation.

During the development of the NCFES for cardiac electrophysiology, a fully-implicit (FI) backward-Euler time-stepping method was considered (Hurtado and Rojas, 2018). While FI schemes have important advantages in delivering a larger time-step stability region in cardiac simulations (Ying et al., 2008; Hurtado and Henao, 2014), they require the solution of a large system of non-linear equations at each time step that can be very costly in computational terms, and may not be well-suited to parallel-computing platforms when compared to other numerical schemes. To improve the computational efficiency, the semi-implicit integration method has been proposed in the literature for solving the semi-discrete equations resulting from standard FE discretizations, showing a relevant decrease in the CT of cardiac simulations, as well as being amenable to HPC platforms (Whiteley, 2006; Pathmanathan et al., 2010). Consequently, the scientific question that motivates this work is: Can we further improve the efficiency-accuracy trade-off in cardiac simulations by combining non-conforming FE spatial discretizations with semi-implicit time-integration schemes? To answer such question, in the following we develop the numerical framework and present an algorithm for the implementation of a semi-implicit non-conforming FE scheme to solve the monodomain electrophysiology equations, and investigate the numerical consequences and potential contributions to cardiac simulations.

2. Methods

2.1. Monodomain Model of Cardiac Electrophysiology

Let Ω ∈ ℝ³ be the heart domain where electrical impulses travel during the time interval [0, T], and V_m:Ω × [0, T] → ℝ be the transmembrane potential. A local statement of current balance yields the monodomain equation (Pullan et al., 2005)

$\begin{array}{ll}{A}_{\text{m}}\left(C_{\text{m}}\frac{\partial V_{\text{m}}}{\partial t}+I_{\text{ion}}\left(V_{\text{m}},r\right)\right)-\text{div}(\sigma \nabla V_{\text{m}})=0,\mathrm{\text{in}}\Omega \times (0,T],& (1)\end{array}$

where A_m, C_m are the surface-to-volume ratio and membrane capacitance, respectively, σ is the conductivity tensor, I_ion is the ionic current depending on the transmembrane potential V_m, and $r:\Omega \times (0,T]\to {ℝ}^{\text{m}}$ is a vector field of state variables that include gating variables and ion concentrations. For convenience, we consider the normalized transmembrane potential field

$\begin{array}{l}\varphi (x,t)=\frac{V_{\text{m}}(x,t)-V_{\text{r}}}{V_{\text{p}}-V_{\text{r}}},\end{array}$

where V_p and V_r are the peak and resting voltages, respectively. Based on this normalization, we obtain the non-dimensional monodomain equation,

$\begin{array}{ll}\frac{\partial \varphi }{\partial t}-\text{div}(D\nabla \varphi )-f(\varphi ,r)=0\mathrm{\text{in}}\Omega \times (0,T],& (2)\end{array}$

where $D=\frac{1}{A_{\text{m}}{C}_{\text{m}}}\sigma$ is the normalized conductivity tensor, and $f(\varphi ,r)=-\frac{I_{\text{ion}}(V_{\text{m}}(\varphi ),r)}{C_{\text{m}}(V_{\text{p}}-V_{\text{r}})}$ is the normalized ionic current. The time evolution of state variables is governed by kinetic equations of the form

$\begin{array}{ll}\frac{\partial r}{\partial t}=g(\varphi ,r)\mathrm{\text{in}}\Omega \times (0,T].& (3)\end{array}$

The expressions for $f(\varphi ,r)$ and $g(\varphi ,r)$ will depend on the choice of ionic model representing the transmembrane ionic current in a single cell. Equations (2, 3) are complemented with Dirichlet and Neumann boundary conditions,

$\begin{array}{ll}\varphi =\bar{\varphi },\mathrm{\text{on}}\partial {\Omega }_{\varphi }\times (0,T],& (4)\end{array}$

$\begin{array}{ll}q·n=\bar{q},\mathrm{\text{on}}\partial {\Omega }_q\times (0,T],& (5)\end{array}$

respectively, as well as initial conditions

$\begin{array}{l}\varphi (x,0)={\varphi }_0(x),x\in \Omega ,\\ r(x,0)=r_0(x),x\in \Omega .\end{array}$

To state the weak form of the cardiac electrophysiology problem, we consider trial spaces ${{S}}^{\varphi },{{S}}^r$ and test spaces ${{V}}^{\varphi },{{V}}^r$ defined as

$\begin{array}{ll}{{S}}^{\varphi }=\left\{\varphi \in L^2\left(\left(0,T\right];H^1\left(\Omega ,ℝ\right)\right):\varphi =\bar{\varphi }\mathrm{\text{on}}\partial {\Omega }_{\varphi }\times \left(0,T\right]\right\}& (6)\end{array}$

$\begin{array}{ll}{{S}}^r=\left\{r\in L^2\left(\left(0,T\right];L^2\left(\Omega ,{ℝ}^m\right)\right)\right\}& \left(7\right)\end{array}$

$\begin{array}{ll}{{V}}^{\varphi }=\left\{\nu \in H^1\left(\Omega ,ℝ\right):\nu =0\mathrm{\text{on}}\partial {\Omega }_{\varphi }\right\}& \left(8\right)\end{array}$

$\begin{array}{ll}{{V}}^r=\left\{\eta \in L^2\left(\Omega ,{ℝ}^m\right)\right\}& \left(9\right)\end{array}$

Multiplying (2) and (3) by appropriate test functions, integrating over Ω and applying the divergence theorem yields the weak equations, and the statement of the weak formulation reads: ∀t ∈ (0, T], find $(\varphi ,r)\in {{S}}^{\varphi }\times {{S}}^r$ such that

$\begin{array}{l}{G}^{\varphi }\left[(\varphi ,r),(\nu ,\eta )\right]:={\displaystyle {\int }_{\Omega }}\nu \frac{\partial \varphi }{\partial t}\mathrm{\text{d}}x+{\displaystyle {\int }_{\Omega }}\nabla \nu ·\mathrm{\text{D}}\nabla \varphi \mathrm{\text{d}}x\\ -{\displaystyle {\int }_{\Omega }}\nu f(\varphi ,r)\mathrm{\text{d}}x+{\displaystyle {\int }_{\partial {\Omega }_q}}\nu \bar{q}\mathrm{\text{d}}s\\ =0,\forall \nu \in {{V}}^{\varphi }& (10)\end{array}$

$\begin{array}{l}{G}^r\left[(\varphi ,r),(\nu ,\eta )\right]:={\displaystyle {\int }_{\Omega }}\eta \left\{\frac{\partial r}{\partial t}-g(\varphi ,r)\right\}\mathrm{\text{d}}x\\ =0,\forall \eta \in {{V}}^r& (11)\end{array}$

2.2. Spatial Discretization Using a Non-conforming Finite-Element Scheme

A Galerkin finite-element scheme to solve the weak formulation of the monodomain problem can be stated as follows. Let ${\Omega }^h={\cup }_{e=1}^{N_{\text{el}}}{\Omega }_e$ be a domain discretization where N_el is the number of elements, and all elements comply with the condition Ω_i∩Ω_j = ∅ for i ≠ j. We construct finite-dimensional subspaces ${{S}}_h^{\varphi }\subset {{S}}^{\varphi }$, ${{S}}_h^r\subset {{S}}^r$ and ${{V}}_h^{\varphi }\subset {{V}}^{\varphi }$, ${{V}}_h^r\subset {{V}}^r$, to solve the following FE problem (Göktepe and Kuhl, 2009; Hurtado and Kuhl, 2014): ∀t ∈ (0, T], find $({\varphi }^h,r^h)\in {{S}}_h^{\varphi }\times {{S}}_h^r$ such that

$\begin{array}{l}{G}^{\varphi }\left[({\varphi }^h,r^h),({\nu }^h,{\eta }^h)\right]=0,\forall {\nu }^h\in {{V}}_h^{\varphi }\\ G^r\left[({\varphi }^h,r^h),({\nu }^h,{\eta }^h)\right]=0,\forall {\eta }^h\in {{V}}^r.\end{array}$

A traditional discretization FE scheme is the hexahedral isoparametric finite-element space,

$\begin{array}{l}{{V}}_h^{\varphi }:=\left\{{\nu }^h\in C^0({\Omega }^h,ℝ):{\nu }^h{|}_{{\Omega }_e}\in Q_k({\Omega }_e),e=1,\dots ,N_{\text{el}}\right\}\end{array}$

where Q_k(Ω_e) represents the space of isoparametric functions resulting from n-tensor product of 1-D Lagrange polynomials of order k, which are defined over the standard (isoparametric) domain $\widehat{\Omega }={\left[-1,1\right]}^n$ and mapping to a hexahedral element. We expand an element ${\nu }^h\in {{V}}_h^{\varphi }$ as

$\begin{array}{l}{\nu }^h(x)={\displaystyle \sum _{A=1}^{N_{\text{dofs}}}}{N}_A(x){\nu }_A,\end{array}$

where {N_A}a = 1, N_dofs are the basis functions, N_dofs is the number of element nodes with unknown degrees of freedom, and {ν_A}a = 1, N_dofs are the nodal coefficients. Using the same element basis functions, we expand the trial functions as

$\begin{array}{ll}{\varphi }^h(x,t)={\displaystyle \sum _{A=1}^{N_{\text{dofs}}}}{N}_A(x)u_A(t)+u_{\text{BC}}(x,t),& (12)\end{array}$

where {u_A(t)}A = 1, N_dofs correspond to the nodal values of the transmembrane potential field, and $u_{\text{BC}}\in {{S}}^{\varphi }$ is a function that satisfies the boundary conditions (4), i.e., $u_{\text{BC}}=\bar{\varphi }$ in ∂Ω_ϕ × (0, T]. For simplicity, and without loss of generality, in the following we assume that u_BC = 0. To construct the elements of ${{V}}_h^r$, we write

$\begin{array}{ll}{\eta }^h(x)={\displaystyle \sum _{e=1}^{N_{\text{el}}}}{\displaystyle \sum _{q=1}^{N_{\text{q}}}}{M}_q^e(x){\eta }_q^e,& (13)\end{array}$

where $M_q^e$ is a characteristic function defined by

$\begin{array}{ll}{M}_q^e(x)=\{\begin{array}{c}1,x\in {\Omega }_{e,q}\\ 0,x\notin {\Omega }_{e,q}\end{array}& (14)\end{array}$

and Ω_e,q ⊂ Ω_e is the subdomain containing the q−quadrature point $x_q$, and is such that ${\bigcup }_{q=1}^{Nq}{\Omega }_{e,q}={\Omega }_e$ and ${\Omega }_{e,q}\cap {\Omega }_{e,q^{\prime }}=\varnothing$ whenever q ≠ q′. Analogously, we expand an element $r^h\in {{S}}_h^r$ as

$\begin{array}{ll}{r}^h(x,t)={\displaystyle \sum _{e=1}^{N_{\text{el}}}}{\displaystyle \sum _{q=1}^{N_{\text{q}}}}{M}_q^e(x)r_q^e(t),& (15)\end{array}$

where $r_q^e:(0,T]\to {ℝ}^{\text{m}}$ represents the time evolution of the state variables at the q-quadrature point.

In this work, we consider a non-conforming spatial-discretization scheme for the monodomain equations (Hurtado and Rojas, 2018). To this end, we rewrite the residuals as

$\begin{array}{l}{G}^{\varphi }[(\varphi ,r),(\nu ,\eta )]={\displaystyle \sum _{e=1}^{N_{\text{el}}}\{{\displaystyle {\int }_{{\Omega }_e}\nu }\frac{\partial \varphi }{\partial t}\text{ d}x+{\displaystyle {\int }_{{\Omega }_e}\nabla }\nu ·\text{D}\nabla \varphi \text{ d}x}\\ -{\displaystyle {\int }_{{\Omega }_e}\nu }f(\varphi ,r)\text{ d}x+{\displaystyle {\int }_{\partial {\Omega }_{e,q}}\nu }\overline{q}\text{ d}s\},& (16)\end{array}$

$\begin{array}{ll}{G}^r\left[(\varphi ,r),(\nu ,\eta )\right]={\displaystyle \sum _{e=1}^{N_{\text{el}}}}\left\{{\displaystyle {\int }_{{\Omega }_e}}\eta \left\{\frac{\partial r}{\partial t}-g(\varphi ,r)\right\}\mathrm{\text{d}}x\right\},& (17)\end{array}$

and note that in such form, integrability of the trial and test functions and their weak derivatives is required only at the element level. We enhance the polynomial basis of ${{V}}_h^{\varphi }$ at the element level by adding polynomial terms not included in Q_k(Ω_e). To this end, we consider the non-conforming space

$\begin{array}{l}{{E}}_h^{\varphi }:=\left\{{\beta }^h:{\beta }^h{|}_{{\Omega }_e}\in P_{k+m}({\Omega }_e)\backslash Q_k({\Omega }_e)\right\}\end{array}$

where m ∈ ℤ₊ and P_k+m(Ω_e) is the space of polynomial functions of degree k+m defined on the standard domain $\widehat{\Omega }$. We then consider enhanced test functions ν^h which we expand as

$\begin{array}{ll}{\nu }^h(x)={\displaystyle \sum _{A=1}^{N_{\text{dofs}}}}{N}_A(x){\nu }_A+{\displaystyle \sum _{e=1}^{N_{\text{el}}}}{\displaystyle \sum _{c=1}^{N_{\text{nc}}}}{W}_c^e(x){\beta }_c^e& (18)\end{array}$

where ${\beta }_c^e\in ℝ$ are coefficients, $W_c^e$ are non-conforming element basis functions, and it holds that $W_c^e=0,x\notin {\Omega }^e$. Analogously, we enhance ${{S}}_h^{\varphi }$ with the time-dependent non-conforming space ${{F}}_h^{\varphi }$, and expand the enhanced trial functions as

$\begin{array}{ll}{\varphi }^h(x,t)=u^h(x,t)+{\alpha }^h(x,t)& (19)\end{array}$

where

$\begin{array}{ll}{u}^h(x,t):={\displaystyle \sum _{B=1}^{N_{\text{dofs}}}}{N}_B(x)u_B(t)& (20)\end{array}$

$\begin{array}{ll}{\alpha }^h(x,t):={\displaystyle \sum _{e=1}^{N_{\text{el}}}}{\displaystyle \sum _{d=1}^{N_{\text{nc}}}}{W}_d^e(x){\alpha }_d^e(t).& (21)\end{array}$

and ${\alpha }_d^e:(0,T]\to ℝ$ is a time-dependent coefficient that scales the non-conforming basis functions $W_d^e$. Substitution of approximations Equations (13 15, 18, and 19) into the residuals Equations (16) and (17) yields the following semi-discrete problem: ∀t ∈ (0, T], find $(u^h,{\alpha }^h,r^h)\in {{S}}_h^{\varphi }\times {{F}}_h^{\varphi }\times {{S}}_h^r$ such that

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{N}_A}\{{\dot{u}}^h+{\dot{\alpha }}^h\}\text{d}x+{\displaystyle {\int }_{\Omega }\nabla }{N}_A·\text{D}\nabla \{u^h+{\alpha }^h\}\text{d}x-{\displaystyle {\int }_{\Omega }{N}_A}f(u^h\\ +{\alpha }^h,r^h)\text{d}x+{\displaystyle {\int }_{\partial {\Omega }_q}{N}_A}\overline{q}\text{ d}s=0,A=1,\dots ,N_{\text{dofs}},& (22)\end{array}$

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }^e}}{W}_c^e\left\{{\dot{u}}^h+{\dot{\alpha }}^h\right\}\mathrm{\text{d}}x+{\displaystyle {\int }_{{\Omega }^e}}\nabla W_c^e·\text{D}\nabla \left\{u^h+{\alpha }^h\right\}\mathrm{\text{d}}x\\ -{\displaystyle {\int }_{{\Omega }^e}}{W}_c^{e}f(u^h+{\alpha }^h,r^h)\mathrm{\text{d}}x=0,e=1,\dots ,N_{\text{el}};c=1,\dots ,N_{\text{nc}},& (23)\end{array}$

$\begin{array}{ll}{\displaystyle {\int }_{{\Omega }^e}}{M}_q^e\left\{{\dot{r}}^h-g(u^h+{\alpha }^h,r^h)\right\}\mathrm{\text{d}}x=0,e=1,\dots ,N_{\text{el}};q=1,\dots ,N_{\text{q}}& (24)\end{array}$

2.3. Semi-Implicit Temporal Discretization

To integrate (22), (23) and (24) in time, we consider partitioning the time interval into [0, …, t_n, t_n+1, …, T], and approximate the time-dependent coefficients □(t_n) ≈ □_n. For a generic time interval [t_n, t_n+1] we define Δt: = t_n+1−t_n. We further group the expansion coefficients into vectors, and write

$\begin{array}{ll}\begin{array}{c}{u}_n={\left[u_{n,1},\dots ,u_{n,N_{\text{dofs}}}\right]}^T,{\alpha }_n^e={\left[{\alpha }_{n,1}^e,\dots ,{\alpha }_{n,N_{\text{nc}}}^e\right]}^T,\\ r_n^e={\left[r_{n,1}^e,\dots ,r_{n,N_{\text{q}}}^e\right]}^T\end{array}& (25)\end{array}$

Following a semi-implicit (SI) time-integration approach (Whiteley, 2006), time derivatives are replaced by the finite-difference approximation

$\begin{array}{ll}\dot{\square }(t_{n+1})\approx \frac{{\square }_{n+1}-{\square }_n}{\Delta t}.& (26)\end{array}$

Diffusive terms in Equations (22) and (23) are evaluated at t = t_n+1 and the reaction terms are evaluated at t = t_n. Evolution Equation (24) were integrated using an explicit Forward-Euler scheme. As a result, the incremental time update for t = t_n+1 reads: Given $u_n,{\left\{{\alpha }_n^e,r_n^e\right\}}_{e=1,\dots ,N_{\text{el}}}$, find $u_{n+1},{\left\{{\alpha }_{n+1}^e,r_{n+1}^e\right\}}_{e=1,\dots ,N_{\text{el}}}$ such that

$\begin{array}{l}{\displaystyle \sum _{B=1}^{N_{\text{dofs}}}}\left\{{\displaystyle {\int }_{\Omega }}\frac{N_A{N}_B}{\Delta t}+{\displaystyle {\int }_{\Omega }}\nabla N_A·\text{D}\nabla N_B\right\}{u}_{n+1,B}\\ +{\displaystyle \sum _{e=1}^{N_{\text{el}}}}{\displaystyle \sum _{d=1}^{N_{\text{nc}}}}\left\{{\displaystyle {\int }_{\Omega }}\frac{N_A{W}_d^e}{\Delta t}+{\displaystyle {\int }_{\Omega }}\nabla N_A·\mathbf{\text{D}}\nabla W_d^e\right\}{\alpha }_{n+1,d}^e\\ -\left\{{\displaystyle {\int }_{\Omega }}\frac{N_A}{\Delta t}\left\{u_n^h+{\alpha }_n^h\right\}+{\displaystyle {\int }_{\Omega }}{N}_{A}f(u_n^h+{\alpha }_n^h,r_n^h)\right\}=0,\\ A=1,\dots ,N_{\text{dofs}},& (27)\end{array}$

$\begin{array}{l}{\displaystyle \underset{N_{en}}{\overset{b=1}{{{\displaystyle \sum }}^{\text{}}}}}\underset{=:L_{cb}^e}{\underbrace{\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{W_c^e{N}_b^e}{\Delta t}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\nabla W_c^e·\text{D}\nabla N_b^e\}}}{u}_{n+1,b}^e\\ +{\displaystyle \underset{N_{nc}}{\overset{d=1}{{{\displaystyle \sum }}^{\text{}}}}}\underset{{}_{=:K_{{\alpha }_{cd}}^e}}{\underbrace{\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{W_c^e{W}_d^e}{\Delta t}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\nabla W_c^e·\text{D}\nabla W_d^e\}}}{\alpha }_{n+1,d}^e\\ -\underset{=:p_{{\alpha }_c}^e}{\underbrace{\left\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{W_c^e}{\Delta t}\right\{u_n^h+{\alpha }_n^h\}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}{W}_c^{e}f(u_n^h+{\alpha }_n^h,r_n^h)\}}}=0\\ e=1,\dots ,N_{el};c=1,\dots ,N_{nc}& (28)\end{array}$

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }^e}}{M}_q^e\left\{{\displaystyle \sum _{s=1}^{N_q}}{M}_s^e\frac{r_{n+1,s}^e-r_{n,s}^e}{\Delta t}-g(u_n^h+{\alpha }_n^h,r_n^h)\right\}\mathrm{\text{d}}x=0,\\ e=1,\dots ,N_{\text{el}};q=1,\dots ,N_q,& (29)\end{array}$

where $N_b^e:={N_B|}_{{\Omega }^e}$ is the restriction of the basis function to the local element domain, and $u_b^e$ is the corresponding nodal value, where lowercase letters indicate the local degree of freedom b corresponding to its global counterpart B. At this point, we note that Equation (28) can be written in matrix form as

$\begin{array}{l}{L}^e{u}_{n+1}^e+K_{\alpha }^e{\alpha }_{n+1}^e-p^e(u_n^e,{\alpha }_n^e,r_n^e)=0,\end{array}$

for e = 1, …, N_el, from where we define the time update for the element non-conforming coefficient vector as

$\begin{array}{l}{\alpha }_{n+1}^{e,*}(u_{n+1}^e;u_n^e,{\alpha }_n^e,r_n^e):={\left\{K_{\alpha }^e\right\}}^{-1}{p}_{\alpha }^e(u_n^e,{\alpha }_n^e,r_n^e)\\ -{\left\{K_{\alpha }^e\right\}}^{-1}{L}^e{u}_{n+1}^e& (30)\end{array}$

which is computed exclusively using element-level variables, given the element vector $u_{n+1}^e$. To update the gating-variable field, we note from Equation (14) that Equation (29) can be solved point-wise at each quadrature point $x_q$ inside an element, and thus is equivalent to writing

$\begin{array}{l}\frac{r_{q,n+1}^e-r_{q,n}^e}{\Delta t}-g(u_n^h(x_q)+{\alpha }_n^h(x_q),r_{q,n}^e)=0,\\ e=1,\dots ,N_{\text{el}};q=1,\dots ,N_q,\end{array}$

from which the (explicit) time update for the gating variables can be solved at the quadrature-point level as

$\begin{array}{ll}{r}_{q,n+1}^{e,*}(u_n^e,{\alpha }_n^e,r_n^e):=r_{q,n}^e+\Delta tg(u_n^h(x_q)+{\alpha }_n^h(x_q),r_{q,n}^e).& (31)\end{array}$

We now turn to residual Equation (27), and note that it can be constructed by assembling element-level nodal contributions defined by

$\begin{array}{l}{R}_a^{u,e}:={\displaystyle \underset{N_{en}}{\overset{b=1}{{{\displaystyle \sum }}^{\text{}}}}}\underset{:=K_{u_{ab}}^e}{\underbrace{\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{N_a{N}_b}{\Delta t}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\nabla N_a·\text{D}\nabla N_b\}}}{u}_{n+1,b}^e\\ +{\displaystyle \underset{N_{en}}{\overset{b=1}{{{\displaystyle \sum }}^{\text{}}}}}\underset{L_{ad}^{eT}}{\underbrace{\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{N_a{W}_d}{\Delta t}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\nabla N_a·\text{D}\nabla W_d\}}}{\alpha }_{n+1,d}^e\\ -\underset{:=p_{u_a}^e}{\underbrace{\left\{{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}\frac{N_a}{\Delta t}\right\{u_n^h+{\alpha }_n^h\}+{{{\displaystyle \int }}^{\text{}}}_{{\Omega }^e}{N}_{a}f(u_n^h+{\alpha }_n^h,r_n^h)\}}},& (32)\end{array}$

which can also be written in matrix form at the element level as

$\begin{array}{ll}{R}^{u,e}=K_u^e{u}_{n+1}^e+L^{eT}{\alpha }_{n+1}^e-p_u^e(u_n^e,{\alpha }_n^e,r_n^e).& (33)\end{array}$

Substituting update Equation (30) into Equation (33), we obtain an element residual that only depends on $u_{n+1}^e$ that reads

$\begin{array}{l}{R}^{u,e}=\underset{A^e}{\underbrace{(K_u^e-L^{eT}{\left\{K_{\alpha }^e\right\}}^{-1}{L}^e)}}{u}_{n+1}^e\\ +\underset{b_n^e(u_n^e,{\alpha }_n^e,r_n^e)}{\underbrace{L^{eT}{\left\{K_{\alpha }^e\right\}}^{-1}{p}_{\alpha }^e(u_n^e,{\alpha }_n^e,r_n^e)-p_u^e(u_n^e,{\alpha }_n^e,r_n^e)}}& (34)\end{array}$

As a consequence, solving residual Equation (27) is equivalent to solving the matrix linear system

$\begin{array}{ll}A{u}_{n+1}+b_n=0& (35)\end{array}$

where A and $b_n$ are the global matrix and vector assembled from element contributions defined in Equation (34). We note that Equation (35) defines the time update for the global potential vector

$\begin{array}{ll}{u}_{n+1}^{*}(u_n,{\left\{{\alpha }_n^e,r_n^e\right\}}_{e=1,\dots ,N_{\text{el}}}):=-A^{-1}{b}_n& (36)\end{array}$

We remark that matrix A does not depend on the coefficient vectors, and therefore will take the same values for all time steps. Thus, it can be computed on a initialization stage, inverted and stored for later use in updating the potential vector. For the sake of clarity, the steps for the solving the semi-implicit scheme are summarized in Algorithm 1.

ALGORITHM 1

Algorithm 1. Solution algorithm

2.4. The Q1NC Element

We materialize the non-conforming scheme defined in the previous section using incompatible-modes basis functions (Wilson et al., 1973; Taylor et al., 1976), which enhance Q1 elements. We recall that the isoparametric basis functions for Q1 3D (solid) elements are

$\begin{array}{l}{\widehat{N}}_1=\frac{1}{8}(1-{\xi }_1)(1-{\xi }_2)(1-{\xi }_3),{\widehat{N}}_2=\frac{1}{8}(1+{\xi }_1)(1-{\xi }_2)(1-{\xi }_3),\\ {\widehat{N}}_3=\frac{1}{8}(1+{\xi }_1)(1+{\xi }_2)(1-{\xi }_3),{\widehat{N}}_4=\frac{1}{8}(1-{\xi }_1)(1+{\xi }_2)(1-{\xi }_3),\\ {\widehat{N}}_5=\frac{1}{8}(1-{\xi }_1)(1-{\xi }_2)(1+{\xi }_3),{\widehat{N}}_6=\frac{1}{8}(1+{\xi }_1)(1-{\xi }_2)(1+{\xi }_3),\\ {\widehat{N}}_7=\frac{1}{8}(1+{\xi }_1)(1+{\xi }_2)(1+{\xi }_3),{\widehat{N}}_8=\frac{1}{8}(1-{\xi }_1)(1+{\xi }_2)(1+{\xi }_3),\end{array}$

where $({\xi }_1,{\xi }_2,{\xi }_3)\in \widehat{\Omega }:={\left[-1,1\right]}^3$, and

$\begin{array}{l}{N}_a^e={\widehat{N}}_a\dot{}{\widehat{x}}^{-1}\end{array}$

with

$\begin{array}{l}\widehat{x}={\displaystyle \sum _{a=1}^8}{\widehat{N}}_a{x}_a^e,\end{array}$

where $x_a^e$ is the vector with nodal coordinates. Incompatible modes enhance the Q1(Ω^e) element basis by adding basis functions ${\left\{M_c^e\right\}}_{c=1,2,3}$, with $M_c^e={\widehat{M}}_c\dot{}{\widehat{x}}^{-1}$, where

$\begin{array}{l}{\widehat{M}}_1=1-{({\xi }_1)}^2,{\widehat{M}}_2=1-{({\xi }_2)}^2,{\widehat{M}}_3=1-{({\xi }_3)}^2\end{array}$

for $({\xi }_1,{\xi }_2,{\xi }_3)\in \widehat{\Omega }$. Table 1 details the number of DOFs used for the 3D elements considered in this work. Integrals have been approximated using Gaussian quadrature on the standard domain. Table 1 reports the quadrature rules employed in the numerical integration of Q1, Q1NC and Q2 element implementations.

TABLE 1

Table 1. Element DOFs and quadrature rules employed in numerical integration of residuals and tangents.

2.5. The Modified Aliev-Panfilov Model for Transmembrane Ionic Current

All simulations considered the modified Aliev-Panfilov model, which accounts for physiological voltage upstroke slopes and conduction velocities (Aliev and Panfilov, 1996; Hurtado et al., 2016), whose expressions are described below for completeness:

$\begin{array}{ll}f(\varphi ,r)=c_1\varphi (\varphi -\alpha )(1-\varphi )-c_{2}r\varphi & (37)\end{array}$

$\begin{array}{ll}g\left(\varphi ,r\right)=\left(\gamma +\frac{{\mu }_{1}r}{{\mu }_2+\varphi }\right)\left(-r-c_2\varphi \left(\varphi -b-1\right)\right)& (38)\end{array}$

where c₁, c₂, α, γ, μ₁, μ₂ and b are constants, whose values are included in Table 2, and are the same employed by Hurtado and Rojas (2018). To account for a steady-state regime, initial values of the recovery value where set to r = 0.1146.

TABLE 2

Table 2. Parameter values for the modified Aliev-Panfilov model.

3. Results

Finite-element simulations using Q1, Q2, and Q1NC element formulations were implemented for the FI and SI time-integration schemes described in the previous section in an enhanced version of FEAP (Taylor, 2014).

3.1. Plane-Wave Tests on CV and CT

A 3D cardiac rod with a total length of 25 mm was discretized using regular hexahedral elements with a uniform element size, with the exception of elements adjacent to the boundary where the size was at times smaller to fit the geometry. To study the effect of the element size, simulations were carried out with mesh sizes ranging from h = 2mm to h = 0.0156mm. A zero-flux boundary condition was assumed for all boundary surfaces, with exception of the left end of the rod which was stimulated with a normalized external current of 20mV/ms, which corresponds to 28, 000μA/cm³, for 2ms to elicit a plane traveling wave along the direction of the rod. A time-step size of 0.001ms was set for all simulations, which is small when compared to standard cardiac simulations using the selected ionic model (Hurtado et al., 2016). Such small time-step size is chosen to minimize the contribution of the temporal discretization error to the overall numerical error. To compute the CV, we tracked the voltage evolution on x₁ = 18mm and x₂ = 22mm and recorded the activation time, which is defined as the time when the ϕ>0.5 for the first time at a certain point. Then, the CV was calculated as the difference between x₂ and x₁ divided by the difference in the activation time. The results for the CV for different element sizes are shown in Figure 1A. All formulations converged to a CV = 36.9cm/s as the mesh size approached h = 0.0156mm. CV monotonically decreased as mesh size was decreased for Q1 and Q2 formulations. The computational effort of simulations in terms of CT is reported in Figure 1B. We observe that the computational demand of simulations monotonically increases as the mesh size decreases, independently of the element formulation. We do observe, however, that the FI time-integration scheme always results in higher CT than the SI scheme for all element formulations considered.

FIGURE 1

Figure 1. CV tests for plane-waves propagating on a 3D bar for FI and SI schemes on different element formulations. (A) Convergence study of CV as a function of the mesh size h. (B) Computational effort in terms of CT as a function of h.

To facilitate the analysis of the accuracy-efficiency trade-off of the different schemes studied, Figure 2 shows the CT vs. the error in CV for the Q1, Q2, and Q1NC formulations for both the implicit and semi-implicit time updates. Since we seek to minimize two objective functions, namely the CT and the CV error, the Pareto frontier, defined as the set of choices that are Pareto-efficient, is included in each subfigure. The subset of the Pareto-efficient cases that correspond to the Q1NC formulation are {1.2, 1.5}[mm] and {1.0, 1.2, 1.5, 2.0}[mm] for the FI and SI cases, respectively.

FIGURE 2

Figure 2. Accuracy-efficiency analysis: Computation time vs. conduction-velocity error for the different spatial discretization schemes using (A) fully-implicit time integration, and (B) semi-implicit time integration. Dashed gray line displays the Pareto frontier, which encompasses optimal cases. Suboptimal combinations are shown in transparent color.

3.2. Benchmark Simulations on a Cardiac Cuboid

We studied the behavior of the SINCFES using as a second test case the benchmark study on a cardiac cuboid developed by Niederer et al. (2011b), and adapted to the case of the Aliev-Panfilov model by Hurtado et al. (2016). To this end, we consider a cuboid domain with dimensions of 20 × 7 × 3mm with cardiac fibers oriented in the longest axis direction. A subdomain with dimensions 1.5 × 1.5 × 1.5mm located at one of the corners of the cuboid was stimulated with an electrical current density of 50, 000/cm³ for 2ms. The normalized longitudinal and transversal conductivities were 0.0952 and 0.0126mm²/ms, respectively. Figure 3A shows the activation map and isochrones obtained on a plane that contains opposite corners in the diagonal, as defined in Niederer et al. (2011b), for a fine (Baseline) and coarse discretization using Q1 elements, and for the same coarse discretization using Q1NC elements. We note that the Q1NC case with mesh size h = 0.8mm resulted in an activation map and isochrones similar to the baseline case, defined as a Q1 model with mesh size h = 0.1mm. In contrast, the activation map delivered by the Q1 coarse-mesh case with mesh size h = 0.8mm largely differed from the baseline case, delivering a less curved wave-front profile. Figure 3 displays the activation time values along the diagonal of the cuboid for the three cases under study. We observe that the Q1NC case closely follows the baseline case, whereas the Q1 coarse-mesh case resulted in shorter activation times at all locations along the diagonal. As a reference, the CT for the Baseline (Q1 fine), Q1NC and Q1 cases were 122, 341, 344, and 184 s, respectively, which is equivalent to a CT ratio of 665 : 2 : 1.

FIGURE 3

Figure 3. Numerical simulations on cuboid benchmark test (A) Meshes and activation maps, and (B) Activation time profile along the cuboid diagonal.

3.3. Biventricular Human Heart Simulations

To study the potential of the Q1NC-SI formulation in whole-heart cardiac simulations, we modeled the propagation of an action potential on an idealized human biventricular domain stimulated at the atrio-ventricular node. The heart biventricular geometry was generated from two truncated ellipsoids (Streeter and Hanna, 1973), and later discretized using non-regular hexahedral elements. For the baseline case, a size-varying mesh with average characteristic length of 0.48mm was employed. A coarse mesh with average element length of 1.0mm was also considered for two additional cases with Q1 and Q1NC element formulations, see left column of Figure 4 for a representation of the biventricular meshes. All three cases considered the same initial boundary conditions and time step size of 0.001ms. The transmembrane potential distribution at different time instants during ventricular activation is depicted in Figure 4. We clearly observe that, as time elapses, the action-potential wave front of the Q1NC case is very similar to the Baseline case, whereas the Q1 case results in a wave front that propagates faster than the other two models due to the artificially high CV. The last column in Figure 4 shows the activation maps, where we observe that isochrones for the Baseline and Q1NC cases are very similar, and they both differ from the Q1 case. Biventricular simulations were ran in a HPC cluster with 128 GB of RAM memory using 32 processors using the parallel implementation of the code FEAP (Taylor, 2014). The CT for the baseline, the Q1NC and the Q1 simulation were 1805, 452 and 154 min respectively, which is equivalent to a CT ratio of 18 : 3 : 1.

FIGURE 4

Figure 4. Numerical simulations on human biventricular idealized geometries. The Q1NC model displays a propagating wave similar to the baseline case during the ventricular activation sequence, whereas the Q1 model hastens the electrical stimulation ahead of the baseline case.

3.4. Spiral Wave Simulations

To assess the performance of the proposed non-conforming scheme in the simulation of spiral waves, we considered a 50 × 50mm cardiac domain excited by means of an S1–S2 stimulation protocol. To this end, we first applied a surface stimulus (S1) of 12mV/(msmm²) for 2ms on the border defined by x = 0 to create a plane wave. After 280ms, we applied a second stimulus (S2) of 15mV/(msmm³) in the quadrant x < 25, y < 25mm for 5ms, which resulted in the formation of a spiral wave (Costabal et al., 2017). This S1–S2 model was solved using three numerical models: a fine mesh with element size h = 0.1mm using Q1 elements (Baseline), a coarse mesh with element size h = 1mm using Q1 elements (Q1), and a coarse mesh with h = 1mm using the proposed non-conforming element formulation (Q1NC). In all cases, we considered a semi-implicit time update with time-step size Δt = 0.005ms. Figure 5 shows the distribution of the transmembrane potential of the three models under study for several time instants. We note that at early times (t = 110ms) the Q1 case displays a wave front that advanced considerably faster than the baseline and Q1NC cases. At t = 400ms a spiral wave formed in the Baseline and Q1NC cases, whereas for the Q1 case a curved wave front propagated in the outward direction but did not create a spiral. At a later instant (t = 600ms), a spiral was steadily rotating in the Baseline and Q1NC cases, constantly reexciting tissue, whereas in the Q1 case cardiac tissue was found under complete rest, and no electrical activity was observed.

FIGURE 5

Figure 5. Spiral generation simulation in a 2D slab. Due to the higher CV, the Q1 case (coarse mesh) cannot capture the genesis of the spiral wave.

4. Discussion

In this work we have studied the features and advantages of a novel SINCFES in the solution of the monodomain model of cardiac electrophysiology. From plane-wave CV tests we note that the FI and SI schemes yield similar results for the conduction velocity for the time-step size employed, see Figure 1A. This is expected, as the time-step size considered here is small compared to standard values employed in numerical simulations (Krishnamoorthi et al., 2013). Interestingly, we observe that in the case of mesh sizes h < 0.6mm, the Q1, Q2, and Q1NC element formulations delivered very similar results in terms of CV error. For the cases where h>0.6mm, the CV error incurred by the Q1 formulation grows at a much faster rate than the Q2 and Q1NC formulations. An interesting result that deserves further study is the convergence trend of the Q1NC formulation, as it is not monotonically convergent in the whole range of mesh sizes studied, and it reverts the sign of the CV error in a bounded interval of mesh sizes. A similar convergence trend has been reported in the literature for standard FE discretizations, in the context of mass-lumping techniques (Pezzuto et al., 2016), which suggest as future work a more detailed study of the effect of NC spatial discretization schemes on the stiffness and mass matrices that govern the dynamics of the problem. To better analyze the accuracy-efficiency trade-off for each scheme, we constructed CT vs CV-error plots, where the Pareto frontier has been identified. We conclude that the SINCFES delivers Pareto-optimal results for cases with mesh size in the range of {1.0, 1.2, 1.5, 2.0}[mm]. For smaller mesh sizes, traditional Q1 formulations deliver better combinations of CT and CV-error than Q1NC and Q2. It is interesting to note that, in general, Q2 elements are less efficient than the Q1 and Q1NC elements from a Pareto-optimality viewpoint for the whole range of mesh sizes studied. We also note that these conclusions are particular to a plane-wave propagation case, where anisotropy of conductivity and curvature of propagating wavefronts are absent.

We further studied the performance of Q1NC elements using a benchmark problem on a cuboid cardiac domain (Niederer et al., 2011b). Our simulations showed that the Q1NC formulation on a coarse mesh (h = 0.8mm) can result in activation maps that are similar to those obtained on fine meshes using Q1 (h = 0.1mm), adequately capturing the anisotropic conduction of the propagating waves, see Figure 3A. An analysis of the activation-time profile along the cuboid diagonal shows that the Q1NC scheme delivers an accurate conduction velocity, which is comparable to Q1 meshes with mesh sizes that are 8 times smaller, see Figure 3B. This result confirms the ability of Q1NC elements to capture the propagation of electrical waves in anisotropic media with good accuracy at significantly reduced CTs. In contrast, Q1 coarse-mesh simulations resulted in markedly higher conduction velocity profiles, and did poorly in capturing the anisotropic propagation of wavefronts when compared to the Q1NC formulation.

Numerical simulations on a biventricular domain showed that our non-conforming scheme can be effectively used in unstructured meshes of idealized anatomical geometries of the heart, see Figure 4. Similarly to the cardiac rod case, a coarse mesh using Q1NC elements performs much better than a simulation using standard Q1 elements on the same discretization level, as it predicts more accurately the wavefront propagation pattern, when compared to the baseline case. This conclusion is also reached from observing the resulting activation maps, where the spatial distribution and curved shape of isochrones in the Q1NC and baseline are similar, whereas the Q1 formulation delivers an isochrone distribution with lower activation-time values. We note here that this study considered an idealized and smooth geometrical representation of the ventricles of the human heart, useful for numerical verification purposes. It is important to note that such idealized domain does not include important anatomical structures such as the intricate endocardial surface, papillary muscles, and Purkinje network, that are currently included in advanced heart models (Ponnaluri et al., 2016; Sahli Costabal et al., 2016). Future work should focus in understanding how non-conforming formulations can handle such fine-scale anatomical details and structures.

The performance of the SINCFES was studied in the simulation of spiral waves. Remarkably, a very coarse mesh using Q1NC elements is capable to correctly produce, and sustain in time, a spiral wave, whereas a standard Q1 formulation using the same mesh size results in no activation of cardiac tissue. The ability of SINCFES to reproduce spiral wave formation and dynamics is a key result of this work, as it shows that the method is physically more accurate than standard FE formulations for coarse discretizations. This result can be explained by the reduced dependance of the CV on the mesh size, and highlights the potential of the SINCFES in the simulation of cardiac arrhythmias, the main clinical focus of cardiac electrophysiology simulations. While spiral patterns and dynamics obtained with the Q1NC formulation are very similar to the baseline results, a time delay is observed for the former, which resulted in differences in the spatial distribution of the transmembrane voltage, see last column of Figure 5. Such delay, which can ultimately be attributed to differences in the local CV, has also been observed in studies employing very high-order space-time formulations (Coudière and Turpault, 2017), confirming that state-of-the-art simulations of spirals using standard values of mesh size and time step are also affected by this time delay. Despite this persistent numerical error, we believe that the focus of future studies should be in recovering the overall dynamical features of spirals, i.e., spiral tip trajectories (Fenton and Karma, 1998; Gizzi et al., 2013).

We close by noting that while whole-heart simulations reported in the literature predominantly employ tetrahedral discretizations, effective methods for generating patient-specific hexahedral meshes are currently available (Lamata et al., 2011). Further, hexahedral meshes have gained great attention in the context of cardiac simulations, as the numerical performance of hexahedral elements is superior to tetrahedral elements when solving mechanics of the heart, particularly under the assumption of incompressible and quasi-incompressible regimes (Hadjicharalambous et al., 2014). As a conclusion, a natural continuation of this work is the application of non-conforming schemes in the solution of electromechanical models of the heart (Nash and Panfilov, 2004). One important reason for mesh-coarsening FE models of the heart is to reduce the number of DOFs, which in the case of electromechanical cardiac models is much larger than in pure electrophysiological simulations, as displacement, fiber stretch/stress variables, and the non-linearity of tissue constitutive models drastically increase the dimensionality and computational effort needed to solve the governing equations (Göktepe and Kuhl, 2010; Hurtado et al., 2017).

Author Contributions

JJ and DH developed the theoretical framework, numerical schemes and computer algorithms. DH designed the numerical experiments. JJ coded the implementation and ran simulations. JJ and DH wrote the manuscript draft. DH reviewed the final version of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This research was financially supported by the Chilean Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) through grant #1180832. The support of the UCT Three-Way PhD Global Partnership Programme funded by the Dr. Leopold und Carmen Ellinger Stiftung is also appreciated. Authors would like to thank the reviewers for the constructive feedback and comments.

References

Aliev, R. R., and Panfilov, A. V. (1996). A simple two-variable model of cardiac excitation. Chaos Solitons Fractals 7, 293–301. doi: 10.1016/0960-0779(95)00089-5

CrossRef Full Text | Google Scholar

Arthurs, C. J., Bishop, M. J., and Kay, D. (2012). Efficient simulation of cardiac electrical propagation using high order finite elements. J. Comput. Phys. 231, 3946–3962. doi: 10.1016/j.jcp.2012.01.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantwell, C. D., Yakovlev, S., Kirby, R. M., Peters, N. S., and Sherwin, S. J. (2014). High-order spectral/hp element discretisation for reaction-diffusion problems on surfaces: application to cardiac electrophysiology. J. Comput. Phys. 257, 813–829. doi: 10.1016/j.jcp.2013.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Costabal, F. S., Concha, F. A., Hurtado, D., and Kuhl, E. (2017). The importance of mechano-electrical feedback and inertia in cardiac electromechanics. Comput. Methods Appl. Mech. Eng. 320, 352–368. doi: 10.1016/j.cma.2017.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Coudière, Y., and Turpault, R. (2017). Very high order finite volume methods for cardiac electrophysiology. Comput. Math. Appl. 74, 684–700. doi: 10.1016/j.camwa.2017.05.012

CrossRef Full Text | Google Scholar

Dupraz, M., Filippi, S., Gizzi, A., Quarteroni, A., and Ruiz-Baier, R. (2015). Finite element and finite volume-element simulation of pseudo-ECGs and cardiac alternans. Math. Methods Appl. Sci. 38, 1046–1058. doi: 10.1002/mma.3127

CrossRef Full Text | Google Scholar

Fenton, F. H., and Karma, A. (1998). Vortex dynamics in three-dimensional continous myocardium with fiber rotation: Filament instablity and fibrillation. Chaos 8, 20–47. doi: 10.1063/1.166311

CrossRef Full Text | Google Scholar

Gizzi, A., Cherry, E. M., Gilmour, R. F., Luther, S., Filippi, S., and Fenton, F. H. (2013). Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue. Front. Physiol. 4:71. doi: 10.3389/fphys.2013.00071

PubMed Abstract | CrossRef Full Text | Google Scholar

Göktepe, S., and Kuhl, E. (2009). Computational modeling of cardiac electrophysiology: a novel finite element approach. Int. J. Num. Methods Eng. 79, 156–178. doi: 10.1002/nme.2571

CrossRef Full Text | Google Scholar

Göktepe, S., and Kuhl, E. (2010). Electromechanics of the heart: a unified approach to the strongly coupled excitation-contraction problem. Comput. Mech. 45, 227–243. doi: 10.1007/s00466-009-0434-z

CrossRef Full Text | Google Scholar

Hadjicharalambous, M., Lee, J., Smith, N. P., and Nordsletten, D. A. (2014). A displacement-based finite element formulation for incompressible and nearly-incompressible cardiac mechanics. Comput. Methods Appl. Mech. Eng. 274, 213–236. doi: 10.1016/j.cma.2014.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurtado, D. E., Castro, S., and Gizzi, A. (2016). Computational modeling of non-linear diffusion in cardiac electrophysiology: a novel porous-medium approach. Comput. Methods Appl. Mech. Eng. 300, 70–83. doi: 10.1016/j.cma.2015.11.014

CrossRef Full Text | Google Scholar

Hurtado, D. E., Castro, S., and Madrid, P. (2017). Uncertainty quantification of two models of cardiac electromechanics. Int. J. Num. Methods Biomed. Eng. 33:e2984. doi: 10.1002/cnm.2894

CrossRef Full Text | Google Scholar

Hurtado, D. E., and Henao, D. (2014). Gradient flows and variational principles for cardiac electrophysiology: Toward efficient and robust numerical simulations of the electrical activity of the heart. Comput. Methods Appl. Mech. Eng. 273, 238–254. doi: 10.1016/j.cma.2014.02.002

CrossRef Full Text | Google Scholar

Hurtado, D. E., and Kuhl, E. (2014). Computational modelling of electrocardiograms: repolarisation and T-wave polarity in the human heart. Comput. Methods Biomech. Biomed. Eng. 17, 986–996. doi: 10.1080/10255842.2012.729582

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurtado, D. E., and Rojas, G. (2018). Non-conforming finite-element formulation for cardiac electrophysiology: an effective approach to reduce the computation time of heart simulations without compromising accuracy. Comput. Mech. 61, 485–497. doi: 10.1007/s00466-017-1473-5

CrossRef Full Text | Google Scholar

Krishnamoorthi, S., Sarkar, M., and Klug, W. S. (2013). Numerical quadrature and operator splitting in finite element methods for cardiac electrophysiology. Int. J. Num. Methods Biomed. Eng. 29, 1243–1266. doi: 10.1002/cnm.2573

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamata, P., Niederer, S., Nordsletten, D., Barber, D. C., Roy, I., Hose, D. R., et al. (2011). An accurate, fast and robust method to generate patient-specific cubic Hermite meshes. Med. Image Anal. 15, 801–813. doi: 10.1016/j.media.2011.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Nash, M. P., and Panfilov, A. V. (2004). Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Prog. Biophys. Mol. Biol. 85, 501–522. doi: 10.1016/j.pbiomolbio.2004.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S., Mitchell, L., Smith, N., and Plank, G. (2011a). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol. 2:14. doi: 10.3389/fphys.2011.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Kerfoot, E., Benson, A. P., Bernabeu, M. O., Bernus, O., Bradley, C., et al. (2011b). Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 369, 4331–4351. doi: 10.1098/rsta.2011.0139

PubMed Abstract | CrossRef Full Text | Google Scholar

Pathmanathan, P., Bernabeu, M. O., Bordas, R., Cooper, J., Garny, A., Pitt-Francis, J. M., et al. (2010). A numerical guide to the solution of the bidomain equations of cardiac electrophysiology. Prog. Biophys. Mol. Biol. 102, 136–155. doi: 10.1016/j.pbiomolbio.2010.05.006

CrossRef Full Text | Google Scholar

Pezzuto, S., Hake, J., and Sundnes, J. (2016). Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology. Int. J. Num. Methods Biomed. Eng. 32:e02762. doi: 10.1002/cnm.2762

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponnaluri, A. V., Perotti, L. E., Liu, M., Qu, Z., Weiss, J. N., Ennis, D. B., et al. (2016). Electrophysiology of heart failure using a rabbit model: from the failing myocyte to ventricular fibrillation. PLoS Comput. Biol. 12:e1004968. doi: 10.1371/journal.pcbi.1004968

PubMed Abstract | CrossRef Full Text | Google Scholar

Pullan, A. J., Cheng, L. K., and Buist, M. L. (2005). Mathematically Modelling the Electrical Activity of the Heart: From Cell to Body Surface and Back Again. Singapore: World Scientific.

Google Scholar

Sahli Costabal, F., Hurtado, D. E., and Kuhl, E. (2016). Generating Purkinje networks in the human heart. J. Biomech. 49, 2455–2465. doi: 10.1016/j.jbiomech.2015.12.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahli Costabal, F., Yao, J., and Kuhl, E. (2018). Predicting drug-induced arrhythmias by multiscale modeling. Int. J. Num. Methods Biomed. Eng. 34:e2964. doi: 10.1002/cnm.2964

PubMed Abstract | CrossRef Full Text | Google Scholar

Streeter, D. D., and Hanna, W. T. (1973). Engineering mechanics for successive states in canine left ventricular myocardium: I. Cavity and wall geometry. Circ. Res. 33, 639–655. doi: 10.1161/01.RES.33.6.639

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., Cai, X., Nielsen, B. F., Mardal, K. A., and Tveito, A. (2006). Computing the Electrical Activity in the Heart. Berlin; Heidelberg: Springer-Verlag.

Google Scholar

Taylor, R. L. (2014). FEAP-Finite Element Analysis Program. Berkeley, CA.

Google Scholar

Taylor, R. L., Beresford, P. J., and Wilson, E. L. (1976). A Nonconforming Element for Stress Analysis. Int. J. Num. Methods Eng. 10, 1211–1219. doi: 10.1002/nme.1620100602

CrossRef Full Text | Google Scholar

Trayanova, N. A., Boyle, P. M., and Nikolov, P. P. (2018). Personalized imaging and modeling strategies for arrhythmia prevention and therapy. Curr. Opin. Biomed. Eng. 5, 21–28. doi: 10.1016/j.cobme.2017.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Vadakkumpadan, F., Rantner, L. J., Tice, B., Boyle, P., Prassl, A. J., Vigmond, E., et al. (2009). Image-based models of cardiac structure with applications in arrhythmia and defibrillation studies. J. Electrocardiol. 42, 157.e1–157.e10. doi: 10.1016/j.jelectrocard.2008.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Vazquez, M., Aris, R., Houzeaux, G., Aubry, R., Villar, P., Garcia-Barnes, J., et al. (2011). A massively parallel computational electrophysiology model of the heart. Int. J. Num. Methods Biomed. Eng. 27, 1911–1929. doi: 10.1002/cnm.1443

CrossRef Full Text | Google Scholar

Whiteley, J. P. (2006). An efficient numerical technique for the solution of the monodomain and bidomain equations. IEEE Trans. Biomed. Eng. 53, 2139–2147. doi: 10.1109/TBME.2006.879425

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, E. L., Taylor, R. L., Doherty, W. P., and Ghaboussi, J. (1973). “Incompatible displacement models,” in Numerical and Computer Methods in Structural Mechanics, eds S. J. Fenves, N. Perrone, A. R. Robinson, and W. C. Schnobrich (Orlando, FL: Academic Press), 668.

Google Scholar

Ying, W., Rose, D. J., and Henriquez, C. S. (2008). Efficient fully implicit time integration methods for modeling cardiac dynamics. IEEE Trans. Biomed. Eng. 55, 2701–2711. doi: 10.1109/TBME.2008.925673

PubMed Abstract | CrossRef Full Text | Google Scholar