A numerical study on the effects of spatial and temporal discretization in cardiac electrophysiology

Millions of degrees of freedom are often required to accurately represent the electrophysiology of the myocardium due to the presence of discretization effects. This study seeks to explore the influence of temporal and spatial discretization on the simulation of cardiac electrophysiology in conjunction with changes in modeling choices. Several finite element analyses are performed to examine how discretization affects solution time, conduction velocity and electrical excitation. Discretization effects are considered along with changes in the electrophysiology model and solution approach. Two action potential models are considered: the Aliev‐Panfilov model and the ten Tusscher‐Noble‐Noble‐Panfilov model. The solution approaches consist of two time integration schemes and different treatments for solving the local system of ordinary differential equations. The efficiency and stability of the calculation approaches are demonstrated to be dependent on the action potential model. The dependency of the conduction velocity on the element size and time step is shown to be different for changes in material parameters. Finally, the discrepancies between the wave propagation in coarse and fine meshes are analyzed based on the temporal evolution of the transmembrane potential at a node and its neighboring Gauss points. Insight obtained from this study can be used to suggest new methods to improve the efficiency of simulations in cardiac electrophysiology.


| INTRODUCTION
Computational models of electrophysiology provide an excellent tool to better understand the inner workings of the heart. They offer a method to complement the vast experimental research that has been devoted to studying the heart in health and pathology. These models quantitatively represent the mechanisms responsible for action potential generation and propagation throughout cardiac tissue in ways that may be challenging to measure experimentally. These mechanisms are often modeled through a set of ordinary differential equations (ODEs) which express the electrical behavior at the cell level and are coupled to a set of partial differential equations (PDEs) to describe the wave propagation in tissue. 1 The ODEs encompassing the cell electrical behavior can be divided into ionic and phenomenological models. Ionic models [2][3][4] keep track of the kinetics of main channels, pumps, exchangers and concentrations that determine the current across the cell membrane. Phenomenological models [5][6][7] are simplified versions of the ionic models which incorporate less biophysical details while maintaining the same general dynamics. 8 The PDEs can take on the form of the bidomain or monodomain equations, whereby the bidomain equations assume the tissue is composed of an intracellular and extracellular space. The monodomain model is a simplification of the bidomain model, where it is assumed that the intracellular and extracellular conductivities are proportional. 8 The systems of differential equations are often solved through the finite difference method, finite element (FE) method and finite volume method. 9 Within the FE framework different treatments are available in the literature regarding storing and interpolating the cell electrophysiology variables. 10 In the case where the variables are stored at the nodes, the options of ionic current interpolation and state variable interpolation are present. 1 On the other hand, for an internal variable approach, the transmembrane potential is considered as a global variable while the ionic state variables are stored at the Gauss points, as in the approach by Göktepe and Kuhl. 11 A more comprehensive description and comparison of these approaches is available in the work of Krishnamoorthi et al. 10 In this manuscript, we will focus on the FE method utilizing the monodomain equations and an internal variable approach.
One challenging aspect of these models is to bridge the massive spatio-temporal scales present in cardiac electrophysiology. The scale ranges from the ionic pump kinetics on cell membranes to patterns observed in electrocardiograms over many heartbeats. Simulating these behaviors numerically leads to the requirement of high resolutions to capture the fast upstroke in time, steep depolarization wavefront in space, microscopic structure and heterogeneities. 12,13 Often times it is necessary to perform simulations using many millions of degrees of freedom for a typical heart geometry. 12,14,15 Furthermore, the discrete propagation velocity of the wavefront can be far from reality if a low resolution mesh is used. Many papers have documented the effect of discretization on conduction velocity (CV) through mesh sensitivity studies, see  to name a few. The importance of mesh sensitivity is highlighted in a benchmark study for verifying cardiac electrophysiology simulators. 19 Since computational demand can be very high, errors of up to 10% in the CV have been considered acceptable. 20 To tackle the computational burden of cardiac electrophysiology, many different methods have been applied. Some researchers have used adaptivity of the discretization in space, 21,22 in time 23 or both. 24 Others have considered p-adaptivity 13,25 which allows an independent change in the order of the shape functions in the elements. Some authors have resorted to modified element formulations, such as non-conforming elements 26 or so-called macro finite elements. 27 The application of isogeometric analysis has also been studied. 28 Further approaches have aimed at improving the solvers of large systems of equations. 9 Finally, operator splitting, where the fast reaction part is separated from the slow diffusion part, 10 has also been a popular method. A review of numerical techniques for solving the bidomain equations is presented by Pathmanathan et al. 29 Regarding different calculation approaches, the time integration scheme and the treatment of the solution of the local system of ODEs can be considered. Different time integration schemes (such as the forward Euler, backward Euler, semi-implicit, Crank-Nicholson methods) have been applied. 30 A comparison of different time integration schemes is presented by Roy et al. 31 Of note also is the use of the Rush-Larsen scheme present in action potential models which makes them easier and enhances their stability. 32,33 Recently, a gating-enhanced semi-implicit approach which treats stiffer variables implicitly and less stiff variables explicitly was proposed by Green and Spiteri. 34 The semi-implicit schemes appear to be favored due to the speed-up that they deliver. 35,36 While an explicit time integration scheme is generally easier to implement, the stability may be severely hindered. In the implicit case, however, the stability is enhanced but a Newton-Raphson iteration is usually required. Moreover, the stability of the implicit time integration is dependent on the nonlinearity of the electrophysiology model and on the chosen model parameters. 37 Taking into account different approaches for solving the Newton-Raphson iterations, the partitioning and method for obtaining the derivatives can be considered. As for the local system of ODEs, one can solve for all variables simultaneously (monolithic approach) or solve for each variable sequentially (staggered approach). Finally, the derivatives for the local and global Newton-Raphson iterations can be calculated exactly (which yields the normal Newton-Raphson iteration) or approximately (which yields a Quasi-Newton approach). 38 We will refer to derivatives in their exact form as analytical tangents and the numerical approximations of the derivatives using a secant as numerical tangents.
Although some of the mentioned methods have been relatively successful in aiding the computational effort, computational load continues to be an important issue in the field of computational electrophysiology. Therefore, the need is still present to further understand discretization effects to develop new and improved methods to reduce the computational load. Discretization lies at the heart of an FE model, since it influences the behavior of the solution algorithm, material response and geometric description. In this manner, the spatio-temporal discretization can underly the stability, accuracy and computational load of a simulation. The goal of this article is to present how discretization influences some of the relevant aspects of cardiac electrophysiology which include: computational times, CV and excitation characteristics. We explore these effects in conjunction with changes in action potential models, solution algorithms, mesh types and model parameters. The solution algorithms depend on the choice of the time integration scheme (implicit or semi-implicit), partitioning approach for the local system of ODEs (monolithic or staggered) and how the derivatives are obtained (analytical or numerical tangent).
This work follows in accordance with other studies that observe the effects of discretization along with changes in the numerical schemes. 1,10,39 Novel contributions of this work include the presentation of the solution time and stability of the various calculation approaches used to solve the action potential model equations. Of these approaches, the staggered method for the local system of ODEs and the use of numerical tangents are new techniques. Further, the impact that material parameters have on the CV convergence is illustrated systematically. Lastly, the excitation in coarse and fine meshes is analyzed based on the spatio-temporal depiction of the transmembrane potential, ionic current and spatial gradient.
This manuscript is organized in the following sections. Section 2 introduces the methods and equations related to FE monodomain cardiac electrophysiology including an introduction of two action potential models and the different calculation approaches considered. The calculation approaches include the time integration schemes, partitioning methods of the local ODE system and the technique for obtaining derivatives. We present numerical experiments in Section 3 to illustrate how discretization affects the solution time and results. Section 4 offers a discussion on the results presented in Section 3. Finally, we provide concluding remarks in Section 5.

