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.


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 wellknown 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 highperformance 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.1 mm (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 timestepping 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 wellsuited 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.

Monodomain Model of Cardiac Electrophysiology
Let ∈ R 3 be the heart domain where electrical impulses travel during the time interval [0, T], and V m : × [0, T] → R be the transmembrane potential. A local statement of current balance yields the monodomain equation (Pullan et al., 2005) (1) 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 : × (0, T] → R m is a vector field of state variables that include gating variables and ion concentrations. For convenience, we consider the normalized transmembrane potential field where V p and V r are the peak and resting voltages, respectively. Based on this normalization, we obtain the non-dimensional monodomain equation, where D = 1 A m C m σ is the normalized conductivity tensor, and f (φ, r) = − I ion (V m (φ),r) C m (V p −V r ) is the normalized ionic current. The time evolution of state variables is governed by kinetic equations of the form The expressions for f (φ, r) and g(φ, 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, respectively, as well as initial conditions To state the weak form of the cardiac electrophysiology problem, we consider trial spaces S φ , S r and test spaces V φ , V r defined as 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 (φ, r) ∈ S φ × S r such that

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 h = ∪ N el e=1 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 to solve the following FE problem (Göktepe and Kuhl, 2009;Hurtado and Kuhl, 2014): 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ˆ = [−1, 1] n and mapping to a hexahedral element. We expand an element 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 where {u A (t)} A=1,N dofs correspond to the nodal values of the transmembrane potential field, and u BC ∈ S φ is a function that satisfies the boundary conditions (4), i.e., u BC =φ 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 r h , we write where M e q is a characteristic function defined by and e,q ⊂ e is the subdomain containing the q−quadrature point x q , and is such that Nq q=1 e,q = e and e,q ∩ e,q ′ = ∅ whenever q = q ′ . Analogously, we expand an element r h ∈ S r h as where r e q :(0, T] → R m represents the time evolution of the state variables at the q-quadrature point. In this work, we consider a non-conforming spatialdiscretization scheme for the monodomain equations (Hurtado and Rojas, 2018). To this end, we rewrite the residuals as 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 at the element level by adding polynomial terms not included in Q k ( e ). To this end, we consider the non-conforming space where m ∈ Z + and P k+m ( e ) is the space of polynomial functions of degree k + m defined on the standard domainˆ . We then consider enhanced test functions ν h which we expand as where β e c ∈ R are coefficients, W e c are non-conforming element basis functions, and it holds that W e c = 0, x / ∈ e . Analogously, we enhance S φ h with the time-dependent non-conforming space F φ h , and expand the enhanced trial functions as where and α e d :(0, T] → R is a time-dependent coefficient that scales the non-conforming basis functions W e d . Substitution of approximations Equations (13 15, 18, and 19) into the residuals Equations (16) and (17) yields the following semi-discrete problem:

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 Following a semi-implicit (SI) time-integration approach (Whiteley, 2006), time derivatives are replaced by the finite-difference approximatioṅ 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 , {α e n , r e n } e=1,...,N el , find u n+1 , {α e n+1 , r e n+1 } e=1,...,N el such that where N e b : = N B e is the restriction of the basis function to the local element domain, and u e b 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 L e u e n+1 + K e α α e n+1 − p e (u e n , α e n , r e n ) = 0, for e = 1, . . . , N el , from where we define the time update for the element non-conforming coefficient vector as which is computed exclusively using element-level variables, given the element vector u e n+1 . 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 , r e q,n ) = 0, e = 1, . . . , N el ; q = 1, ..., N q , from which the (explicit) time update for the gating variables can be solved at the quadrature-point level as r e, * q,n+1 (u e n , α e n , r e n ) : = r e q,n + t g(u h n (x q ) + α h n (x q ), r e q,n ). (31) We now turn to residual Equation (27) which can also be written in matrix form at the element level as R u,e = K e u u e n+1 + L eT α e n+1 − p e u (u e n , α e n , r e n ).
Substituting update Equation (30) into Equation (33), we obtain an element residual that only depends on u e n+1 that reads + L eT K e α −1 p e α (u e n , α e n , r e n ) − p e u (u e n , α e n , r e n ) b e n (u e n ,α e n ,r e n ) As a consequence, solving residual Equation (27) is equivalent to solving the matrix linear system where A and b n are the global matrix and vector assembled from element contributions defined in Equation (34). We note that Equation (35) 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.

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 arê Compute K e α , K e u and L e (Equations (28) and (32)) and store Compute A e (Equation (34)) and assemble contribution to A end Compute A −1 and store / * time integration loop * / for n = 0 to N steps do for e = 1 to N el do Compute b e (u e n , α e n , r e n ) (Equation (34)) and assemble contribution to b n end Update u n+1 = u * n+1 (u n , {α e n , r e n } e=1,...,N el ) = −A −1 b n for e = 1 to N el do Update α e n+1 = α e, * n+1 (u e n+1 ; u e n , α e n , r e n ) (see Equation Update r e q,n+1 = r e, * q,n+1 (u e n , α e n , r e n ) (see Equation 31) end end end for (ξ 1 , ξ 2 , ξ 3 ) ∈ˆ . 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.

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 where c 1 , c 2 , α, γ , µ 1 , µ 2 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.

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

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 = 2 mm to h = 0.0156 mm. 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 20 mV/ms, which corresponds to 28, 000 µA/cm 3 , for 2 ms to elicit a plane traveling wave along the direction of the rod. A time-step size of 0.001 ms was set for all simulations, which is small when compared to standard cardiac simulations using the selected ionic model . 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 1 = 18 mm and x 2 = 22 mm 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 2 and x 1 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.9 cm/s as the mesh size approached h = 0.0156 mm. 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.
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.

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 × 3 mm with cardiac fibers oriented in the longest axis direction. A subdomain with dimensions 1.5 × 1.5 × 1.5 mm located at one of the corners of the cuboid was stimulated with an electrical current density of 50, 000/cm 3 for 2 ms. The normalized longitudinal and transversal conductivities were 0.0952 and 0.0126 mm 2 /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.8 mm resulted in an activation map and isochrones similar to the baseline case, defined as a Q1 model with mesh size h = 0.1 mm. In contrast, the activation map delivered by the Q1 coarse-mesh case with mesh size h = 0.8 mm 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.

Biventricular Human Heart Simulations
To study the potential of the Q1NC-SI formulation in wholeheart 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.48 mm was employed. A coarse mesh with average element length of 1.0 mm 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.001 ms. 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.

Spiral Wave Simulations
To assess the performance of the proposed non-conforming scheme in the simulation of spiral waves, we considered a 50 × 50 mm cardiac domain excited by means of an S1-S2 stimulation protocol. To this end, we first applied a surface stimulus (S1) of 12 mV/(ms mm 2 ) for 2 ms on the border defined by x = 0 to create a plane wave. After 280 ms, we applied a second stimulus (S2) of 15 mV/(ms mm 3 ) in the quadrant x < 25, y < 25 mm for 5 ms, 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.1 mm using Q1 elements (Baseline), a coarse mesh with element size h = 1 mm using Q1 elements (Q1), and a coarse mesh with h = 1 mm using the proposed non-conforming element formulation (Q1NC). In all cases, we considered a semi-implicit time update with timestep size t = 0.005 ms. 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 = 110 ms) the Q1 case displays a wave front that advanced considerably faster than the baseline and Q1NC cases. At t = 400 ms 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 = 600 ms), 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.

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.6 mm, the Q1, Q2, and Q1NC element formulations delivered very similar results in terms of CV error. For the cases where h > 0.6 mm, 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 CVerror 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.8 mm) can result in activation maps that are similar to those obtained on fine meshes using Q1 (h = 0.1 mm) , 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 patientspecific 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 nonlinearity 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.