# Cardiac Electrophysiology Meshfree Modeling through the Mixed Collocation
Method

Konstantinos A. Mountris
Aragón Institute of Engineering Research, IIS Aragón
CIBER-BBN
University of Zaragoza
Spain, Zaragoza, ZGZ 50018
<EMAIL_ADDRESS>
&Esther Pueyo
Aragón Institute of Engineering Research, IIS Aragón,
CIBER-BBN
University of Zaragoza
Spain, Zaragoza, ZGZ 50018
<EMAIL_ADDRESS>
[mail]<EMAIL_ADDRESS>[url] https://www.mountris.org

###### Abstract

We present the meshfree Mixed Collocation Method (MCM) to solve the monodomain
model for numerical simulation of cardiac electrophysiology. We apply MCM to
simulate cardiac electrical propagation in 2D tissue sheets and 3D tissue
slabs as well as in realistic large-scale biventricular anatomies.
Capitalizing on the meshfree property of MCM, we introduce an immersed grid
approach for automated generation of nodes in the modeled cardiac domains. We
demonstrate that MCM solutions are in agreement with FEM solutions, thus
confirming their suitability for simulation of cardiac electrophysiology both
in healthy and disease conditions, such as left-bundle-branch block (LBBB) and
myocardial infarction. Despite the fact that the computational time for MCM
calculations is longer than for FEM, its efficiency in dealing with domains
presenting irregularity, nonlinearity and discontinuity make MCM a promising
alternative for heart’s electrical investigations.

_K_ eywords meshfree $\cdot$ mixed collocation method $\cdot$ MCM $\cdot$
cardiac electrophysiology $\cdot$ monodomain model

##  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 [1, 2]. Novel in silico models of increasing complexity are
continuously being developed to simulate the electrophysiology [3] and
mechanics of the heart [4] from the cell [5] to the whole-organ level [6]. At
the tissue and organ levels, electrophysiology is simulated by using the well-
known bidomain [7] and monodomain [8] 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 [9].

Commonly, state-of-the-art simulators [10, 11] 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 [12].

Alternative meshfree methods can alleviate the mesh requirement and have been
proposed for both cardiac electrophysiological [13, 14, 15] and mechanical
[12, 16] 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 [17]. 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 [18].

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) [19]. To overcome these problems, the total Lagrangian
formulation of SPH [20] and gradient normalization [19] were applied in [13,
12] 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 [21, 22]. 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 [23], the accuracy is ameliorated in
MCM due to the interpolation of the field function’s gradient. In the seminal
work on MCM [24], the moving least squares (MLS) approximation [25] was used
as the trial function. Recently, the radial point interpolation (RPI) [26] was
proposed as an alternative to MLS in MCM [27]. 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) [28] 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:

$\begin{array}[]{ll}\partial V(t)/\partial
t=\bm{\nabla}\cdot(\bm{D\nabla}V(t))-I_{ion}(V(t))/C&\textrm{ in }\Omega\\\
\bm{n}\cdot(\bm{D\nabla}V(t))=0&\textrm{ on }\partial\Omega\end{array}$ (1)

where $\partial V/\partial t$ is the time derivative of the transmembrane
potential $V$, $I_{ion}$ is the total ionic current and $C$ is the cell
capacitance per unit surface area. $\Omega$ and $\partial\Omega$ denote the
domain of interest and its boundary, respectively, and $\bm{n}$ is the outward
unit vector normal to the boundary. $\bm{D}$ is the diffusion tensor given by:

$\bm{D}=d_{0}\left[(1-\rho)\bm{f}\otimes\bm{f}+\rho\bm{I}\right]$ (2)

where $d_{0}$ denotes the diffusion coefficient along the cardiac fibers,
$\rho\leq 1$ is the transverse-to-longitudinal conductivity ratio, $\bm{f}$
denotes the cardiac fiber direction vector, $\bm{I}$ is the identity matrix
and $\otimes$ is the tensor product operation.