| METHODS
In this section we introduce the main equations along with different calculation approaches to represent the electrophysiology of the myocardium. We begin with the PDE defining the propagation of the transmembrane potential in time and space and proceed with the introduction of both a phenomenological and an ionic model to represent the source term. To keep this section brief, the main equations of the ionic model and the FE approximation are omitted in this work, but can be found in the study by Wong et al. 23 The various calculation approaches (time integration schemes, local Newton-Raphson iteration options and differentiation techniques) are also described.

| Governing differential equations of electrophysiology
The monodomain equations can be represented in the following form 19 The first equation is written in terms of the surface-to-volume ratio of cells χ, the specific membrane capacitance C m , the membrane voltage or transmembrane potential Φ, the transmembrane current density I ion and the conductivity tensor σ. In the case of isotropic conductivity, the conductivity tensor is defined through σ = σI, with I being the second order identity tensor. The scalar conductivity value σ can be expressed in terms of the intracellular σ i and extracellular σ e conductivities, respectively, through σ = σ i σ e (σ i + σ e ) −1 . The function f for the transmembrane current density has a set of internal variables y as an argument, which in turn are defined by a system of kinetic ODEs g. The form that f and g take is dependent on the action potential model. Additionally, the notation _ Φ = ∂Φ=∂t is used, while rΦ is defined as the spatial gradient of Φ. Further, the divergence operator is represented through div{−}.
Consolidating Equation (1) 1 , the PDE can be represented as in terms of the electrical flux vector related to transmembrane potential diffusion q = 1 χC m σÁrΦ and a nonlinear current source term f Φ = − I ion C m . The electrical flux can also be expressed as follows where D is the conductivity tensor defined as with isotropic d iso and anisotropic d ani conductivity material parameters. Additionally, f represents the vector aligned with the fiber direction. We define the source term f Φ , internal variables y and kinetic ODEs g either through the Aliev-Panfilov (A-P) model 5 or through the ten Tusscher-Noble-Noble-Panfilov (TNNP) model. 40 Furthermore, Equation (2) is complemented by boundary and initial conditions. For a body ℬ of interest, the boundary ∂ℬ can be decomposed into ∂ℬ Φ and ∂ℬ q based on the boundary conditions. For Dirichlet boundary conditions, Φ = Φ is applied on ∂ℬ Φ . For Neumann boundary conditions, rΦÁn = h is applied on ∂ℬ q . Here, n represents the unit outward normal to ∂ℬ. We applied homogeneous Neumann boundary conditions rΦÁn = h = 0 along with initial conditions (Φ = Φ 0 and y = y 0 in ℬ at t = 0) representing the resting state of the transmembrane potential.

| Aliev-Panfilov model
In the A-P model, the source term f Φ is described through the approach proposed by Göktepe and Kuhl. 11 Whereby, the transmembrane potential is treated as a global variable and the recovery variable is treated as a local variable. Emulating the strategy in the literature for phenomenological cardiac electrophysiology, the transmembrane potential Φ and the time t are expressed through their respective dimensionless counterparts ϕ and τ as with the scaling factors β ϕ , β t and the potential difference δ ϕ . The current source term is based on the Fitz-Hugh Nagumo type excitation equation modified by Aliev and Panfilov 5 to describe the action potential of myocardial cells with an externally applied electrical stimulus I and the material parameters c, α, γ, μ 1 , μ 2 and b. The term F ϕ is the global source term for the fast acting depolarization, while F r represents the source term for the local repolarization behavior responsible for the evolution of the local recovery variable r. This results in the set of internal variables y = {r}. The physical part of Equation (6) 1 is given by with its sensitivity with respect to Φ being

