## 1 Introduction

In the last decades, computational modeling and simulation has taken a growing role as a method to deepen the understanding of cardiac function in health and disease [lopez2015three, niederer2019computational]. Novel in silico models of increasing complexity are continuously being developed to simulate the electrophysiology [sampedro2020characterization] and mechanics of the heart [rama2018towards] from the cell [pueyo2016experimentally] to the whole-organ level [chabiniok2016multiphysics]. At the tissue and organ levels, electrophysiology is simulated by using the well-known bidomain [tung1978bidomain] and monodomain [keener2009mathematical] models. The latter is a simplified version of the former under the assumption of equal anisotropy ratios for the intracellular and extracellular spaces. The monodomain model is more computationally efficient than its bidomain counterpart and is able to produce accurate transmembrane potential values in the absence of extracellularly applied currents [potse2006comparison].

Commonly, state-of-the-art simulators [mirams2013chaste, vigmond2003carp] employ the Finite Element Method (FEM) to solve either the bidomain or the monodomain model for the simulation of cardiac electrophysiology. Despite the fact that FEM is a mature and robust numerical method, its requirement for a good-quality mesh poses challenges to generate realistic heart models with reasonable computational cost especially when cardiac mechanics are also accounted for and large deformations are aimed to be simulated [lluch2019breaking].

Alternative meshfree methods can alleviate the mesh requirement and have been proposed for both cardiac electrophysiological [lluch2017smoothed, zhang2012meshfree, mountris2019novel] and mechanical [lluch2019breaking, legner2014studying] simulations. Among the different proposed meshfree solutions, models based on the Element-Free Galerkin (EFG) method offer high convergence rate and high resolution of localized steep gradients [belytschko1994element]. Nevertheless, special treatment for the imposition of essential boundary conditions is required since the approximation functions do not possess the Kronecker delta property. Recently, Cell-based Maximum Entropy (CME) approximants were used in EFG to alleviate this limitation [mountris2020cell].

CME possesses the weak Kronecker delta property where approximation functions of internal nodes vanish on the boundaries. Therefore, essential boundary conditions can be imposed directly as in FEM. However, the CME approximants give rise to complex integrals requiring a large number of quadrature points for accurate integration that may lead to increased computational cost. Similarly, the computational cost of methods based on the Smoothed Particle Hydrodynamics (SPH) may be significantly higher than that of mesh-based methods. Furthermore, the standard SPH formulation may imply inaccurate computation of gradients of constant and linear fields (first-order incomplete approximation) [chen1999corrective]. To overcome these problems, the total Lagrangian formulation of SPH [bonet1999variational] and gradient normalization [chen1999corrective] were applied in [lluch2017smoothed, lluch2019breaking] to accurately simulate the propagation of the electrical impulse in the heart and cope with large deformations in the context of cardiac mechanics.

In the present study, we propose the Mixed Collocation Method (MCM) as an alternative to mesh-based and previously mentioned meshfree methods for cardiac electrophysiology simulation. MCM is based on the Meshfree Local Petrov-Galerkin (MLPG) method [atluri2004meshless, atluri2005basis]

. MLPG implies quadrature over locally-defined domains providing the flexibility to select the trial and test functions from different spaces. In the mixed collocation variant of the MLPG method, the Dirac function is used as test function and both the field function and its gradient are interpolated by the trial function. As a result, the computational cost is decreased since the local integration is reduced to nodal summation. Moreover, while collocation methods suffer from inaccurate imposition of Neumann boundary conditions

[libre2008], the accuracy is ameliorated in MCM due to the interpolation of the field function’s gradient. In the seminal work on MCM [atluri2006], the moving least squares (MLS) approximation [lancaster1981surfaces] was used as the trial function. Recently, the radial point interpolation (RPI) [wang2002point] was proposed as an alternative to MLS in MCM [mountris2020radial]. It was demonstrated that accuracy is improved when RPI trial functions replace MLS. This was mainly attributed to the delta Kronecker property of RPI that allows the direct imposition of essential boundary conditions in contrast to MLS.Here we investigate the application of MCM with interpolating trial functions for the solution of the monodomain equation in a series of cardiac electrophysiology simulations in 2D and 3D domains. We evaluate RPI as well as moving Kriging interpolation (MKI) [gu2003moving] and we compare the obtained solutions with state-of-the-art FEM solution. The structure of the paper is the following. In section 2, we derive the form of the cardiac monodomain equation using the MCM method. In section 3, we describe the mathematical formulation of RPI and MKI interpolations. In section 4, we evaluate the solution of the monodomain model with the MCM method in several 2D and 3D problems and we report on the accuracy and efficiency of the method in comparison to FEM. Finally, in section 5, we discuss some concluding remarks.

