The Finite Element Method (FEM) is widely-used for simulating cardiac electrophysiology mainly due to its robustness and accuracy [heidenreich2010adaptive]. However, accuracy is substantially deteriorated when the mesh undergoes a large deformation or if it does not satisfy specific quality criteria [eppstein2001global]. For this reason, in recent years there has been growing interest in alternative meshfree methods that partially or completely alleviate the mesh-related limitations of FEM. Several meshfree methods have been investigated so far in the context of cardiac electrophysiology.
The element-free Galerkin (EFG) method offers high convergence rate [belytschko1994element] and it has been used successfully in several applications [belytschko1996meshless, liu2003mesh]. It has been used in cardiac electrophysiology to solve the cardiac monodomain equation using a meshfree representation of the Auckland heart model [zhang2012meshfree]. EFG employs the Moving Least Squares (MLS) approximation [lancaster1981surfaces] for the solution of the Galerkin weak form. Therefore, the imposition of essential boundary conditions requires special treatment due to the lack of the Kronecker delta property in MLS. To solve this issue, Maximum Entropy (MaxEnt) approximants [arroyo2007local, millan2015cell] have been proposed as an alternative to MLS in EFG [mountris2020cell]. Due to the weak Kronecker delta property of MaxEnt, the influence of internal nodes on the boundary of the domain of interest is eliminated and essential boundary conditions can be imposed directly as in FEM. However, the MaxEnt approximation functions are significantly more complex than MLS. As a result, a large number of quadrature points is required to deal with the quadrature accuracy problem of meshfree methods [dolbow1999numerical, babuvska2008quadrature].
On the other hand, the Mixed Collocation Method (MCM) [mountris2020radial, mountris2021cardiac] is a purely meshfree method that avoids the quadrature accuracy problem. It has been used for the simulation of the cardiac monodomain model demonstrating results in good agreement with FEM [mountris2019novel, mountris2020next]. MCM is a variation of the Meshless Local Petrov-Galerkin method [atluri2006meshless]
where the Dirac delta distribution is used as test function and interpolation is applied both on the field function and its gradient[zhang2014meshless, li2008mlpg, atluri2006meshless]. As a result, the integrals in the Petrov-Galerkin weak form are replaced with nodal summation. In addition, essential boundary conditions are imposed directly through collocation. However, being a collocation method, it usually requires the construction of support domains with a large number of collocation points to avoid inaccuracy during the imposition of natural boundary conditions [Libre2008] and may not be as efficient as FEM in large-scale 3D problems.
Very recently, a novel meshfree technique going by the name of Fragile Points Method (FPM) has been added to the artillery of meshfree methods [dong2019new]. FPM uses local, simple, polynomial, discontinuous functions [liszka1980finite] as trial and test functions for the solution of the Galerkin weak form. Using these polynomial trial and test functions, the integration in the Galerkin weak form becomes trivial. Due to the simplicity of the test and trial functions, single-point integration is sufficient to compute integrals with high accuracy in FPM. We refer to table 1 for a comparison of integration in FPM and other meshfree methods. Moreover, both essential and natural boundary conditions can be imposed as in FEM. However, due to the discontinuity of the trial and test functions, the assembly of the FPM point stiffness matrices leads to an inconsistent global stiffness matrix. To remedy this issue, the numerical flux corrections, which are widely used in Discontinuous Galerkin methods [arnold2002unified], are employed in FPM to obtain a consistent, sparse and symmetric global stiffness matrix.
Despite being only recently introduced, FPM has been demonstrated to be an accurate and efficient method with many desired properties (i.e. accurate integration, exact imposition of boundary conditions) and has been already employed to solve linear elasticity [yang2021simple], heat conduction [guan2020heat, guan2021meshless] and flexoelectric problems with crack propagation [guan2020new].
In the present study, we employ FPM for the solution of the cardiac monodomain model. Our motivation is to provide a meshfree alternative to FEM for cardiac electrophysiology simulation, while maintaining accuracy and efficiency. The paper is structured as follows. In section 2, we describe the theoretical aspects of FPM and the derivation of the cardiac monodomain models. In section 3, we present several numerical examples where FPM is applied to simulate electrophysiology in both 2D and 3D benchmark problems, as well as in a large scale biventricular geometry under myocardial infarction conditions. Finally, in section 4, we provide our concluding remarks.
2 Cardiac monodomain model FPM formulation
2.1 The cardiac monodomain equation
We consider the cardiac monodomain model for the simulation of electrical impulse propagation in the human heart, which is governed by the following reaction-diffusion partial differential equation (PDE)[keener2009mathematical]:
where denotes the transmembrane voltage time derivative, the total ionic current, the cell capacitance per unit surface area and
the diffusion tensor.and are the domain of interest and its boundary and
is the outward unit vector normal to the boundary.
The diffusion tensor is given by:
where denotes the diffusion coefficient along the cardiac fiber direction, the cardiac fiber direction vector, the transverse-to-longitudinal conductivity ratio,
the identity matrix andthe tensor product operator.
We employ the operator splitting technique to obtain the decoupled system of Equation (1):
The decoupled system can be solved efficiently by applying either the Godunov (first-order) or Strang (second-order) methods [schroll2007accuracy]. In the following, we consider the ionic currents of equation (3) normalized by the capacitance .
2.2 Computation of trial and test functions in FPM
FPM is a meshfree method where trial and test functions are established on arbitrarily distributed points in the domain of interest . Unlike mesh-based methods, such as FEM, connectivity information is not required. Simple, local, polynomial trial functions are defined in compact support domains, formed by partitioning the domain in conforming and nonoverlapping subdomains around each point. The partition is not unique and can include subdomains of arbitrary shape, such as polygons that can be obtained by the Voronoi diagram partition.
We define the trial functions at each subdomain in terms of the transmembrane voltage and its gradient . The trial function for subdomain that contains the point is given by:
where is the coordinate vector of a point in , is the value of at , is the coordinate vector of and (i.e. the voltage gradient at ) is unknown. To compute , we employ the Generalized Finite Difference (GFD) method [liszka1980finite] to minimize the weighted discrete L norm:
where denotes the coordinate vector of , denotes the value of at and denotes the value of the weight function at (), with being the number of points in the support domain of . We should note that a compact support domain for is defined by the points in the first ring of adjacent subdomains to the subdomain of (Figure 1). Assuming constant weight functions, we obtain the transmembrane voltage gradient at by:
Using Equation (6), the trial function can be obtained by:
where denotes the shape function of in :
Since the shape function is defined in each subdomain independently, it can be discontinuous at the internal boundaries. The Galerkin weak form in FPM is established by constructing both the trial and test functions using Equation (8). It should be noted that due to the discontinuity of the trial and test functions, the Galerkin weak form will lead to an inconsistent matrix and inaccurate results. Therefore, numerical flux corrections, which are common in Discontinuous Galerkin Finite Element Method [arnold2002unified], are introduced in FPM to remedy this issue.
2.3 Numerical flux corrections to remedy the inconsistency
We address the inconsistency issue by employing the Interior Penalty (IP) numerical flux corrections [mozolevski2007hp]. We start by writing the Galerkin weak form of Equation (3) for each subdomain , which is given by:
where and are the test and trial functions, respectively, is the boundary of the subdomain , denotes the outward normal vector to , and denotes the set of internal and external boundaries, i.e., , where is the set of internal boundaries.
Next, the jump operator and average operator are used to sum Equation (10) over all subdomains [guan2021meshless]:
where the jump operator and average operator are given, , by:
When , is a unit vector normal to and pointing outward from .
Additionally, and when is the exact solution since there is no jump in the internal boundaries. As a result, and we can replace the term in Equation (11) with without affecting the accuracy of the formula.
Finally, the internal penalty numerical flux is applied on with penalty parameter to obtain the consistent FPM formula [guan2020heat]:
where is a boundary-dependent parameter with the unit of length. In this work, we define as the distance between points , when . The penalty parameter , , has the same units as and is independent of the boundary size. It should be noted that the penalty parameter should be large enough to ensure stability, but excessively large values should be avoided since they may cause a condition number problem. An extensive discussion on recommended values for the penalty parameter can be found in [guan2020heat]. In this work, we use given by:
where is a penalty coefficient, denotes the volume of the cell containing the of the points in the support domain, and is the mean value of the diagonal entries in the diffusion tensor of the point.
2.4 Numerical implementation
The formula of FPM can be written in matrix form:
where and denote the global normalized capacity and diffusion matrices, respectively, and is the unknown vector collecting the nodal values of the transmembrane potential.
To assemble the global matrices and , we substitute the shape function for and and matrix for and in equation (14) to obtain the point normalized capacity matrix , the point diffusion matrix , and the internal boundary diffusion matrix :
The assembly of global matrices and is performed as in FEM, where is obtained by assembling all the individual capacity point matrices , and is obtained by assembling all the point diffusion and internal boundary diffusion matrices . It should be noted that the resulting global diffusion matrix is symmetric, sparse, and positive definite.
3 Numerical examples
Simulations for the numerical examples that are presented in the following were performed on a laptop with Intel® Core™i7-4720HQ CPU and 16 GB of RAM. The efficiency of FPM was evaluated by comparing the execution time of FPM simulations with the execution time of FEM using linear elements. The execution time for both FPM and FEM for the considered numerical examples is summarized in Table 2.
|Electrical propagation in a 2D ventricular tissue (subsection 3.1)|
|Acetylcholine-induced effects in human atrial electrical activity (subsection 3.2)|
|Electrical propagation in a benchmark 3D cuboid geometry (subsection 3.3)|
|Simulation of electrical activation in 3D biventricular infarction model (subsection 3.4)|
3.1 Electrical propagation in a 2D ventricular tissue
We considered a 4x4 cm human ventricular tissue with fibers aligned parallel to the x-axis. The longitudinal diffusion coefficient was cm/ms and the transverse-to-longitudinal conductivity ratio was . Electrophysiology was modeled by the O’Hara et al. [ohara2011simulation] human ventricle action potential (AP) model for epicardial cells. Periodic stimuli of duration ms and amplitude equal to twice the diastolic threshold were applied on the left side of the tissue ( cm) at a frequency Hz. AP propagation was simulated for a total time s, after achievement of steady-state, using the dual adaptive explicit time integration method (DAETI) [mountris2021dual] with time step ms.
Solutions obtained by FPM with different penalty coefficients, , were compared to a solution obtained by FEM with bilinear isoparametric elements (Figure 2). We considered four different nodal discretizations with spacing mm, each of which is presented in one of the panels of Figure 2. Differences between FEM and FPM solutions were evaluated in terms of conduction velocity () and AP duration (APD) at 90% repolarization.
The maximum percentage difference between FPM and FEM in terms of was 116.9%, while in terms of it was 2.7% for nodal spacing mm and . The minimum percentage difference between FPM and FEM in terms of was 7.1%, while in terms of it was 0.0% for nodal spacing mm and . For a given nodal discretization, the best agreement between FEM and FPM solutions was always obtained for . For penalty factor values and , the values obtained by FPM with coarse nodal discretizations mm and mm were in better agreement with the value obtained by FEM for a dense nodal discretization ( mm).
To evaluate the effect of the penalty coefficient, we performed a convergence analysis where we evaluated the relative error () in for the solutions obtained by FPM with with respect to the FEM solution for nodal spacing mm:
where is the value obtained by FPM or FEM for each of the tested nodal discretizations and is the value obtained by FEM for a dense nodal discretization with mm. The relative error, , for FPM with and was smaller than for FEM for coarse discretizations, but it was larger for finer discretizations (Figure 3). For FPM with and , decreased monotonically and its values were smaller than for FEM for practically all nodal spacing values. In the following numerical examples, we used .
3.2 Acetylcholine-induced effects in human atrial electrical activity
Acetylcholine-induced APD shortening in atrial myocytes facilitates the initation and perpetuation of atrial fibrillation [bayer2019acetylcholine, smeets1986wavelength], which is a major risk factor for ischemic stroke. The parasympathetic neurotransmitter acetylcholine (ACh) shortens APD by activating the ACh-sensitive inward rectifier potassium current, I. Simulation of ACh-induced alterations in atrial electrophysiology is of interest in atrial fibrillation research [celotto2021location].
Here, we considered a 2D atrial tissue of cm with fibers aligned parallel to the X-axis in which ACh was distributed homogeneously throughout the tissue. The atrial myocyte electrophysiology was described using the Maleckar et al. model [maleckar2008mathematical]. The longitudinal diffusion coefficient was set to cm/ms and the transverse-to-longitudinal conductivity ratio was set to . The same periodic stimulation protocol as in subsection 3.1 was simulated for s and the DAETI method with time step ms was used to numerically solve AP propagation.
Simulations were performed using FPM with penalty coefficient considering a regular nodal distribution with mm and were compared with FEM simulations using bilinear isoparametric elements. The induced APD shortening was recorded at the center of the tissue for different ACh concentrations, M (Figure 4). APD shortening was found to range from 19 ms for the lowest ACh concentration to 38 ms for the highest ACh concentration. FPM results were found in good agreement with FEM, with differences of up to 2 ms. These results are in line with previous studies reporting APD shortening of up to 40 ms for M [bayer2019acetylcholine, celotto2019calcium, celotto2020sk].
3.3 Electrical propagation in a benchmark 3D cuboid geometry
In this example, we solved the standard benchmark problem for verification of cardiac tissue electrophysiology simulators described in [niederer2011verification]. The problem considered the electrical stimulation of a 3D cuboid of human ventricular tissue with dimensions mm and cardiac fibers parallel to the Z axis. Epicardial cell electrophysiology was described by the Ten Tusscher et al. model [ten2006alternans]. The longitudinal diffusion coefficient was set to cm/ms and the transverse-to-longitudinal ratio was set to . A periodic stimulus with frequency Hz, amplitude and duration ms was delivered at a cubic region with dimensions mm located at corner P1 (Figure 5).
The activation time at the corners (P1 – P8) and the center (C) of the cuboid were recorded for three different nodal discretizations, mm, for a simulation using FPM with penalty coefficient . Integration was performed using the DAETI method with time step ms. The activation times obtained by FPM were compared with those of FEM (Table 3) previously reported in [mountris2021dual] and validated with the reported activation times in [niederer2011verification].
The activation times for FPM and FEM for nodal spacing mm were found to be in good agreement. For coarse nodal discretizations, the FEM solution led to larger activation times as compared to FPM, especially at points P5–P8. Activation times at P7 for FEM simulations were found to be larger for the spacing mm than for the spacing mm. On the other hand, this difference was found to be for FPM simulations. For a nodal spacing mm, activation times were notably closer to those of mm when using FPM than when using FEM. These results demonstrate the higher convergence of FPM compared to FEM and the improved accuracy for coarse nodal discretizations, in good agreement with the findings in subsection 3.1.
|FPM activation time (ms)|
|FEM activation time (ms) as in [mountris2021dual]|
3.4 Simulation of electrical activation in a 3D biventricular infarction model
We employed the FPM method to simulate AP propagation in a 3D biventricular model under myocardial infarction conditions and we compared the activation pattern obtained by FPM with the one obtained by FEM. We considered a simplified model where zero conduction was assumed for the scar tissue, while border zone effects were not taken into account.
The biventricular anatomy was constructed using ex vivo diffusion weighted imaging (DWI) of a porcine heart. The two ventricles of the heart were manually segmented and a scar was introduced at the septal-anterior wall of the left ventricle (Figure 6). A tetrahedral mesh (nodes: 70521, elements: 311150) was generated from the segmented data using iso2mesh [fang2009tetrahedral]. FPM cells were generated for the nodes of the tetrahedral mesh using a dual polyhedral mesh generation algorithm [garimella2014polyhedral, kim2014efficient]. Fiber direction was determined by computing the diffusion tensors using an algorithm based on Riemannian distances [barmpoutis2010unified].
Cell electrophysiology was modeled by the O’Hara et al. model [ohara2011simulation] considering endocardium:midmyocardium:epicardium with 50:20:30 ratio. Pacing was performed by applying a periodic stimulus with amplitude equal to twice the diastolic threshold and frequency Hz at two tested stimulation regions: one at the base of the model located at the anterior wall of the right ventricular base near the septum and the other one at the apex of the left ventricle (Figure 6).
The local activation time (LAT) at each node of the model was computed using FPM with penalty coefficient and it was compared with the LAT value obtained from a FEM simulation. Mean LAT for basal pacing was 170 ms for FPM and 173 ms for FEM, while mean LAT for apical pacing was 151 ms for FPM and 148 ms for FEM. The LAT histograms presented in figure 7 demonstrate the good agreement between FPM and FEM simulations, with FPM rendering a valuable alternative to FEM for large scale cardiac electrophysiology simulations.
4 Concluding remarks
In this work, we presented the Fragile Points Method (FPM) for in silico cardiac electrophysiology applications. Despite being common for meshfree methods to sacrifice efficiency for accuracy and vice versa, FPM was proven to achieve similar accuracy and efficiency to FEM (Table 2).
We found that solutions obtained by FPM were in good agreement with FEM for both 2D and 3D scenarios, allowing to simulate action potential propagation with high accuracy under different physiological and pathological conditions. FPM demonstrated higher convergence than FEM (Table 3), in line with previous findings. This may have an important role in large scale applications where FPM could increase time efficiency without losing accuracy. By adjusting the penalty coefficient improved solutions could be obtained even for coarse discretizations. However, for large values of the penalty coefficient () accuracy could be deteriorated for fine discretizations. We found to lead to accurate results for all discretizations.
Finally, the ability of FPM to provide similar accuracy and efficiency to FEM without requiring mesh connectivity information renders the method an interesting alternative to FEM, particularly for personalized image-based modeling applications.
This work was supported by MCIN/AEI/10.13039/501100011033 (Spain) through project PID2019-105674RB-I00, by the European Research Council under grant agreement ERC-StG 638284, by the European Union’s H2020 Program under grant agreement No. 874827 (BRAV3) and by European Social Fund (EU) and Aragón Government through BSICoS group (T39_20R). Computations were performed by the ICTS NANBIOSIS (HPC Unit at University of Zaragoza).