| Ten Tusscher-Noble-Noble-Panfilov model
The TNNP model 40 describes the evolution of the transmembrane potential based on 15 ionic currents, 4 ion concentrations and 13 gating variables. We follow the implementation and notation presented by Wong et al 23 which fits very well in an FE framework. The TNNP model serves to provide the source term f Φ expressed in terms of the ionic currents I crt which are in turn functions of the transmembrane potential Φ, gating variables g gate and ion concentration variables c ion = − I Na + I bNa + I NaK + I NaCa + I K1 + I Kr + I Ks + I pK + I t0 + I CaL + I bCa + I pCa + I stim À Á : Further, the sensitivity of the source term f Φ with respect to Φ is obtained through The currents are expressed as functions and organized as a set of variables I crt through I crt = I crt Φ,g gate , c ion , I crt = I Na , I bNa , I NaK , I NaCa , I K1 , I Kr , I Ks , I pK , I t0 , I CaL ,I bCa , I pCa , I leak , I up , I rel À Á : The set I crt includes the fast sodium current I Na , background sodium current I bNa , sodium potassium pump current I NaK , sodium calcium exchanger current I NaCa , inward rectifier current I K1 , rapid delayed rectifier current I Kr , slow delayed rectifier current I Ks , plateau potassium current I pK , transient outward current I t0 , L-type calcium current I CaL , background calcium current I bCa , plateau calcium current I pCa , leakage current I leak , sarcoplasmic reticulum (SR) uptake current I up , and the SR release current I rel . A stimulus current I stim is introduced to cause depolarizations but it is not included in I crt because it does not depend on any variables.
Following the suggestion by Wong et al, 23 we consider the gate variables in two different parts: one part contains the 10 gating variables for which their evolution in time depends only on themselves and the transmembrane potential Φ, whereas the second part contains three gating variables that depend additionally on the concentrations c ion . These gating variables are expressed through a set of ODEs and can be organized in sets of variables as follows The ion concentrations c ion consisting of the free intracellular ion concentrations of sodium (Na + ), potassium (K + ) and calcium (Ca 2+ ) and the free ion concentration of calcium in the SR (Ca 2 + sr ) are expressed through a set of ODEs and can be organized as a set of variables where f c ion is a weighted sum of the currents related to a given ion. The total set of internal variables is thus y = {g gate , c ion }. The evolution equations for the concentrations are given as where C is the membrane capacitance per unit area, V is the cytoplasmic volume, V sr is the volume of the SR and F is the Faraday constant. Furthermore, the variables γ Ca and γ sr Ca are weighting factors to account for the total calcium and SR calcium concentrations, respectively.
Due to its strong effect on depolarization, we show the equation for the fast sodium current where C max Na is the maximal conductance of the fast sodium current, g m is the sodium activation gate, g h is the fast sodium inactivation gate, g j is the slow sodium inactivation gate and Φ Na is the concentration dependent Nernst or reversal potential. We introduce a slight change for the sodium activation gate g m by which the value of −56.86 in g ∞ m is replaced by a new material parameter Ξ. This parameter is introduced because it alters the excitation threshold, i. e., the value of the transmembrane potential at which a depolarization will be induced in an un-excited myocyte. Such changes in the kinetics of the fast sodium current have been observed due to variations in the intra-or extracellular concentrations of Ca 2+ and Mg 2+ or in the case of genetic modifications. 41 Nonetheless, we will focus more on the numerical effects accompanying the change in parameter Ξ, which can help to understand how the activation threshold affects the dependence of the CV on the spatial and temporal discretization. As a result, we obtain For the sake of compactness, the rest of the equations for the ionic currents I crt and the equations related to the gate variables g gate are omitted in this work. For a more comprehensive presentation of the model, we refer to the work of ten Tusscher et al 40 and Wong et al. 23

| Calculation approaches
In the following, we explain the different methods that we apply when solving the electrophysiology problem in the FE setting. Two different time integration schemes are considered. A monolithic and a staggered approach when solving the local system of ODEs are outlined. Lastly, the derivatives for the Newton-Raphson iterations are treated in an exact (analytical) or approximated (numerical) manner.

| Time integration schemes: Implicit and semi-implicit
We consider two different time integration schemes for variables related to the ionic current f Φ based on a finite difference approximation in time. The first being an implicit backward Euler scheme and the second an explicit forward Euler scheme. Since the diffusion term in Equation (2) is treated in an implicit manner, using an explicit scheme for the ionic current variables leads to a semi-implicit solution overall. Note that in the semi-implicit case, a local Newton-Raphson iteration is not required because the local variables can be calculated explicitly. In the following presentation of the time integration schemes, we omit the subscript n + 1 for clarity.
Regarding the A-P model, the implicit scheme is defined by while the explicit scheme is obtained by Note that another option can be considered where r is explicit but ϕ is implicit, namely For the TNNP model, the gating variables g gate and concentration variables c ion are integrated in the implicit scheme via To obtain the source term f Φ for a given value of Φ, we need to solve this system of 17 ODEs related to the concentration and gating variables at the material level. The explicit scheme for the TNNP model is obtained through