## 2 Mixed collocation form of the monodomain equation

The propagation of an electrical impulse in the heart is modeled through the monodomain model by the reaction-diffusion equation:

(1) |

where is the time derivative of the transmembrane potential , is the total ionic current and is the cell capacitance per unit surface area. and denote the domain of interest and its boundary, respectively, and

is the outward unit vector normal to the boundary.

is the diffusion tensor given by:

(2) |

where denotes the diffusion coefficient along the cardiac fibers, is the transverse-to-longitudinal conductivity ratio, denotes the cardiac fiber direction vector,

is the identity matrix and

is the tensor product operation.The transmembrane potential

is characterized by a rapid upstroke phase followed by a slow repolarization period, thus requiring to solve a system of stiff ordinary differential equations (ODEs) to reproduce it accurately

in silico. Therefore, it is common for large scale tissue simulations to decouple Equation (1) using the operator-splitting method [qu1999advanced] to improve computational efficiency. By applying operator-splitting, the decoupled monodomain equation is given by:(3) |

### 2.1 Deriving the mixed collocation form

To derive the mixed collocation form of the monodomain model we consider the discretization of the domain into a cloud of arbitrarily distributed field nodes. We write the diffusion term from Equation (3) for each field node in the set of nodes as:

(4) |

where

(5) |

denotes the transmembrane potential flux at the field node . Interpolating the transmembrane potential and the transmembrane potential flux we obtain:

(6) |

(7) |

where is the number of field nodes in the local support domain of the node . is the component of the vector containing the basis function of the meshfree approximation at each of the nodes in the local support domain of . denotes the transmembrane potential at node of the local support domain of and the corresponding transmembrane potential flux vector.

By introducing Equation (6) into Equation (5), we can express the transmembrane potential flux in terms of the nodal transmembrane potential as follows:

(8) |

or in matrix form and removing the dependence on for simplicity:

(9) |

where contains the spatial partial derivatives of the meshfree basis function at the nodes in the support domain of node scaled by .

By grouping the equations for all field nodes , with , the following equation in matrix form can be obtained:

(10) |

Finally, introducing Equations (6–8) into Equation (4), the mixed collocation formulation of the monodomain model’s diffusion term can be obtained in terms of the transmembrane potential:

(11) |

By grouping the equations for all nodes and removing the dependence on , as above, Equation (11) can be written in the equivalent matrix form:

(12) |

where is the sparse matrix collecting the basis functions and denotes the stiffness matrix.

### 2.2 Boundary conditions imposition

For the monodomain model, the domain is assumed to be isolated in the sense that no current can flow in or out of the boundary . To model electrical isolation, we enforce the Neumann boundary conditions (BC) in mixed collocation using the penalty method described in [overvelde2012moving]. From Equations (1) and (5), the Neumann BC imposition on the nodes of the Neumann boundary at a given time , where , can be written in matrix form as:

(13) |

where is the vector collecting the transmembrane potential fluxes at nodes and is the matrix containing the normal vectors given by:

(14) |

The Neumann BCs are enforced at the nodes by multiplying Equation (13) with the penalty factor and adding it to Equation (10) to obtain:

(15) |

By rearranging terms, Equation (15) can be written as:

(16) |

where is the identity matrix and . Combining Equations (12) and (16), the matrix form of the monodomain model’s diffusion term is given by:

(17) |

where . Here, the superscripts and connote the row entries of the matrices and that correspond to the nodes on and the nodes in , respectively, such that . The value of the penalty factor should be sufficiently large to ensure accurate imposition of the boundary condition, as instability issues may arise if is too large. In this study, we used as it was found to be the optimal value in [mountris2020radial].