The transmembrane potential $V$ 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 [29] to improve computational
efficiency. By applying operator-splitting, the decoupled monodomain equation
is given by:

$\begin{array}[]{ll}\partial V(t)/\partial t=-I_{ion}(V(t))/C&\textrm{ in
}\Omega\\\ \partial V(t)/\partial
t=\bm{\nabla}\cdot(\bm{D\nabla}V(t))&\textrm{ in }\Omega\\\
\bm{n}\cdot(\bm{D\nabla}V(t))=0&\textrm{ on }\partial\Omega\end{array}$ (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 $N$ arbitrarily distributed field
nodes. We write the diffusion term from Equation (3) for each field node $I$
in the set of $N$ nodes as:

$\partial V_{I}(t)/\partial t=\bm{\nabla}\cdot\bm{q}_{I}(t)\textrm{ in
}\Omega\\\ $ (4)

where

$\bm{q}_{I}(t)=\bm{D}_{I}\bm{\nabla}V_{I}(t)$ (5)

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

$V_{I}(t)=\sum_{i=1}^{n}\phi_{I}^{i}V_{I}^{i}(t)$ (6)
$\bm{q}_{I}(t)=\sum_{i=1}^{n}\phi_{I}^{i}\bm{q}_{I}^{i}(t)$ (7)

where $n$ is the number of field nodes in the local support domain of the node
$I$. $\phi_{I}^{i}$ is the $i^{th}$ component of the vector $\bm{\phi}_{I}$
containing the basis function of the meshfree approximation at each of the $n$
nodes in the local support domain of $I$. $V_{I}^{i}(t)$ denotes the
transmembrane potential at node $i$ of the local support domain of $I$ and
$\bm{q}_{I}^{i}(t)$ 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:

$\bm{q}_{I}(t)=\bm{D}_{I}\sum_{i=1}^{n}\bm{\nabla}\phi_{I}^{i}V_{I}^{i}(t),\quad
I=1,2,\dotso,N$ (8)

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

$\bm{q}_{I}=\bm{K}_{a,I}\bm{V}_{I},\quad I=1,2,\dotso,N.$ (9)

where $\bm{K}_{a,I}$ contains the spatial partial derivatives of the meshfree
basis function at the nodes in the support domain of node $I$ scaled by
$\bm{D}_{I}$.

By grouping the equations for all field nodes $I$, with $I=1,2,\dotso,N$, the
following equation in matrix form can be obtained:

$\bm{q}=\bm{K_{a}}\bm{V}.$ (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:

$\sum_{i=1}^{n}\phi_{I}^{i}\partial V_{I}^{i}(t)/\partial
t-\sum_{i=1}^{n}\divergence{(\bm{D}_{I}\gradient{\phi_{I}^{i}}V_{I}^{i}(t))}=0,\quad
I=1,2,\dotso,N.$ (11)

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

$\bm{M}\dot{\bm{V}}+\bm{KV}=\bm{0},\quad\bm{K}=\bm{K_{s}K_{a}}$ (12)

where $\bm{M}$ is the sparse matrix collecting the basis functions and
$\bm{K}$ denotes the stiffness matrix.

###  2.2 Boundary conditions imposition

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

$\bm{N}_{bc}\bm{q}_{bc}=\bm{0},$ (13)

where $\bm{q}_{bc}$ is the vector collecting the transmembrane potential
fluxes at $\gamma_{bc}$ nodes and $\bm{N}_{bc}$ is the matrix containing the
normal vectors given by:

$\bm{N}_{bc}=\begin{bmatrix}\bm{n}_{1}&&0\\\ &\ddots&\\\
0&&\bm{n}_{\gamma_{bc}}\end{bmatrix}$ (14)

The Neumann BCs are enforced at the $\gamma_{bc}$ nodes by multiplying
Equation (13) with the penalty factor $\alpha\bm{N}_{bc}^{T}$ and adding it to
Equation (10) to obtain:

$\bm{q}_{bc}+\alpha\bm{N}_{bc}^{T}\bm{N}_{bc}\bm{q}_{bc}=\bm{K}_{a}^{bc}\bm{V}^{bc}.$
(15)

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

$\bm{q}_{bc}=\\{\bm{I}+\alpha\bm{N}_{bc}^{T}\bm{N}_{bc}\\}^{-1}\\{\bm{K_{a}}^{bc}\bm{V}^{bc}\\}=\bm{Q}^{-1}\\{\bm{K}_{a}^{bc}\bm{V}^{bc}\\},$
(16)

where $\bm{I}$ is the identity matrix and
$\bm{Q}=\bm{I}+\alpha\bm{N}_{bc}^{T}\bm{N}_{bc}$. Combining Equations (12) and
(16), the matrix form of the monodomain model’s diffusion term is given by:

$\bm{M}\dot{\bm{V}}+\bm{K^{\prime}V}=\bm{0},$ (17)

where
$\quad\bm{K^{\prime}}=\left[\bm{K_{s}^{bc}Q^{-1}K_{a}^{bc}}\quad,\quad\bm{K_{s}^{in}K_{a}^{in}}\right]$.
Here, the superscripts $bc$ and $in$ connote the row entries of the matrices
$K_{s}$ and $K_{a}$ that correspond to the $\gamma_{bc}$ nodes on
$\partial\Omega$ and the $\gamma_{in}$ nodes in $\Omega$, respectively, such
that $\gamma_{bc}\cup\gamma_{in}=\\{1,\cdots,N\\}$. The value of the penalty
factor $\alpha$ should be sufficiently large to ensure accurate imposition of
the boundary condition, as instability issues may arise if $\alpha$ is too
large. In this study, we used $\alpha=10^{6}$ as it was found to be the
optimal value in [27].

##  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 $\bm{\phi}$. In this work, we
consider only trial functions that possess the delta Kronecker property,
namely the radial point interpolation (RPI) [26] and the moving Kriging
interpolation (MKI) [28].

###  3.1 Radial point interpolation

The RPI trial function $\bm{\phi}_{I}$ for a field node $I$ is obtained by:

$\bm{\phi}_{I}=\\{\bm{r}_{I}\;\;\bm{p}_{I}\\}\bm{G}^{-1},$ (18)

where $\bm{r}_{I}$ is the radial basis function (RBF) for node $I$. For a node
$i$ in the support domain of $I$ (with $I$ included in its support domain),
$\bm{r}_{i}$ is given by:

$\bm{r}_{i}=\left[r_{i1}\;r_{i2}\;\dotso\;r_{in}\right].$ (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:

$r_{Ii}=\left(d^{i2}_{I}+r_{c}^{2}\right)^{q}=\left[(x_{I}-x_{i})^{2}+(y_{I}-y_{i})^{2}+(z_{I}-z_{i})^{2}+r_{c}^{2}\right]^{q}$
(20)

where $d^{i}_{I}$ denotes the Euclidean distance between nodes $i$ and $I$,
$i=1,\dotso,n$ and $r_{c}$ and $q$ are positive-valued shape parameters of the
MQ-RBF. For spherical support domains, the shape parameter $r_{c}$ is given
by:

$r_{c}=\alpha_{c}d_{c}$ (21)

where $d_{c}$ denotes the radius of the support domain of node $I$ and
$\alpha_{c}$ is a dimensionless constant. RBF fail to reconstruct exactly a
linear polynomial field, therefore, the RPI is enriched with the linear
polynomial basis $\bm{p_{I}}$ to ensure $C^{1}$ continuity. For 3D problems,
$\bm{p}_{I}$ is given by:

$\bm{p}_{I}=[1\;x_{I}\;y_{I}\;z_{I}].$ (22)

Finally the matrix $\bm{G}$ is given by:

$\bm{G}=\left[\begin{matrix}\bm{R}&\bm{P}\\\
\bm{P}^{T}&\bm{0}\end{matrix}\right]_{(n+m\times n+m)},$ (23)

where $m$ is the number of components of the polynomial basis ($m=4$ for
linear $\bm{p}_{I}$ in 3D). $\bm{R}$ and $\bm{P}$ denote the RBF and
polynomial basis moment matrices:

$\bm{R}=\left[\begin{matrix}r_{11}&r_{12}&\ldots&r_{1n}\\\
r_{21}&r_{22}&\ldots&r_{2n}\\\ \ldots&\ldots&\ldots&\ldots\\\
r_{n1}&r_{n2}&\ldots&r_{nn}\\\ \end{matrix}\right]_{(n\times
n)},\;\;\bm{P}=\left[\begin{matrix}1&x_{1}&y_{1}&z_{1}\\\
1&x_{2}&y_{2}&z_{2}\\\ \ldots&\ldots&\ldots&\ldots\\\ 1&x_{n}&y_{n}&z_{n}\\\
\end{matrix}\right]_{(n\times m)}.$ (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 $C^{1}$
continuity. The trial function $\bm{\phi}$ at node $I$ is given by:

$\bm{\phi}_{I}=\bm{p}_{I}\bm{A}+\bm{c}_{I}\bm{B}$ (25)

where $\bm{p}_{I}$ is the linear polynomial basis defined in Equation (22) and
$\bm{c}_{I}$ denotes the correlation function for node $I$. For a node $i$ in
the support domain of $I$, $\bm{c}_{i}$ is given by:

$\bm{c}_{i}=\left[c_{i1}\;c_{i2}\;\dotso\;c_{in}\right]$ (26)

where $c_{ij}$ is the value of the correlation function $\bm{c}_{i}$ at the
$j^{th}$ node of the support domain. In this work, we use the MQ-RBF as the
correlation function (Equation 20). The matrices $\bm{A}$ and $\bm{B}$ are
obtained by:

$\displaystyle\begin{split}\bm{A}&=(\bm{P}^{T}\bm{C}^{-1}\bm{P})^{-1}\bm{P}^{T}\bm{C}^{-1}\\\
\bm{B}&=\bm{C}^{-1}(\bm{I}-\bm{P}\bm{A})\end{split}$ (27)

where $\bm{I}$ is the $n\times n$ identity matrix, $\bm{P}$ is the $n\times m$
moment matrix of the linear polynomial basis given by Equation (28) and
$\bm{C}$ is the $n\times n$ correlation matrix for the $n$ nodes in the
support domain of $I$ given by:

$\bm{C}=\left[\begin{matrix}c_{11}&c_{12}&\ldots&c_{1n}\\\
c_{21}&c_{22}&\ldots&c_{2n}\\\ \ldots&\ldots&\ldots&\ldots\\\
c_{n1}&c_{n2}&\ldots&c_{nn}\\\ \end{matrix}\right]_{(n\times n)}.$ (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
$\alpha_{c}=1.03$ and $q=1.42$ in 2D simulations. In 3D, they are constructed
using $\alpha_{c}=1.03$ and $q=1.82$, 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 [31] with adaptive time step $dt=0.1$ ms. Human ventricular
cellular electrophysiology is represented by the O’Hara et al. cell model
[32]. 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 $5\times 5$ 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
$\bm{f}$ is considered to be parallel to the $x$-axis. We use a diffusion
coefficient value $d_{0}=0.0013$ cm2/ms and a transverse-to-longitudinal
conductivity ratio $\rho=0.2$. Stimuli with amplitude twice the diastolic
threshold, period $t_{T}=1$ s and duration $t_{d}=1$ ms are applied on the
left side of the tissue ($x=0$ cm). The propagation of the action potential
(AP) is simulated for a total time $t_{s}=3$ 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
$h=\\{0.2,\;0.1,\;0.05,\;0.025\\}$ cm. The considered support domains in the
meshfree approximation have size $s_{d}=\alpha_{sd}h$, with $\alpha_{sd}=2.8$.
The generated APs at the center of the tissue sheet ($x=2.5$ cm, $y=2.5$ cm)
in the time interval $t=[0-3]$ 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 $h=\\{0.2,\;0.1,\;0.05,\;0.025\\}$
cm. The efficiency of each simulation is evaluated in terms of execution time
in Figure 3a.

Figure 1: Voltage traces at the center ($x=2.5$ cm, $y=2.5$ cm) of a $5\times
5$ cm ventricular tissue sheet calculated using FEM (continuous black), MCM
with MKI (dotted green) and MCM with RPI (dashed red). The nodal
discretization spacing is: a) $h=0.2$ cm, b) $h=0.1$ cm, c) $h=0.05$ cm, and
d) $h=0.025$ cm.

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

We investigate the effect of the support domain’s dilatation coefficient
$\alpha_{sd}$ by computing the normalized root mean square (NRMS) error
between MCM and FEM solutions for a $3\times 3\times 3$ cm tissue slab. The
tissue is assumed to be composed of epicardial ventricular cells. Stimuli of
amplitude twice diastolic threshold, period $t_{T}=1$ s and duration $t_{d}=1$
ms are applied onto the left side of the tissue slab ($x=0$ cm). The tissue
slab is discretized with $h=0.05$ cm and varying dilatation coefficient
$\alpha_{sd}\in=\\{2.25,\;2.50,\;2.75,\;3.00,\;3.25,\;3.50\\}$. The NRMS error
is computed using the formula:

$\textrm{NRMS}=\frac{\sqrt{\displaystyle\frac{\displaystyle\sum_{I=1}^{N}\left(V_{MCM_{I}}-V_{FEM_{I}}\right)^{2}}{N}}}{\displaystyle\max_{I}V_{FEM_{I}}-\displaystyle\min_{I}V_{FEM_{I}}}$
(29)

where $V_{MCM_{I}}$ and $V_{FEM_{I}}$ denote the transmembrane potential value
at node $I$ computed with MCM and FEM, respectively. The NRMS error
convergence plot is given in Figure 2.

Figure 2: Convergence rate for varying support domain’s dilatation coefficient
$\alpha_{sd}\in\\{2.25,\;2.50,\;2.75,\;3.00,\;3.25,\;3.50\\}$ for MCM with MKI
(continuous blue) and RPI (slashed red) trial functions. Convergence is
evaluated by computing the NRMS error with respect to a FEM simulation.

The maximum NRMS error of the MCM solutions with RPI or MKI trial functions is
obtained for $\alpha_{sd}=2.25$ and is equal to 0.262 and 0.263, respectively.
The minimum NRMS error is obtained for $\alpha_{sd}=3.5$ 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.

Figure 3: Execution time for a) varying spacing in a $5\times 5$ cm tissue
sheet and b) varying support dilatation coefficient in a $3\times 3\times 3$
cm tissue slab. Execution time is reported for FEM (black), MCM with MKI trial
functions (blue) and MCM with RPI trial functions (red).