| Local Newton-Raphson iteration: Monolithic and staggered approach
For applying the Newton-Raphson method on systems of more than one variable, we can decide whether we should compute all variables simultaneously (monolithic) or sequentially (staggered). For the A-P model, no difference exists between a staggered or a monolithic approach because the model has only one local ODE. On the other hand, the TNNP model has a system of nonlinear equations that needs to be solved at each time step. For the monolithic case in the TNNP model, we utilize the algorithm suggested by Wong et al. 23 This algorithm is a fully implicit local update of the equations. Basically, the gate variables g gate and the ionic currents I crt are calculated before entering a local Newton-Raphson iteration responsible for calculating the concentrations. In the local Newton-Raphson iteration, the residuals of the ion concentrations and the linearization of the residuals with respect to each concentration are used to update the concentrations. Then, the gating variables g II gate and ionic currents I crt are updated. This is performed until the Euclidean norm of the residual vector is below a certain tolerance value. At the end of this iteration, the source term f Φ can be calculated based on the ionic currents I crt . The residual vector is composed of the residual form of the evolution of the concentrations presented in Equation (14) as follows The derivatives of the residuals with respect to each of the ion concentrations can be expressed in a non-symmetric matrix form Note the use of the notation ∂ x y = ∂y/∂x. The iterative solution depends on the following linearization Solving for the change in the concentrations Δc yields such that the concentrations can be updated in the following way In the staggered solution, we perform four Newton-Raphson iterations sequentially. Each Newton-Raphson iteration solves for a single concentration and is exited when the residual of that concentration is below a tolerance value. This means that the linearizations and updates take on the following form

| Derivatives for the Newton-Raphson iterations: Analytical and numerical tangent
The derivatives required for the Newton-Raphson iterations, when the implicit time integration scheme is used, can be calculated analytically or numerically. We will refer to derivatives in their exact form as analytical tangents and the numerical approximations of the derivatives using a secant as numerical tangents. Analytical (algorithmic) tangents are determined via symbolic computation (i.e., by applying the laws of calculus for differentiation). To calculate the numerical tangent of a function of interest with respect to a differentiating variable, we first compute the value of the function using the present values of the global and state variables. Then, a second value of the function is calculated using the present global and state variables, except that a perturbation is introduced for the differentiating variable. Lastly, the ratio between the differences in the function and the differentiating variable provides the numerical tangent. For example, for the derivative of the residual for the sodium concentration with respect to the calcium concentration where R Ã c Na and c Ã Ca represent the values where the perturbation c Ca,p is introduced. When applying the numerical tangent approach, we use this approximation for the derivatives with respect to the concentrations and for the sensitivity of the source term with respect to Φ. As a simplification, the same perturbation value c ion,p is used for all of the concentrations, i.e. c ion,p = c K,p = c Na,p = c Ca,p = c sr Ca,p .

| NUMERICAL EXPERIMENTS
In this section, we present three sets of numerical experiments which explore the effects of discretization along with changes in the calculation approaches and material parameters. The first set is a comparison of the calculation approaches in terms of computational time, stability and CV convergence. The convergence behavior of the CV for changes in the material parameters is explored in the second set. In the last part, we try to explain the influence of the spatial discretization on the local excitation and its implications on wave propagation. In all numerical experiments, excluding the first set, a fully implicit time integration scheme, analytical tangent approach and monolithic treatment of the local ODE system are used. Before the results are presented, the parameters, initial conditions and geometries for the FE analyses are introduced. The initial conditions are obtained by stimulating the single cell action potential models every 800 ms for 1000 cycles. Unless noted otherwise, the parameters and initial conditions present in Table 1 are used for the results related to the A-P model. Likewise, the parameters for the TNNP model are given in Table 2 complemented by the initial conditions shown in Table 3. Unless noted otherwise, all simulations used a rectangular cuboid geometry, as depicted in Figure 1, meshed with linear hexahedral elements. The CV is calculated through the equation where (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ) are points belonging to the geometry and t act (x i , y i , z i ) is the activation time of point (x i , y i , z i ). The points used in the cuboid geometry are (x 1 , y 1 , z 1 ) = (8,0,0) mm and (x 2 , y 2 , z 2 ) = (18,0,0) mm. The first point is chosen to be far enough from the stimulus region to avoid any interference from the stimulus. The second point is chosen to be far away from the first point to more accurately calculate the CV. We define electrical activation time t act as the time in which the transmembrane potential has depolarized above −40 mV. The t act values are obtained at the Gauss points and are then projected to the nodes to be used for the computation of the CV.