## 3 Interpolating meshfree approximants

One of the advantages of MCM, being a meshfree method, is the flexibility that offers on the choice of the trial function . In this work, we consider only trial functions that possess the delta Kronecker property, namely the radial point interpolation (RPI) [wang2002point] and the moving Kriging interpolation (MKI) [gu2003moving].

### 3.1 Radial point interpolation

The RPI trial function for a field node is obtained by:

(18) |

where

is the radial basis function (RBF) for node

. For a node in the support domain of (with included in its support domain), is given by:(19) |

Different RBFs, such as multiquadric, Gaussian, etc., can be used. In this work, we used the multiquadric RBF (MQ-RBF). The value of the MQ-RBF for the 3D case, is calculated as:

(20) |

where denotes the Euclidean distance between nodes and , and and are positive-valued shape parameters of the MQ-RBF. For spherical support domains, the shape parameter is given by:

(21) |

where denotes the radius of the support domain of node and is a dimensionless constant. RBF fail to reconstruct exactly a linear polynomial field, therefore, the RPI is enriched with the linear polynomial basis to ensure continuity. For 3D problems, is given by:

(22) |

Finally the matrix is given by:

(23) |

where is the number of components of the polynomial basis ( for linear in 3D). and

denote the RBF and polynomial basis moment matrices:

(24) |

### 3.2 Moving Kriging Interpolation

The moving Kriging interpolation (MKI) has similar interpolation properties to RPI but it does not require polynomial enrichment to ensure continuity. The trial function at node is given by:

(25) |

where is the linear polynomial basis defined in Equation (22) and denotes the correlation function for node . For a node in the support domain of , is given by:

(26) |

where is the value of the correlation function at the node of the support domain. In this work, we use the MQ-RBF as the correlation function (Equation 20). The matrices and are obtained by:

(27) |

where is the identity matrix, is the moment matrix of the linear polynomial basis given by Equation (28) and is the correlation matrix for the nodes in the support domain of given by:

(28) |

## 4 Numerical examples

In this section, we investigate the accuracy and efficiency of MCM using both RPI and MKI as trial functions. We consider regular and irregular nodal distributions in 2D tissue sheets and 3D tissue slabs as well as in realistic anatomical models and we perform a comparison of MCM and FEM simulation results. MQ-RBFs for both RPI and MKI trial functions are constructed using and in 2D simulations. In 3D, they are constructed using and , as these combinations of parameters are found to minimize the difference with FEM. In all examples, time integration is performed explicitly using the dual adaptive explicit time integration (DAETI) algorithm [mountris2020dual] with adaptive time step ms. Human ventricular cellular electrophysiology is represented by the O’Hara et al. cell model [ohara2011simulation]. Simulations are performed on a laptop with Intel^{®} Core™i7-4720HQ CPU and 16 GB of RAM.

### 4.1 Electrical propagation in a 2D tissue sheet

We consider a cm human ventricular tissue sheet, where transmural heterogeneities are included by defining endocardial, midmyocardial and epicardial regions in a 5:2:3 ratio. The cardiac fiber direction vector is considered to be parallel to the -axis. We use a diffusion coefficient value cm/ms and a transverse-to-longitudinal conductivity ratio . Stimuli with amplitude twice the diastolic threshold, period s and duration ms are applied on the left side of the tissue ( cm). The propagation of the action potential (AP) is simulated for a total time s.

We compare the MCM solution with RPI and MKI trial functions against FEM simulation results using bilinear isoparametric elements. We consider regular nodal discretizations and quadrilateral meshes with nodal spacing cm. The considered support domains in the meshfree approximation have size , with . The generated APs at the center of the tissue sheet ( cm, cm) in the time interval s for the different nodal spacing values are shown in Figure 1. We quantify the differences between MCM and FEM solutions in terms of mean transmembrane potential difference (TPD). Mean TPD between FEM and MCM with MKI trial functions was TPD = {3.111, 0.339, 0.401, 0.583} mV while mean TPD between FEM and MCM with RPI trial functions was TPD = {3.112, 0.340, 0.400, 0.582} mV for nodal spacing cm. The efficiency of each simulation is evaluated in terms of execution time in Figure 3a.