###  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 [33, 34] and
is part of an LBBB dataset for experimental studies of cardiac
resynchronization therapy [35, 36]. We compute the direction of the myocardial
fibers using a rule-based method [37].

We define a diffusion coefficient value $d_{0}=0.002$ cm2/ms in the
longitudinal direction of the myocardial fibers and a transverse-to-
longitudinal conductivity ratio $\rho=0.2$. 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 [38].

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 [39]. We applied
stimuli with $t_{d}=1$ ms, $t_{T}=1$ 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.

Figure 4: Local activation time maps for FEM (left), MCM-RPI (center), and
MCM-MKI (right) simulations on a biventricular model. Connective tissue is
represented in yellow.

The mean local activation time is computed for internal field nodes (LATin)
and surface field nodes (LATs) separetely. Mean LATin is found to be 23.2 ms
for MCM-RPI and MCM-MKI under healthy conditions, while mean LATin for FEM is
20.9 ms. Mean LATs is found to be 27.3 ms for MCM-RPI and MCM-MKI and 21.4 ms
for FEM. For LBBB conditions, mean LATin is 47.9 ms for MCM-RPI and MCM-MKI
and 40.4 ms for FEM. Mean LATs 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 [40].

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, $h=0.9$ mm. A point
containment test algorithm [41] 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).