| Comparison of calculation approaches
In light of the different calculation approaches described in Section 2.2, several simulations are performed to compare the computational time of the approaches. The simulations are performed using an Intel® Xeon® CPU E5-2667 v4 @ 3.20GHz processor. For all simulations, we use a total time of 400 ms. Beginning with the A-P model, a total of five different calculation approaches are used: 1. fully implicit with analytical tangent, 2. semi-implicit, 3. fully implicit with numerical tangent, 4. explicit in r while implicit in ϕ with analytical tangent, 5. explicit in r while implicit in ϕ with numerical tangent. Only the first two of these approaches are worth comparing. The other three options give identical activation time distributions and similar computational time when compared to the fully implicit with analytical tangent case. As expected, the solution time increases with decreasing mesh size and decreasing time step, see Figure 2A,B. Of note, the semi-implicit case is always faster than the fully implicit case. The influence of the spatial and temporal discretization on the CV is depicted in Figure 3. Regarding the time step dependence, we note that the CV for the semi-implicit case approaches a converged value from below (larger time steps give a lower CV). On the other hand, the implicit cases approach the converged CV from above (larger time steps give a higher CV). As for the dependence of the CV on the element size, both the implicit and the semi-implicit case approach the converged CV from above (coarser meshes give a higher CV).
In a like manner, we consider five different calculation approaches for the TNNP model: 1. fully implicit monolithic with analytical tangent, 2. fully implicit monolithic with numerical tangent, 3. fully implicit staggered with analytical tangent, 4. fully implicit staggered with numerical tangent, 5. semi-implicit.
To determine the optimal values for the perturbations Φ p and c ion,p , a numerical analysis is carried out. Basically, we obtain the Euclidean norm of the difference between the derivatives obtained using analytical tangents and the derivatives obtained using numerical tangents at every time step in a single cycle with a single paced cell. The derivatives considered are the local tangent matrix K and the sensitivity ∂ Φ f Φ . The Euclidean norms are integrated over a single cycle of 800 ms to compare the effect of the perturbations. Therefore, the integrals can be represented as follows where * AT and * NT are the values obtained using the analytical and numerical tangent approach, respectively. The lowest value of the integral is used to determine the optimal value of the perturbation. The optimal values exhibited in Figure 4A,B, are Φ p = 10 −7 mV and c ion,p = 10 −8 , respectively. As seen in Figures 5A,B, the monolithic and staggered approaches give similar computational times using the same differentiation approach. The staggered with numerical F I G U R E 1 Cuboid geometry (dimensions in mm) used for the numerical experiments. The volume in orange represents the region where a stimulus is applied tangent option is less stable than the other fully implicit cases (note the absent data points for large time steps in Figure 5B). Moreover, the simulations with the numerical tangent approach take longer to compute. The semi-implicit case is much more unstable and requires a much smaller time step to converge, see the data point in Figure 5B. Nevertheless, at the smallest time step, the semi-implicit approach requires only about half the time to compute as the fully implicit monolithic with analytical tangent approach. The stability of the semi-implicit approach appears to depend on both the spatial and temporal discretization. For coarse meshes, a finer time step is required. This can be observed in Figure 6B, where the absent data points are due to failures in convergence. All simulations utilizing a fully implicit time integration scheme give identical activation time distributions and convergence from above (lower resolution gives a higher value) for the CV as the temporal and spatial resolution are refined. In Figure 6A, the CV is shown to converge from above (lower resolution gives a higher CV) for the fully implicit case. The semi-implicit case does not show any strong dependence of the CV on the time step but does show a dependence of the CV on the mesh size similar to that of the fully implicit approach, as depicted in Figure 6B.

| Conduction velocity convergence
This section explores how the convergence of the CV relative to the spatial and temporal discretization is dependent on material parameters. The material parameters can affect the conductivity, upstroke velocity and excitability. As  42 these are some of the main determining factors for the CV. We present results in the simple form of a percentage difference between the CV at a given discretization and the CV for the finest discretization CV finest . The CV percent difference CV dif is thus calculated Although we only consider the change of the CV relative to the finest discretization, altering the conductivity, upstroke velocity and excitability leads to a change in the magnitude of the CV. For the A-P model, the CV is increased for a larger conductivity d iso and for lower values of parameter α (related to the excitability). For the TNNP model, the CV is increased for a larger conductivity d iso , a larger maximal sodium channel conductance C max Na (related to the upstroke velocity) and a smaller parameter Ξ (related to the excitability).
Regarding the A-P model, Figure 7 shows how the convergence of the CV due to discretization refinement changes for increasing conductivity. As the conductivity increases, the CV dif reduces. Changes in the conductivity do not appear to have a large effect on the time step dependence of the CV. As for changes in α, the excitability parameter, increasing α leads to a decrease in the temporal dependence and no recognizable change in the spatial dependence of the CV dif , see Figure 8. As for the TNNP model, the time step has only a small effect on the CV, since the time step required for stability is already fine. Therefore, we consider only the spatial dependency of the CV. A time step of Δt = 0.05 ms is applied. As seen in Figure 9A, an increase in the conductivity causes a decrease in the mesh dependency. The contours display a linear appearance with regards to the logarithm of the conductivity d iso and the mesh size h. Next, considering the effect of the maximal conductance of the fast sodium current C max Na , a pronounced increase in the mesh dependency of coarse meshes can be perceived for larger values of C max Na , see Figure 9B. The mesh dependency appears to be unaltered for fine meshes, however, as seen in the plotted contours. Additionally, non monotonic convergence is exhibited for lower values of C max Na . Finally, the relationship between the excitability parameter Ξ, element size h and CV dif is depicted in Figure 10. Surprisingly, the CV mesh dependency is reversed for higher values of Ξ. That is, for lower values of Ξ, the CV converges from above, whereas for higher values of Ξ, the CV converges from below. As seen when comparing Figure 10A,B, the effects of Ξ are different for the hexahedral and tetrahedral meshes. For some values of Ξ, the hexahedral mesh approaches the CV from below, while the tetrahedral mesh approaches the CV from above. Figures 10A,B also show that conduction block (shown as values of −100%) occurred for higher values of Ξ. Furthermore, a non monotonic convergence is seen for a few values of Ξ in the hexahedral case.
For the sake of completeness, the benchmark of Niederer et al 19 is performed for both the A-P model and the TNNP model. In brief, a cuboid of dimensions 20x7x3 mm 3 is stimulated in one corner and the activation times are determined in the rest of the domain based on the propagation of the transmembrane potential. The results, presented in Figure 11, display the activation times on a line connecting the points (0,0,0) to (20,7,3). Even at the finest discretization, the activation times do not appear to be completely converged.