### 4.2 Electrical propagation in a 3D tissue slab

We investigate the effect of the support domain’s dilatation coefficient by computing the normalized root mean square (NRMS) error between MCM and FEM solutions for a cm tissue slab. The tissue is assumed to be composed of epicardial ventricular cells. Stimuli of amplitude twice diastolic threshold, period s and duration ms are applied onto the left side of the tissue slab ( cm). The tissue slab is discretized with cm and varying dilatation coefficient . The NRMS error is computed using the formula:

(29) |

where and denote the transmembrane potential value at node computed with MCM and FEM, respectively. The NRMS error convergence plot is given in Figure 2.

The maximum NRMS error of the MCM solutions with RPI or MKI trial functions is obtained for and is equal to 0.262 and 0.263, respectively. The minimum NRMS error is obtained for and it is equal to 0.056 for RPI and 0.057 for MKI trial functions. The execution time for the simulations with varying dilatation coefficient is summarized in Figure 3b.

### 4.3 Electrical propagation in a 3D biventricular geometry under left bundle branch block conditions

We simulate electrical propagation in a porcine cardiac biventricular model under healthy (baseline) and left bundle branch block (LBBB) conditions. The biventricular model is represented by a tetrahedral mesh (273,919 nodes, 1,334,218 elements). It is provided by the CRT-EPiggy19 challenge [camara2019best, pop2020statistical] and is part of an LBBB dataset for experimental studies of cardiac resynchronization therapy [rigol2013development, iglesias2016quantitative]. We compute the direction of the myocardial fibers using a rule-based method [doste2019rule].

We define a diffusion coefficient value cm/ms in the longitudinal direction of the myocardial fibers and a transverse-to-longitudinal conductivity ratio . The electrophysiology of the ventricular myocardial tissue is represented using the O’Hara et al. model, as in previous examples. For a portion of connective tissue within the ventricles, the diffusion coefficient is reduced by a factor of 3 and the electrophysiology is modeled using the MacCannell et al. active fibroblast model [maccannell2007mathematical].

To identify a set of points with the earliest activation across the ventricles, we coupled the biventricular model with a network of Purkinje fibres generated using a fractal-tree generation algorithm [costabal2016generating]. We applied stimuli with ms, s and amplitude twice the diastolic threshold onto the Purkinje-Myocardial Junctions (PMJs) identified from the terminal nodes of the Purkinje network. Electrical impulse propagation is simulated using MCM with RPI and MKI approximants as well as FEM.

The nodal support domains in MCM simulations are constructed using the nearest-neighbor approach with the 150 nearest nodes of each field node included in its support domain. We adopt this approach due to the irregular distribution of the nodes in the tetrahedral mesh. For such a distribution, constructing dilated support domains requires using a large value of the dilatation coefficient to avoid numerical instability, which leads to a very large number of neighboring nodes. Here, we use 150 nearest nodes to accurately capture the steep voltage gradients of the monodomain model.

In Figure 4, we provide the MCM and FEM simulation results in terms of local activation time (LAT) maps for both baseline and LBBB conditions.

The mean local activation time is computed for internal field nodes (LAT) and surface field nodes (LAT) separetely. Mean LAT is found to be 23.2 ms for MCM-RPI and MCM-MKI under healthy conditions, while mean LAT for FEM is 20.9 ms. Mean LAT is found to be 27.3 ms for MCM-RPI and MCM-MKI and 21.4 ms for FEM. For LBBB conditions, mean LAT is 47.9 ms for MCM-RPI and MCM-MKI and 40.4 ms for FEM. Mean LAT is 51.9 ms for MCM-RPI and MCM-MKI and 41.1 ms for FEM. The execution time under healthy and LBBB conditions is 15.7 and 16.2 mins for MCM-RPI, 15.3 and 16.1 mins for MCM-MKI and 7.7 and 7.9 mins for FEM. For MCM-RPI and MCM-MKI methods, an additional time of 14.7 and 11.6 mins, respectively, is required for the calculation of the trial functions and their gradients.