Figure 5: Immersed grid model generation. a) The boundary surface mesh is
extracted from the tetrahedral biventricular mesh. b) Equidistant nodes with
spacing twice the mean circumference of the triangles in the surface mesh are
distributed inside its bounding box. c) Equidistant nodes located outside the
surface mesh are discarded using a point containment test algorithm and the
final immersed grid model is given by the union of the surface mesh nodes and
the equidistant nodes enclosed by the surface.

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 LATin is 29.8 ms for MCM-RPI and MCM-MKI and 21.6 for FEM. Mean LATs 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 LATin is 25.2 ms for MCM-RPI and MCM-MKI, while the mean LATs 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.

Figure 6: Local activation time maps for FEM (left), MCM-RPI (center), and
MCM-MKI (right) simulations on an immersed grid of the biventricular model.
Non-conductive connective tissue is represented in yellow.

##  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).

## References

  * [1] Alejandro Lopez-Perez, Rafael Sebastian, and Jose M Ferrero.  Three-dimensional cardiac computational modelling: methods, features and applications.  Biomedical engineering online, 14(1):35, 2015.
  * [2] Steven A Niederer, Joost Lumens, and Natalia A Trayanova.  Computational models in cardiology.  Nature reviews cardiology, 16(2):100–111, 2019.
  * [3] David Adolfo Sampedro-Puente, Fabien Raphel, Jesus Fernandez-Bes, Pablo Laguna, Damiano Lombardi, and Esther Pueyo.  Characterization of spatio-temporal cardiac action potential variability at baseline and under $\beta$-adrenergic stimulation by combined unscented kalman filter and double greedy dimension reduction.  IEEE Journal of Biomedical and Health Informatics, 2020.
  * [4] Ritesh R Rama and Sebastian Skatulla.  Towards real-time cardiac mechanics modelling with patient-specific heart anatomies.  Computer methods in applied mechanics and engineering, 328:47–74, 2018.
  * [5] Esther Pueyo, CE Dangerfield, OJ Britton, László Virág, Kornél Kistamás, Norbert Szentandrássy, Norbert Jost, András Varró, Péter Pál Nánási, Kevin Burrage, et al.  Experimentally-based computational investigation into beat-to-beat variability in ventricular repolarization and its response to ionic current inhibition.  PloS one, 11(3):e0151461, 2016.
  * [6] Radomir Chabiniok, Vicky Y Wang, Myrianthi Hadjicharalambous, Liya Asner, Jack Lee, Maxime Sermesant, Ellen Kuhl, Alistair A Young, Philippe Moireau, Martyn P Nash, et al.  Multiphysics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics.  Interface focus, 6(2):20150083, 2016.
  * [7] Leslie Tung.  A bi-domain model for describing ischemic myocardial DC potentials.  PhD thesis, Massachusetts Institute of Technology, 1978.
  * [8] James P Keener and James Sneyd.  Mathematical physiology: systems physiology. II.  Springer, 2009.
  * [9] Mark Potse, Bruno Dubé, Jacques Richer, Alain Vinet, and Ramesh M Gulrajani.  A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart.  IEEE transactions on biomedical engineering, 53(12):2425–2435, 2006\.
  * [10] Gary R Mirams, Christopher J Arthurs, Miguel O Bernabeu, Rafel Bordas, Jonathan Cooper, Alberto Corrias, Yohan Davit, Sara-Jane Dunn, Alexander G Fletcher, Daniel G Harvey, et al.  Chaste: an open source c++ library for computational physiology and biology.  PLoS computational biology, 9(3):e1002970, 2013.
  * [11] E. J. Vigmond, M. Hughes, G. Plank, and L. J. Leon.  Computational tools for modeling electrical activity in cardiac tissue.  Journal of electrocardiology, 36 Suppl:69–74, 2003.
  * [12] Èric Lluch, Mathieu De Craene, Bart Bijnens, Maxime Sermesant, Jérôme Noailly, Oscar Camara, and Hernán G Morales.  Breaking the state of the heart: meshless model for cardiac mechanics.  Biomechanics and modeling in mechanobiology, pages 1–13, 2019.
  * [13] Eric Lluch, Rubén Doste, Sophie Giffard-Roisin, Alexandre This, Maxime Sermesant, Oscar Camara, Mathieu De Craene, and Hernán G Morales.  Smoothed particle hydrodynamics for electrophysiological modeling: an alternative to finite element methods.  In International conference on functional imaging and modeling of the heart, pages 333–343. Springer, 2017.
  * [14] Heye Zhang, Huajun Ye, and Wenhua Huang.  A meshfree method for simulating myocardial electrical activity.  Computational and mathematical methods in medicine, 2012, 2012.
  * [15] Konstantinos A Mountris, Carlos Sanchez, and Esther Pueyo.  A novel paradigm for in silico simulation of cardiac electrophysiology through the mixed collocation meshless petrov-galerkin method.  In 2019 Computing in cardiology (CinC), pages Page–1. IEEE, 2019\.
  * [16] D Legner, S Skatulla, J MBewu, RR Rama, BD Reddy, C Sansour, NH Davies, and T Franz.  Studying the influence of hydrogel injections into the infarcted left ventricle using the element-free galerkin method.  International journal for numerical methods in biomedical engineering, 30(3):416–429, 2014.
  * [17] Ted Belytschko, Yun Yun Lu, and Lei Gu.  Element-free galerkin methods.  International journal for numerical methods in engineering, 37(2):229–256, 1994.
  * [18] Konstantinos A Mountris, George C Bourantas, Daniel Millán, Grand R Joldes, Karol Miller, Esther Pueyo, and Adam Wittek.  Cell-based maximum entropy approximants for three-dimensional domains: Application in large strain elastodynamics using the meshless total lagrangian explicit dynamics method.  International journal for numerical methods in engineering, 121(3):477–491, 2020.
  * [19] JK Chen, JE Beraun, and TC Carney.  A corrective smoothed particle method for boundary value problems in heat conduction.  International journal for numerical methods in engineering, 46(2):231–252, 1999.
  * [20] Javier Bonet and T-SL Lok.  Variational and momentum preservation aspects of smooth particle hydrodynamic formulations.  Computer methods in applied mechanics and engineering, 180(1-2):97–115, 1999.
  * [21] Satya N. Atluri.  The meshless method (MLPG) for domain and BIE discretizations.  Tech Science Press, 2004.
  * [22] Satya N. Atluri and Shengping Shen.  The basis of meshless domain discretization: the meshless local Petrov-Galerkin (MLPG) method.  Advances in computational mathematics, 23(1-2):73–93, jul 2005\.
  * [23] Nicolas Libre, Arezoo Emdadi, Edward Kansa, Mohammad Rahimian, and Mohammad Shekarchi.  A stabilized rbf collocation scheme for neumann type boundary value problems.  Computer modeling in engineering & sciences, 24(1):61–80, jan 2008.
  * [24] Satya N. Atluri, Hugh T. Liu, and Zdravko Dovedan Han.  Meshless local petrov-galerkin (mlpg) mixed collocation method for elasticity problems.  Computer modeling in engineering & sciences, 4(3):141, 2006.
  * [25] P Lancaster and K Salkauskas.  Surfaces generated by moving least squares methods.  Mathematics of computation, 37(155):141–158, 1981.
  * [26] JG Wang and GRs Liu.  A point interpolation meshless method based on radial basis functions.  International journal for numerical methods in engineering, 54(11):1623–1648, 2002.
  * [27] Konstantinos A Mountris and Esther Pueyo.  The radial point interpolation mixed collocation method for the solution of transient diffusion problems.  Engineering analysis with boundary elements, 121:207–216, 2020\.
  * [28] L Gu.  Moving kriging interpolation and element-free galerkin method.  International journal for numerical methods in engineering, 56(1):1–11, 2003.
  * [29] Zhilin Qu and Alan Garfinkel.  An advanced algorithm for solving partial differential equation in cardiac conduction.  IEEE transactions on biomedical engineering, 46(9):1166–1168, 1999\.
  * [30] Johannes TB Overvelde.  The moving node approach in topology optimization.  Master’s thesis, TU Delft, Delft University of Technology, 2012.
  * [31] Konstantinos A Mountris and Esther Pueyo.  A dual adaptive explicit time integration algorithm for efficiently solving the cardiac monodomain equation.  arXiv preprint arXiv:2011.04747, 2020.
  * [32] Thomas O’Hara, László Virág, András Varró, and Yoram Rudy.  Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation.  PLoS computational biology, 7(5), 2011.
  * [33] Oscar Camara.  Best (and worst) practices for organizing a challenge on cardiac biophysical models during ai summer: the crt-epiggy19 challenge.  In International workshop on statistical atlases and computational models of the heart, pages 329–341. Springer, 2019.
  * [34] Mihaela Pop, Maxime Sermesant, Oscar Camara Rey Camara, Xiahai Zhuang, Shuo Li, Alistair Young, Tommaso Mansi, and Avan Suinesiaputra.  Statistical atlases and computational models of the heart: multi-sequence CMR segmentation, CRT-EPiggy and LV full quantification challenges: 10th international workshop, STACOM 2019, held in conjunction with MICCAI 2019, Shenzhen, China, October 13, 2019, revised selected papers, volume 12009.  Springer Nature, 2020.
  * [35] Montserrat Rigol, Núria Solanes, Juan Fernandez-Armenta, Etelvino Silva, Adelina Doltra, Nicolas Duchateau, Aina Barcelo, Luigi Gabrielli, Bart Bijnens, Antonio Berruezo, et al.  Development of a swine model of left bundle branch block for experimental studies of cardiac resynchronization therapy.  Journal of cardiovascular translational research, 6(4):616–622, 2013.
  * [36] David Soto Iglesias, Nicolas Duchateau, Constantine Butakoff Kostantyn Butakov, David Andreu, Juan Fernández-Armenta, Bart Bijnens, Antonio Berruezo, Marta Sitges, and Oscar Camara.  Quantitative analysis of electro-anatomical maps: application to an experimental model of left bundle branch block/cardiac resynchronization therapy.  IEEE journal of translational engineering in health and medicine, 5:1–15, 2016.
  * [37] Ruben Doste, David Soto-Iglesias, Gabriel Bernardino, Alejandro Alcaine, Rafael Sebastian, Sophie Giffard-Roisin, Maxime Sermesant, Antonio Berruezo, Damian Sanchez-Quintana, and Oscar Camara.  A rule-based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts.  International journal for numerical methods in biomedical engineering, 35(4):e3185, 2019.
  * [38] K Andrew MacCannell, Hojjat Bazzazi, Lisa Chilton, Yoshiyuki Shibukawa, Robert B Clark, and Wayne R Giles.  A mathematical model of electrotonic interactions between ventricular myocytes and fibroblasts.  Biophysical journal, 92(11):4121–4132, 2007.
  * [39] Francisco Sahli Costabal, Daniel E Hurtado, and Ellen Kuhl.  Generating purkinje networks in the human heart.  Journal of biomechanics, 49(12):2455–2465, 2016.
  * [40] Konstantinos A Mountris and Esther Pueyo.  Next-generation in silico cardiac electrophysiology through immersed grid meshfree modeling: Application to simulation of myocardial infarction.  In 2020 Computing in Cardiology, pages 1–4. IEEE, 2020.
  * [41] Jianfei Liu, YQ Chen, José M Maisog, and George Luta.  A new point containment test algorithm based on preprocessing and determining triangles.  Computer-aided design, 42(12):1143–1150, 2010.