| Excitation stages in coarse and fine meshes
Signal propagation in FE electrophysiology occurs due to the sequential increase in the transmembrane potential at the nodes. Therefore, analysing the excitation in a specific part of the geometry can help to understand the excitation wave propagation across the domain. In this section, we investigate how the excitation evolves at a node of interest in the A-P model and the TNNP model in the case of coarse (h= 2.0 mm) and fine (h= 0.2 mm) meshes with a time step of Δt = 0.01 ms.
We investigate excitation at a node in the cuboid geometry with a linear hexahedral mesh. Since the governing equation of cardiac electrophysiology [Equation (2)] is spatio-temporal in nature, the excitation at a node is also dependent on space and time. To make meaningful observations, we reduce the dimensionality and, thus, the amount of information to display. For the given FE simulation, the transmembrane potential propagates in the form of planar waves with the x direction being normal to the plane. Due to this planar waveform, the spatial gradient in the y and z directions are negligible and the problem can be best regarded from a one dimensional perspective. Further reducing the information, we consider only three points in space: one node of interest (labeled N) and two neighboring Gauss points used for the numerical integration. We label these Gauss points based on their position relative to the wavefront. The Gauss point closest to the incoming wavefront is labeled GP1 and the Gauss point furthest from the incoming wavefront is labeled GP2. The node of interest is located on one edge of the rectangular cuboid at x = 10 mm, y = 0 mm and z = 0 mm. F I G U R E 9 Percent difference of the CV compared to the CV at the finest discretization CV dif using the TNNP model. Discretization is refined only in space, while (A) the conductivity d iso or (B) the maximal sodium current conductance C max Na are varied. The range of percentages covered by the colors is shown in the legend, while the contours are placed at 0.5% intervals. Values above 100% are ignored We display the evolution in time of the transmembrane potential at the node of interest and at its neighboring Gauss points for the A-P model in Figure 12 and the TNNP model in Figure 13. These figures also display the graphs of the ionic current f Φ and the spatial gradient in the x direction ∂ x Φ at the Gauss points. As can be seen in Figures 12 and  13, we divide the excitation into seven stages which are outlined in Table 4. These stages are categorized based on the effects of the ionic current f Φ and the diffusion term div(q) on the nodal value of the transmembrane potential Φ N . The first stage designates the case where GP1 is un-excited. Next, in Stage 2, due to a neighboring increase in the transmembrane potential, Φ GP1 increases as a result of conduction. During Stage 3, Φ GP1 reaches the activation threshold Φ act and thus now contributes to Φ N with a positive value of f Φ . The value of Φ GP2 has not yet reached the activation threshold Φ act . Further, for Stage 4, Φ GP1 and Φ GP2 have now reached the activation threshold and both contribute to raising Φ N . The diffusion term still has a positive contribution to Φ N . On the other hand, for Stage 5, GP1 and GP2 have a positive source term contribution but now the diffusion takes away from Φ N . In Stage 6, Φ GP1 and Φ N have reached the excited state. However, Φ GP2 will not be completely excited until the next node is excited and thus no diffusion is taking away from Φ N . Finally, GP1, N and GP2 are excited in Stage 7. The A-P model excitation in Figure 12 displays similarities and differences between the coarse and fine mesh. For both the coarse and fine meshes, the duration of Stage 2 is short, due to a low activation threshold. The switch from Stage 4 to Stage 5 can be seen to occur when the spatial gradient at GP1 is equal to the spatial gradient at GP2. This leads to the same influx as outflux of Φ, meaning that the diffusion has no contribution to the excitation at the node in this point of time. During Stage 3, the coarse mesh shows GP1 has an increasing ionic current f Φ contribution while GP2 has no contribution. By observing the time at which the excitation occurs in the coarse and fine mesh, one can see that the coarse mesh has an earlier excitation and, thus, a higher CV. Oscillations occur in Φ for the coarse mesh before C1 F1 C2 F2

C3 F3
F I G U R E 1 2 Graphs of (1) the transmembrane potential Φ, (2) the ionic current f Φ and (3) the spatial gradient in the x direction ∂ x Φ using the A-P model. The graphs for the coarse mesh are denoted with (C) while the graphs for the fine mesh are shown with (F). The evolution in time of the variables is depicted for a node of interest and the neighboring Gauss points GP1 and GP2. The variable Φ is defined at the nodes. At the Gauss points, Φ is obtained by linear interpolation between the nodes. For f Φ and ∂ x Φ, only the values at GP1 and GP2 are shown, since these variables are calculated only at the Gauss points. The stages outlined in Table 4 are depicted using their respective colors. The shades of the colors are altered depending on whether they represent GP1, N or GP2 and after depolarization. In the coarse mesh, the spatial gradient is greatly underestimated and contributes to changing the transmembrane potential over a longer period of time. Furthermore, the evolution of Φ and f Φ in the coarse mesh is not the same at GP1 and GP2. Finally, observing the excitation for the TNNP model in Figure 13, it can be seen that the excitation behavior is similar to that of the A-P model. The coarse mesh with the TNNP model also shows oscillations in Φ and a non uniformity in the response at each Gauss point for Φ and f Φ . The spatial gradient also contributes to changing the transmembrane C1 F1 C2 F2

C3 F3
F I G U R E 1 3 Graphs of (1) the transmembrane potential Φ, (2) the ionic current f Φ and (3) the spatial gradient in the x direction ∂ x Φ using the TNNP model. The graphs for the coarse mesh are denoted with (C) while the graphs for the fine mesh are shown with (F). The evolution in time of the variables is depicted for a node of interest and the neighboring Gauss points GP1 and GP2. The variable Φ is defined at the nodes. At the Gauss points, Φ is obtained by linear interpolation between the nodes. For f Φ and ∂ x Φ, only the values at GP1 and GP2 are shown, since these variables are calculated only at the Gauss points. The stages outlined in Table 4 are depicted using their respective colors. The shades of the colors are altered depending on whether they represent GP1, N or GP2 potential for a longer period of time and has a smaller magnitude, when comparing the coarse mesh to the fine mesh. The fine mesh is much more uniform with only a short time devoted to Stages 3 and 4.