### 4.4 Electrical propagation in a 3D biventricular geometry under myocardial infarction conditions

We evaluate the MCM method for the simulation of myocardial infarction conditions using the biventricular anatomy of section 4.3. We introduce an infarct scar at the anterior wall of the left ventricle and we assume that the conductivity of the scarred tissue is zero. For these simulations, we do not model the border zone between the infarct and non-infarct areas. We compare the solution provided by MCM-RPI and MCM-MKI with the FEM solution in terms of LAT. The MCM is solved by representing the biventricular anatomy with a tetrahedral mesh, also used for FEM, as well as by using an immersed grid model [mountris2020next].

The immersed grid (238,519 nodes) is generated by distributing equidistant nodes inside the bounding box of the biventricular anatomy’s boundary surface mesh. The nodal spacing of the immersed grid nodes is selected as twice the mean circumference of the triangles in the surface mesh, mm. A point containment test algorithm [liu2010new] is used to discard nodes that span outside of the surface’s boundary. Consequently, the final immersed grid is composed of the equidistant nodes inside the biventricular model and the nodes on the boundary surface mesh of the model (Figure 5).

We use the same tissue properties and electrical stimulation protocol as in section 4.3 to generate the LAT maps for MCM-RPI, MCM-MKI and FEM for the tetrahedral and immersed grid models of the biventricular anatomy with scarred tissue. The LAT maps are shown in Figure 6. For the tetrahedral model, the mean LAT is 29.8 ms for MCM-RPI and MCM-MKI and 21.6 for FEM. Mean LAT is 36.5 ms for MCM-RPI and MCM-MKI and 22.8 for FEM. The execution time is 16.1 mins for MCM-RPI, 15.3 mins for MCM-MKI and 7.8 mins for FEM. The required time to compute the RPI trial function and gradient is 14.5 mins and for the MKI trial function and gradient it is 11.6 mins. For the immersed grid model, the mean LAT is 25.2 ms for MCM-RPI and MCM-MKI, while the mean LAT 33.4 ms for MCM-RPI and MCM-MKI. The execution time is 13.8 mins for MCM-RPI and 13.5 mins for MCM-MKI. Additional time of 12.6 and 9.8 mins is required for the calculation of the RPI and MKI trial functions and gradients for the immersed grid model.

## 5 Concluding remarks

In this study, we derived the mixed collocation method (MCM) to solve the monodomain model for the numerical simulation of cardiac electrophysiology. We considered two different interpolating trial functions, the radial point interpolation (RPI) and the moving Kriging interpolation (MKI). We solved several numerical examples in 2D and 3D domains comparing the MCM solution with a solution obtained by FEM. The accuracy of MCM solutions was very similar for both RPI and MKI approximants. However, MKI was found to be more efficient since, in contrast to RPI, it does not require polynomial enrichment.

In all the numerical examples in 2D tissue sheets and 3D tissue slabs, good agreement was found between MCM and FEM solutions. A convergence analysis in 3D demonstrated that the MCM solution improves for larger support domains as the number of included collocation points is increased. For large scale problems, including the 150 nearest neighbours of each field node in its support domain was found to be an optimal choice balancing accuracy and memory footprint.

In large scale simulations of a biventricular anatomy model, MCM was able to produce LAT maps in good agreement with FEM under healthy and disease (LBBB, myocardial infarction) conditions. The largest differences in LAT between MCM and FEM were mainly found at the surface nodes of the model. This limitation was attributed to the negative effect of non-smooth changes of normal vectors in the imposition of Neumann boundary conditions. Introducing the immersed grid model approach, we were able to obtain improved results both in terms of accuracy and efficiency.

MCM is proved to be a promising alternative to FEM for cardiac electrophysiology simulation since its meshfree nature alleviates the need for the generation of a mesh and could thus allow fast model generation in a clinical setting.

## Acknowledgements

This work was supported by the European Research Council under grant agreement ERC-StG 638284, by Ministerio de Ciencia e Innovación (Spain) through project PID2019-105674RB-I00 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).

Comments

There are no comments yet.