| DISCUSSION
Following the aim of this paper, several simulation results are presented considering discretization in conjunction with other modeling aspects. The results show a strong dependence on the discretization and will be interpreted further in individual sections.

| Calculation approaches
The results for the different options for calculations at the material level help to get an idea of the advantages and drawbacks of the different methods. Considering the time integration scheme, the benefit of using an implicit solution is that it can be much more stable than the semi-implicit solution. As seen in the results, the semi-implicit scheme can fail at larger time steps. For the TNNP model, the ODE system is very stiff and requires a much smaller time step when a semi-implicit case is used. Roy et al 31 show that the TNNP model with the semi-implicit case is stable when Δt ≈ 0.001 ms and h ≈ 0.06 mm. This agrees with the results present in this work, where for Δt = 0.001 ms and 0.01 mm <h < 0.2 mm, the semi-implicit TNNP model is stable. However, the semi-implicit case is faster compared to the implicit case because it does not require a Newton-Raphson iteration at the local level. Jilberto and Hurtado report a significant speed-up when using the semi-implicit case as opposed to the implicit case with a modified A-P model. 43 In this work, the A-P model only displayed a minor speed-up, while, for the TNNP model, the speed-up was significant. Seeing as the semi-implicit A-P model was not much faster than the implicit case, this may imply that the speed-up is dependent on the chosen material parameters. When choosing which time integration scheme is the most optimal, the stability and the solution time need to be considered together for a given action potential model. Regarding the CV convergence, Hoermann et al 25 display a convergence behavior for the CV using a semi-implicit approach similar to what we exhibit in Figure 3B (convergence from below for temporal discretization and convergence from above for the spatial discretization). While, the CV convergence behavior for the fully implicit case agrees with the results from Wong et al 17 (convergence from above for both spatial and temporal discretization). It is worth noting that the semi-implicit case considered in this work is simple (using the classic forward Euler method) and that more sophisticated semi-implicit schemes can also provide enhanced stability and improved computational efficiency. Furthermore, implicit schemes also suffer from time step constraints for stability and may even, in some cases, be comparable to some semi-implicit approaches. The monolithic and staggered treatments when solving the local system of ODEs give very similar results in terms of the solution time. The monolithic solution is generally more stable since it accounts for the sensitivities arising from the coupling of variables and thus provides a more accurate update in the local Newton-Raphson iteration. The staggered solution is simpler since less sensitivities need to be computed but the stability of the solution can be diminished. The viability of the staggered approach for the TNNP model is connected to the fact that not all concentrations are coupled to each other. At large time steps, issues arise with the local convergence related to the lack of coupling.
The use of a numerical tangent can provide an approximation of the Jacobian matrix in cases when analytical differentiation is unwanted. Similar to other studies that utilize numerical tangents, 44,45 we found an optimal value exists at around 10 −8 for the perturbation. This value results from the tradeoff between accuracy and precision. Ideally, the smaller the perturbation the more accurate the approximation. Numerically, however, a limited numerical precision is used. As a result, if the perturbation is too small, not enough precision is available to make an accurate calculation. Although the simulations with the numerical tangent cases have larger computational times compared to the analytical tangent cases, using numerical tangents greatly reduces the complexity of the code by not requiring exact derivatives.

| Conduction velocity convergence
In the benchmark study by Niederer et al, 19 the dependence of the activation time on the spatial and temporal discretization is also explored. In their study, a 3D cuboid geometry, a transversely isotropic conductivity, the TNNP model and similar material parameters as those in this work are employed. The converged CV value from the benchmark study is 0.5 m/s (derived from the length of the specimen and the converged activation time). This value is slightly lower than the CV value obtained in Figure 6A (0.61 m/s) because only one dimension of the problem is considered, specifically, the dimension with the higher conductivity. In the benchmark results presented in Figure 11B, the activation time of the finest discretization at the far end of the cuboid (39.3 ms) are similar to the converged result of about 43 ms (8.6% different). The dependence of the CV and t act on the spatial discretization have a similar trend as that of the first research group model (or code) considered in the benchmark study (coarser meshes provide a higher CV and a lower t act ).
The parameters of the action potential models are shown to have a strong effect on the discretization dependency of the CV. This fact is noted in the literature. For example, Vincent et al propose the use of a metric for convergence, the Thiele modulus, that considers the element size, conductivity and upstroke velocity. 46 As such, the conductivity, upstroke velocity and excitation threshold are important factors which should be considered when comparing the discretization dependency between action potential models.
The conductivity appears to have a length scaling effect. That is, a change in conductivity shifts the relationship of the CV convergence with respect to the element size. For instance, when a higher conductivity is used, the CV converges at a larger element size. This is observed for both the A-P model and the TNNP model. This is most likely due to the conductivity being related to the spatial gradient that can develop. For a high conductivity, the transmembrane potential gets diffused easily and a large spatial gradient cannot develop. When the spatial gradient is low, the element size required to faithfully represent the wavefront through interpolation can be larger than when the spatial gradient is high. In fact, a limit exists for the magnitude of the spatial gradient that can be represented by linear elements. This limit arises from the electrophysiology and element size and is equivalent to the maximum difference in transmembrane potential divided by the element size h.
Increasing the upstroke velocity through parameter C max Na in the TNNP model leads to a larger mesh dependency for coarse meshes. This behavior is in accordance with what is observed for changes in conductivity. Increasing the upstroke velocity leads to a higher spatial gradient and, thus, a greater mesh dependency.
The activation threshold is shown to have a strong effect on the discretization dependency of the CV. The activation threshold influences the relationship between the diffusion, excitation and geometry. After all, the activation threshold defines how much diffusion must take place before excitation occurs. Moreover, it affects the way interpolation will overestimate or underestimate the transmembrane potential at the Gauss points. For the A-P model, the activation threshold parameter α reduces the time step dependency of the CV. This behavior has, to the best of our knowledge, not been reported before. Interestingly, the CV dif plot in Figure 8C has a very similar distribution as that of Figure 6A. This is likely due to the fact that the excitation threshold obtained for the value of α = 0.25 is similar to the excitation threshold present in the TNNP model. In the TNNP model, we explore the effects of changing parameter Ξ on the CV. High activation thresholds cause the signal propagation to rely more on diffusion which is underestimated in coarse meshes. This causes the coarse meshes with high activation threshold to have a lower CV for coarse meshes compared to fine meshes. In agreement with the mesh dependency findings, Pezzuto et al 39 demonstrates that the CV error convergence behavior is strongly dependent on the threshold potential.
The difference in response between the tetrahedral and hexahedral meshes is most likely related to the position of the Gauss points. In the tetrahedral elements, we use a single Gauss point in the middle of the element. For the hexahedral elements, we use 8 Gauss points, which means two Gauss points are present in each direction (a 2x2x2 configuration). Due to linear interpolation, a change in the location of the Gauss points corresponds to a change in the value of T A B L E 4 Stages of excitation Stage Defining condition the transmembrane potential at the Gauss points. Furthermore, since the value of the transmembrane potential is important for the excitation, the excitation properties are altered for different Gauss point locations. In light of the effects of Gauss point location on CV error, we are currently exploring the idea of modifying the quadrature by shifting the position of the Gauss points to reduce the CV error in coarse meshes. It is worth noting that Ξ possibly affects the upstroke velocity. Only the steady state variable of the sodium activation gate g ∞ m is affected by this parameter, while the other constants and gating variables for the fast sodium current I Na are left unchanged. The g m gate is related to the I Na activation while the g j and g h gates lead to inactivation. As such, this could lead to changes in the duration and magnitude of the fast sodium current I Na meaning the upstroke velocity could be different when Ξ is changed.

| Excitation stages in coarse and fine meshes
In the A-P model and TNNP model, the spatial gradient is greatly underestimated for the coarse mesh. An underestimated spatial gradient leads to a reduction in the magnitude of the diffusion of the transmembrane potential. However, this does not mean that the CV will be underestimated in coarse meshes. After all, the CV is higher in the coarse mesh for the A-P and TNNP models. In the coarse mesh, the response at GP1 is much different than the response at GP2. Along with the presence of oscillations, this leads to noticeable changes in the shape of the wavefront. The lack in uniformity and presence of oscillations also has effects at the material level where the internal variables (gating variables and concentrations) will evolve differently and present a different time history behavior. Regarding the excitation, stages 2 (starting with GP1 being un-excited and ending with GP1 being excited) and 3 (starting with GP1 being excited and ending with GP2 being excited) represent important moments in the excitation. Stage 2 determines how much diffusion must occur before excitation. While, stage 3 forces an excitation to occur at the node.

| Limitations
The cuboid geometry utilized for the simulations in this work is very simplified and results in wave propagation which is essentially one-dimensional. For a more comprehensive 3D study, the benchmark problem by Niederer et al 19 is more suitable. However, the benchmark problem requires isotropic adjustments of mesh size which can greatly increase the computational demand of the problem, especially when considering the finest mesh sizes used in this work which are an order of magnitude smaller than those in the normal benchmark problem. Tetrahedral meshes also require isotropic adjustments of the mesh size. As a result, the mesh used for the results in Figure 10B was coarse and the accuracy for the finest discretization was still poor. Lastly, the use of the monodomain model is a simplification, since it assumes the intracellular and extracellular conductivities are proportional and can lead to different activation distributions compared to the bidomain model. Furthermore, the monodomain model does not consider extracellular potentials, which are often of clinical relevance.

| CONCLUSION
The experiments considering the various calculation approaches show the stability and computational time required for each approach. In terms of stability, the fully implicit, monolithic and analytical tangent approach provides the most stable solution. The semi-implicit case, on the other hand, can demand a much smaller spatial and temporal resolution. With regards to the computational time, the semi-implicit time integration scheme provides a faster solution compared to the fully implicit case. This speed-up is more pronounced when using the TNNP model. Moreover, for the TNNP model, the analytical tangent approach is faster than the numerical tangent approach. Lastly, the solution time required when using the monolithic treatment of the local system of ODEs is comparable to that required using the staggered approach.
The dependence of the CV on the spatial and temporal discretization is shown to be strongly influenced by the material parameters. The conductivity is shown to have a length scale effect. In other words, a higher conductivity allows the CV to converge at a larger mesh size. The parameter C max Na in the TNNP model affects the upstroke velocity. The upstroke velocity can strongly affect the spatial gradient that develops. Lastly, the excitation threshold affects the behavior of the CV error convergence. A high activation threshold can lead to an underestimation of the CV for coarse meshes while a low activation threshold causes an overestimation of the CV. Additionally, an increase in the activation threshold is shown to lower the dependency of the CV on the time step for the fully implicit A-P model. CV error in linear elements arises due to two main factors which are more prevalent in coarse meshes than fine meshes. First, coarse meshes with linear elements cannot capture steep spatial gradients, which leads to an underestimation of the CV. Second, conduction is affected by an overestimation of the transmembrane potential at the Gauss points due to linear interpolation. This overestimation can avoid signal block in coarse meshes but is the reason why the coarse meshes overestimate the CV when compared to fine meshes.