1Laboratory for Mathematics in Biology and Medicine, Department of
Mathematics, Università di Trento, via Sommarive 14, 38123, Trento, Italy;
1Faculty of Mathematics and Informatics, FernUni Schweiz, Schinerstrasse 18,
3900, Brig, Switzerland G. Rosilho de Souza et al.

# Explicit stabilized multirate methods for the monodomain model in cardiac
electrophysiology

Giacomo Rosilho de Souza<EMAIL_ADDRESS>Euler
Institute, Università della Svizzera italiana, via G. Buffi 13, 6900 Lugano,
Switzerland;
,  Marcus J. Grote<EMAIL_ADDRESS>Department of Mathematics,
University of Basel, Rheinsprung 21, 4051 Basel, Switzerland;
,  Simone Pezzuto<EMAIL_ADDRESS>and  Rolf Krause

(Date: dates will be set by the publisher)

###### Abstract.

Fully explicit stabilized multirate (mRKC) methods are well-suited for the
numerical solution of large multiscale systems of stiff ordinary differential
equations thanks to their improved stability properties. To demonstrate their
efficiency for the numerical solution of stiff, multiscale, nonlinear
parabolic PDE’s, we apply mRKC methods to the monodomain equation from cardiac
electrophysiology. In doing so, we propose an improved version, specifically
tailored to the monodomain model, which leads to the explicit exponential
multirate stabilized (emRKC) method. Several numerical experiments are
conducted to evaluate the efficiency of both mRKC and emRKC, while taking into
account different finite element meshes (structured and unstructured) and
realistic ionic models. The new emRKC method typically outperforms a standard
implicit-explicit baseline method for cardiac electrophysiology. Code
profiling and strong scalability results further demonstrate that emRKC is
faster and inherently parallel without sacrificing accuracy.

###### Key words and phrases:

multirate explicit stabilized methods, Rush–Larsen, electrophysiology,
monodomain model, ionic model, local time-stepping

###### 1991 Mathematics Subject Classification:

65L04, 65L06, 65L10, 65L20

This work was supported by the European High-Performance Computing Joint
Undertaking EuroHPC under grant agreements No 955495 (MICROCARD) and No 955701
(TIME-X), co-funded by the Horizon 2020 programme of the European Union (EU)
and the Swiss State Secretariat for Education, Research and Innovation. We
also acknowledge the CSCS-Swiss National Supercomputing Centre (project no.
s1074)

Les méthodes multirate stabilisées explicites (mRKC) sont bien adaptées à la
résolution numérique de grands systèmes multi-échelles d’équations
différentielles ordinaires rigides grâce à leur meilleure stabilité. Afin de
démontrer leur efficacité pour la résolution numérique d’EDP paraboliques non-
linéaires raides et multi-échelles, nous appliquons les méthodes mRKC à
l’équation monodomaine de l’électrophysiologie cardiaque. Ce faisant, nous
proposons une version améliorée, spécifiquement adaptée au modèle monodomaine,
qui conduit à la méthode exponentielle multirate stabilisée (emRKC) explicite.
Plusieurs expériences numériques permettent d’évaluer l’efficacité des
méthodes mRKC et emRKC, tout en tenant compte de différents maillages
d’éléments finis (structurés et non-structurés) et des modèles ioniques
réalistes. La nouvelle méthode emRKC surpasse généralement une méthode de base
implicite-explicite standard pour l’électrophysiologie cardiaque. Le profilage
du code et les résultats de scalabilité démontrent en outre que emRKC est plus
rapide et intrinsèquement parallèle sans sacrifier la précision.

##  1\. Introduction

Spatial discretizations of parabolic partial differential equations (PDEs),
possibly nonlinear, typically lead to large systems of stiff ordinary
differential equations (ODEs). The degree of stiffness, be it mild or severe,
depends in particular on the rate of change of individual components or the
mesh size. For time discretization, there is generally a trade-off between
explicit and implicit time-marching schemes. Explicit schemes are easier to
implement and also cheaper per time-step, but require more steps due to
stability constraints. In contrast, implicit schemes require sophisticated
(parallel) solvers for the needed linear or nonlinear systems and thus are
more expensive per time-step, but they also require fewer time-steps. Choosing
between explicit and implicit time integration is a delicate process which
depends on several factors, such as the stiffness of the system or the
efficiency of the linear or nonlinear solvers available. Indeed, an
unfortunate choice can increase the computational cost by orders of magnitude
due to abrupt changes in any methods’ performance.

Explicit stabilized methods, such as Runge–Kutta–Chebyshev (RKC) methods [63,
59], strike a balance between explicit and implicit schemes while their
computational cost smoothly depends on the degree of stiffness. On the one
hand, they are fully explicit and thus simple to implement in parallel, even
in the presence of nonlinearity; on the other hand, they adapt to the current
degree of stiffness by modifying the number of stages accordingly. Since their
stability domain grows quadratically with the number of stages, they are
significantly more efficient than standard explicit methods. Compared to
implicit methods, they can also be more efficient in particular for large-
scale three-dimensional linear, or nonlinear, problems[4, 1, 5, 18]. In
summary, explicit stabilized schemes combine the advantages of both explicit
and implicit methods while avoiding any abrupt performance deterioration when
transitioning from a non-stiff to stiff regime, or from linear to nonlinear
problems.

When stiffness is induced only by a few components, as in the presence of
spatial local mesh refinement or coupled chemical reactions, the efficiency of
RKC methods can deteriorate. To overcome that crippling effect due to a few
severely stiff components, multirate Runge–Kutta–Chebyshev (mRKC) methods were
proposed in [3]; they remain fully explicit, and thus easy to parallelize,
though with a stability condition independent of those few severely stiff
components. Explicit multirate RKC methods thus permit to circumvent any
overly stringent stability constraint due to local mesh refinement; hence,
they also extend local time-stepping (LTS) methods [27, 29, 30, 28] to
parabolic PDE’s.

Here we shall consider the monodomain model for cardiac electrophysiology,
which consists of a parabolic PDE for the diffusion of the electric potential
and a system of ODE’s for the ionic model of the intracellular and membrane
dynamics of the cardiac cell [13]. After spatial discretization, the
monodomain model can be written as

$y^{\prime}=f_{F}(y)+f_{S}(y)+f_{E}(y),\qquad y(0)=y_{0},$ (1)

where $f_{F}$ is cheap but stiff (“F” for “fast”), $f_{S}$ is non-stiff but
expensive (“S” for “slow”), and $f_{E}$ is in general a severely stiff term of
the form

$f_{E}(y)=\Lambda(y)(y-y_{\infty}(y)),$ (2)

with $\Lambda(y)$ a diagonal matrix and $y_{\infty}(y)$ a nonlinear function.
Typically, $f_{F}$ represents the discrete diffusion operator while $f_{S}$,
$f_{E}$ depend on the ionic model; here, $f_{E}$ corresponds to the gating
variables and $f_{S}$ to the remaining ones.

A common and efficient approach for the time integration of the non-diffusive
component of 1 is based on the Rush–Larsen scheme [53, 14]. The Rush–Larsen
scheme is a first-order splitting method employing explicit Euler for $f_{S}$
and exponential Euler for $f_{E}$ (hence the “E” notation for the
“exponential”)– see [39] for a recent review paper on the possible variants of
this approach. The diffusive term $f_{F}$ is treated either explicitly or
implicitly, depending on the degree of stiffness. For the monodomain model,
$f_{F}$ is mildly stiff and, as a consequence, the choice between an explicit
or an implicit approach somewhat challenging. Clearly, explicit stabilized
schemes, such as RKC methods, can efficiently treat the $f_{S}$ and $f_{F}$
terms, but often perform poorly for the full monodomain model, due to its
multiple components and inherent multiscale nature.

In contrast, multirate explicit stabilized methods such as mRKC [3, 16] are
more appropriate for the monodomain model because they can adapt the number of
stages locally to individual components depending on their respective degree
of stiffness. Here following [53, 14], we shall further take full advantage of
the special structure 2 of $f_{E}$ by integrating it with an exponential
method, while integrating $f_{F}$ and $f_{S}$ explicitly with the mRKC scheme.
We denote this combined exponential–mRKC time integrator by _exponential
multirate RKC_ (emRKC). The emRKC method is cheap because it is explicit while
utilizing an exponential integrator for the most severely stiff yet diagonal
terms. Albeit only first-order, it is more accurate than other commonly used
first-order schemes, as its special structure, inherited from mRKC, guarantees
a robust coupling between the different model components. Although the emRKC
method was developed specifically for the monodomain equation (for both
ventricular or atrial geometries, or even the entire heart), it can be applied
to any physical model that can be written as 1 and 2.

In Section 2, we first present the monodomain model from cardiac
electrophysiology. Next, in Section 3, we recall the RKC [63, 59] and mRKC
methods from [3, 16], summarize their properties, and describe them from an
algorithmic point of view. In Section 4, we first adapt mRKC to the particular
structure 1 of the monodomain model, which yields the new emRKC method; then
we analyze its accuracy and stability properties. In Section 5, we compare
both mRKC and emRKC to a standard baseline method on a sequence of numerical
experiments. Both two- and three-dimensional problems are considered, either
with structured or unstructured meshes, including realistic ionic models.
Concluding remarks are presented in Section 6.

##  2\. Monodomain model

The monodomain model is commonly used to simulate electric potential
propagation in cardiac tissue. Though a simplification of the more complex
bidomain model, it remains important in practice. In this section, we provide
a brief introduction to the model; for a more complete description, we refer
to [10]. The monodomain model is given by the following reaction-diffusion
system of equations:

$\displaystyle C_{m}\frac{\partial V_{m}}{\partial t}$
$\displaystyle=\chi^{-1}\nabla\cdot\bigl{(}\mathbf{D}\nabla
V_{m}\bigr{)}+I_{\text{stim}}(t,\bm{x})-I_{\text{ion}}(V_{m},\bm{z}_{E},\bm{z}_{S}),\qquad$
in $\Omega\times(0,T]$, (3a) $\displaystyle\frac{\partial\bm{z}_{E}}{\partial
t}$
$\displaystyle=\bm{g}_{E}(V_{m},\bm{z}_{E})=\bm{\alpha}(V_{m})\bigl{(}1-\bm{z}_{E}\bigr{)}-\bm{\beta}(V_{m})\bm{z}_{E},$
in $\Omega\times(0,T]$, (3b) $\displaystyle\frac{\partial\bm{z}_{S}}{\partial
t}$ $\displaystyle=\bm{g}_{S}(V_{m},\bm{z}_{E},\bm{z}_{S}),$ in
$\Omega\times(0,T]$, (3c) $\displaystyle-\mathbf{D}\nabla V_{m}\cdot\bm{n}$
$\displaystyle=0,$ $\displaystyle\text{on }\partial\Omega\times(0,T],$ (3d)
$\displaystyle V_{m}(0,\bm{x})$
$\displaystyle=V_{m,0},\quad{\bm{z}}_{E}(0,\bm{x})={\bm{z}}_{E,0},\quad{\bm{z}}_{S}(0,\bm{x})={\bm{z}}_{S,0},$
$\displaystyle\text{in $\Omega$},$ (3e)

where parameters and variables are described in Table 1. Electric conduction
in the myocardium is anisotropic, with a conductivity tensor of the form

$\mathbf{D}(\bm{x})=\sigma_{\ell}\bm{a}_{\ell}(\bm{x})\otimes\bm{a}_{\ell}(\bm{x})+\sigma_{t}\bm{a}_{t}(\bm{x})\otimes\bm{a}_{t}(\bm{x})+\sigma_{n}\Bigl{(}\mathbf{I}-\bm{a}_{\ell}(\bm{x})\otimes\bm{a}_{\ell}(\bm{x})-\bm{a}_{t}(\bm{x})\otimes\bm{a}_{t}(\bm{x})\Bigr{)},$
(4)

where
$\sigma_{\star}=\sigma_{\star}^{i}\sigma_{\star}^{e}/(\sigma_{\star}^{i}+\sigma_{\star}^{e})$
and $\sigma^{i}_{\star}$, $\sigma^{e}_{\star}$ are also given in Table 1. The
vector fields $\bm{a}_{\ell}(\bm{x})$ and $\bm{a}_{t}(\bm{x})$ are orthonormal
and are aligned with the local fiber and sheet direction, respectively.

The ionic current $I_{\text{ion}}$ models the total transmembrane currents
through a set of gating and auxiliary variables, respectively denoted by
${\bm{z}}_{E}$ and ${\bm{z}}_{S}$. The gating variables ${\bm{z}}_{E}$ are
probabilities that control the opening and closing of ion channels in the cell
membrane, whereas the remaining variables ${\bm{z}}_{S}$ represent other
quantities, such as ion concentration, or intracellular calcium dynamic. Many
different ionic models exist in the literature, depending on the type of cell
and the pathophysiological state, but in most cases the overall structure is
fixed as in 3 and follows the so-called Hodgkin–Huxley formalism.
Specifically, the total current $I_{\text{ion}}$ is the sum of several ionic
currents, each depending on its own set of gating and auxiliary variables. An
archetypal instance is the celebrated Hodgkin–Huxley model for electric
propagation in axons:

$I_{\text{ion}}(V_{m},{\bm{z}}_{E})=g_{\text{Na}}m^{3}h(V_{m}-V_{\text{Na}})+g_{\text{K}}n^{4}(V_{m}-V_{\text{K}})+g_{\text{leak}}(V_{m}-V_{\text{leak}}),$

where $\bm{z}_{E}=(m,h,n)$ and there is no auxiliary variable ${\bm{z}}_{S}$.
State-of-the-art cardiac ionic models are more advanced in the sense that they
include several more currents and the intracellular calcium dynamic, thus
significantly increasing the computational complexity and temporal stiffness.
In this respect, it is common to take advantage of the Hodgkin–Huxley
structure of the system to avoid implicit numerical schemes or excessively
small time steps. In fact, 3b can be rewritten as:

$\frac{\partial\bm{z}_{E}}{\partial
t}=\bm{\Lambda}_{E}(V_{m})\bigl{(}\bm{z}_{E}-{\bm{z}}_{\infty}(V_{m})\bigr{)},$
(5)

where $\bm{\Lambda}_{E}(V_{m})=-(\bm{\alpha}(V_{m})+\bm{\beta}(V_{m}))$ is a
diagonal matrix and
${\bm{z}}_{E,\infty}(V_{m})=-\bm{\alpha}(V_{m})(\bm{\alpha}(V_{m})+\bm{\beta}(V_{m}))^{-1}$
is the steady state value of ${\bm{z}}_{E}$. In general,
$\bm{\Lambda}_{E}(V_{m})$ is very stiff for a physiological range of $V_{m}$,
but it can be efficiently integrated with exponential methods (hence the “E”
notation) since $\bm{\Lambda}_{E}(V_{m})$ is diagonal. In contrast,
$\bm{g}_{S}(V_{m},\bm{z}_{E},\bm{z}_{S})$ is not stiff (“S” for “slow”) and
rather expensive to evaluate; thus, it is usually integrated with the explicit
Euler scheme. This split approach of integrating ${\bm{z}}_{E}$ and
${\bm{z}}_{S}$ with an exponential and standard explicit method, respectively,
was first proposed in the seminal paper of Rush and Larsen [53] and remains
very popular even to this day.

Variable | Description | Value | Units
---|---|---|---
$V_{m}(t,\bm{x})$ | Transmembrane potential | — | $\mathrm{mV}$
$\bm{z}_{E}(t,\bm{x})$ | Vector of gating variables | — | Probability
$\bm{z}_{S}(t,\bm{x})$ | Vector of auxiliary ionic variables | — | Depends
$I_{\text{ion}}(V_{m},\bm{z}_{E},\bm{z}_{S})$ | Transmembrane current density | — | $\mathrm{\SIUnitSymbolMicro A}\text{\,}{\mathrm{mm}}^{-2}$
$I_{\text{stim}}(t,\bm{x})$ | Stimulus current density | — | $\mathrm{\SIUnitSymbolMicro A}\text{\,}{\mathrm{mm}}^{-2}$
$\chi$ | surface area-to-volume ratio of cardiac cells | $140$ | ${\mathrm{mm}}^{-1}$
$C_{m}$ | cell membrane electrical capacitance | $0.01$ | $\mathrm{\SIUnitSymbolMicro F}\text{\,}{\mathrm{mm}}^{-2}$
$\sigma_{\ell}^{i}$, $\sigma_{t}^{i}$, $\sigma_{n}^{i}$ | intracellular electric conductivities | $[0.17,0.019,0.019]$ | $\mathrm{mS}\text{\,}{\mathrm{mm}}^{-1}$
$\sigma_{\ell}^{e}$, $\sigma_{t}^{e}$, $\sigma_{n}^{e}$ | extracellular electric conductivities | $[0.62,0.24,0.24]$ | $\mathrm{mS}\text{\,}{\mathrm{mm}}^{-1}$
$\bm{\alpha}(V_{m})$, $\bm{\beta}(V_{m})$ | transition rates (diagonal matrices) | — | ${\mathrm{ms}}^{-1}$
Table 1. Description of variables and parameters for the monodomain model. The
values of the coefficients are from [42].

###  2.1. Space discretization

The spatial discretization of the monodomain equations 3 using conforming
piecewise linear finite elements leads to the semi-discrete system

$\displaystyle
C_{m}\mathbf{M}\frac{\operatorname{d\\!}\bm{V}_{m}}{\operatorname{d\\!}t}$
$\displaystyle=\chi^{-1}\mathbf{A}\bm{V}_{m}+\mathbf{M}\bm{I}_{\text{stim}}(t)-\mathbf{M}I_{\text{ion}}(t,\bm{V}_{m},\bm{z}_{E},\bm{z}_{S}),$
(6a) $\displaystyle\frac{\operatorname{d\\!}\bm{z}_{E}}{\operatorname{d\\!}t}$
$\displaystyle=g_{E}(\bm{V}_{m},\bm{z}_{E}),$ (6b)
$\displaystyle\frac{\operatorname{d\\!}\bm{z}_{S}}{\operatorname{d\\!}t}$
$\displaystyle=g_{S}(\bm{V}_{m},\bm{z}_{E},\bm{z}_{S}),$ (6c)

with appropriate initial conditions and where $\mathbf{A}$ and $\mathbf{M}$
are respectively the stiffness and mass matrices. (For the sake of simplicity,
we keep using the vector notation for ${\bm{z}}_{S}$ and ${\bm{z}}_{E}$.)
Equations 6b and 6c are evaluated point-wise, due to the absence of
differential operators in space in the equations. The nonlinear term
$I_{\text{ion}}$ is evaluated component-wise and then interpolated, whereas
$\bm{I}_{\text{stim}}(t)$ is obtained by linear interpolation of
$I_{\text{stim}}(t,\bm{x})$. This is commonly done in the cardiac modeling
community and called Ionic Current Interpolation (ICI), see [44, 45].
Alternatively, one could first evaluate $\bm{V}_{m}$, $\bm{z}_{E}$,
$\bm{z}_{S}$ at the quadrature nodes and then evaluate $I_{\text{ion}}$ there,
an approach known as State Variable Interpolation (SVI). Here, we opt for ICI
due to its lower computational cost, as SVI requires interpolating every state
variable at the quadrature nodes, which is expensive for ionic models with a
large number $N$ of state variables. Finally, the mass matrix $\mathbf{M}$ is
always lumped to take advantage of explicit time integrators.

The semi-discrete equation 6 can be written as 1 for
$y=(\bm{V}_{m},\bm{z}_{E},\bm{z}_{S})^{\top}$ by defining

$\displaystyle f_{F}(t,y)=\begin{pmatrix}(\chi
C_{m}\mathbf{M})^{-1}\mathbf{A}\bm{V}_{m}\\\ \bm{0}\\\
\bm{0}\end{pmatrix},\qquad
f_{S}(t,y)=\begin{pmatrix}C_{m}^{-1}(\bm{I}_{stim}(t)-I_{ion}(t,\bm{V}_{m},\bm{z}_{E},\bm{z}_{S}))\\\
\bm{0}\\\ g_{S}(\bm{V}_{s},\bm{z}_{E},\bm{z}_{S})\end{pmatrix},$ (7a)
$\displaystyle f_{E}(t,y)=\begin{pmatrix}\bm{0}\\\
g_{E}(\bm{V}_{m},\bm{z}_{E})\\\
\bm{0}\end{pmatrix}=\Lambda(y)(y-y_{\infty}(y)),$ (7b) with
$\Lambda(y)=\begin{pmatrix}0&0&0\\\ 0&\Lambda_{E}(\bm{z}_{E})&0\\\
0&0&0\end{pmatrix},\quad y_{\infty}(y)=\begin{pmatrix}\bm{0}\\\
\bm{z}_{E,\infty}(\bm{V}_{m})\\\ \bm{0}\end{pmatrix}.$ (7c)

Note that $f_{F}$ is stiff but very cheap to evaluate, since it only involves
a matrix-vector multiplication thanks to mass-lumping. The $f_{S}$ term, which
contains the most complicated terms of the ionic model, is expensive yet non-
stiff since the stiff terms of the ionic model are contained in $g_{E}$.
Clearly, $f_{E}$ is very stiff but straightforward to integrate exponentially,
since $\Lambda(y)$ is diagonal.

The time stepping methods presented in the sequel are not limited to the
present spatial discretization. Still, the resulting mass matrix must be
diagonal, or very cheap to invert, regardless of the spatial discretization
used, to preserve the efficiency of any explicit method. For standard
conforming triangular or tetrahedral finite elements, mass-lumping is
available up to a certain polynomial order [12, 24]. Alternative spatial
discretizations that lead to diagonal or block diagonal mass matrices at any
order are, for instance, higher order finite elements with orthogonal basis
functions, spectral finite element methods [6], finite differences, and
discontinuous Galerkin methods [43].

##  3\. Multirate Explicit Stabilized Methods

In this section, we recall the multirate explicit stabilized mRKC method from
[3] for

$y^{\prime}=f(y)=f_{F}(y)+f_{S}(y),\qquad y(0)=y_{0},$ (8)

where $f_{F}$ is fast (stiff) but cheap, while $f_{S}$ slow (non-stiff or
mildly stiff) yet expensive to evaluate.

To simplify the notation, we restrict the discussion to autonomous problems.
Still, the methods proposed here also apply to non-autonomous problems while
all listed pseudo-codes apply to the general non-autonomous case. For the sake
of brevity, we concentrate here on algorithmic aspects of the mRKC method from
a practical point of view. First, in Section 3.1, we recall standard explicit
stabilized Runge–Kutta–Chebyshev (RKC) methods without any multirate time-
stepping strategy. Then, in Section 3.2, we summarize the mathematical
derivation and main properties of the mRKC method. The full mRKC Algorithm is
then listed in Section 3.3. Finally, in Section 3.4 we apply the mRKC method
to the monodomain model (3).

###  3.1. Explicit stabilized methods

The coefficients of standard explicit Runge–Kutta (RK) methods for

$y^{\prime}=f(y),\qquad y(0)=y_{0},$ (9)

are generally optimized for high accuracy. Those high-order methods, however,
suffer from severe stability constraints on the time-step. In contrast,
explicit stabilized (ES) methods fix the order $p$ (usually $p=1,2$) while
utilizing an increasing number of stages $s\geq p$ and the corresponding free
coefficients to optimize for stability. The resulting stability polynomial,
$R_{s}(z)$, of ES methods then satisfies $|R_{s}(z)|\leq 1$ for $z\in[-\beta
s^{2},0]$, where $s$ is the number of stages and $\beta>0$ depends on the
method. Therefore, the stability domain grows quadratically with $s$, and
hence with the number of function evaluations, so that ES methods remain
efficient without sacrificing explicitness. Moreover, the number of stages $s$
can easily be adapted every time-step to accommodate the problem’s
instantaneous stiffness, measured in terms of the spectral radius $\rho$ of
the right-hand side’s Jacobian. Since ES methods are defined via a recursive
relation, they can be implemented using only three vectors, regardless of $s$,
and are therefore very efficient memory-wise.

Various explicit stabilized methods, such as RKC [59], RKL [41], RKU [50] and
ROCK [4, 1], are available. Both first- and second-order versions exist for
RKC, RKL, and RKU methods, while ROCK methods are second- or fourth-order
accurate. Here, we shall only consider the classical first-order RKC method
from Section 3.1, see [2, 63, 59, 64, 65, 66] for further details.

1: function RKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f$, $\rho$) function
RKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f$, $\rho$)

2: function RKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f$, $\rho$)
$\varepsilon=0.05$ $\varepsilon=0.05$

3: $\varepsilon=0.05$ $s,\ell_{s}=\ $get_stages($\Delta t$, $\rho$,
$\varepsilon$) $s,\ell_{s}=\ $get_stages($\Delta t$, $\rho$, $\varepsilon$)

4: $s,\ell_{s}=\ $get_stages($\Delta t$, $\rho$, $\varepsilon$) $y_{n+1}=\
$RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$, $f$, $s$, $\varepsilon$)
$y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$, $f$, $s$,
$\varepsilon$)

5: $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$, $f$, $s$,
$\varepsilon$) return $y_{n+1}$ return $y_{n+1}$ return $y_{n+1}$

6: return $y_{n+1}$ return $y_{n+1}$ function get_stages($\Delta t$, $\rho$,
$\varepsilon$) return $y_{n+1}$ return $y_{n+1}$ function get_stages($\Delta
t$, $\rho$, $\varepsilon$)

7: return $y_{n+1}$ return $y_{n+1}$ function get_stages($\Delta t$, $\rho$,
$\varepsilon$) $\beta=2-4\varepsilon/3$ $\beta=2-4\varepsilon/3$

8: $\beta=2-4\varepsilon/3$ $s=\lceil\sqrt{\Delta t\,\rho/\beta}\rceil$,
$\ell_{s}=\beta s^{2}$ $s=\lceil\sqrt{\Delta t\,\rho/\beta}\rceil$,
$\ell_{s}=\beta s^{2}$

9: $s=\lceil\sqrt{\Delta t\,\rho/\beta}\rceil$, $\ell_{s}=\beta s^{2}$ return
$s$, $\ell_{s}$ return $s$, $\ell_{s}$ return $s$, $\ell_{s}$

10: return $s$, $\ell_{s}$ return $s$, $\ell_{s}$ function RKC_Iteration($t$,
$y$, $\Delta t$, $f$, $s$, $\varepsilon$) return $s$, $\ell_{s}$ return $s$,
$\ell_{s}$ function RKC_Iteration($t$, $y$, $\Delta t$, $f$, $s$,
$\varepsilon$)

11: return $s$, $\ell_{s}$ return $s$, $\ell_{s}$ function RKC_Iteration($t$,
$y$, $\Delta t$, $f$, $s$, $\varepsilon$) $\mu_{j},\nu_{j},\kappa_{j},c_{j}=\
$get_coeffs($s$, $\varepsilon$) $\mu_{j},\nu_{j},\kappa_{j},c_{j}=\
$get_coeffs($s$, $\varepsilon$)

12: $\mu_{j},\nu_{j},\kappa_{j},c_{j}=\ $get_coeffs($s$, $\varepsilon$)
$g_{0}=y$ $g_{0}=y$

13: $g_{0}=y$ $g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$
$g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$

14: $g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$ $g_{1}=g_{0}+\mu_{1}\Delta
tf(t,g_{0})$ $g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$ for $j=2,\ldots,s$ do
$g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$ $g_{1}=g_{0}+\mu_{1}\Delta
tf(t,g_{0})$ for $j=2,\ldots,s$ do

15: $g_{1}=g_{0}+\mu_{1}\Delta tf(t,g_{0})$ $g_{1}=g_{0}+\mu_{1}\Delta
tf(t,g_{0})$ for $j=2,\ldots,s$ do
$g_{j}=\nu_{j}g_{j-1}+\kappa_{j}g_{j-2}+\mu_{j}\Delta tf(t+c_{j-1}\Delta
t,g_{j-1})$ $g_{j}=\nu_{j}g_{j-1}+\kappa_{j}g_{j-2}+\mu_{j}\Delta
tf(t+c_{j-1}\Delta t,g_{j-1})$
$g_{j}=\nu_{j}g_{j-1}+\kappa_{j}g_{j-2}+\mu_{j}\Delta tf(t+c_{j-1}\Delta
t,g_{j-1})$

16: $g_{j}=\nu_{j}g_{j-1}+\kappa_{j}g_{j-2}+\mu_{j}\Delta tf(t+c_{j-1}\Delta
t,g_{j-1})$ return $g_{s}$ return $g_{s}$ return $g_{s}$

17: return $g_{s}$ return $g_{s}$ function get_coeffs($s$, $\varepsilon$)
return $g_{s}$ return $g_{s}$ function get_coeffs($s$, $\varepsilon$)

18: return $g_{s}$ return $g_{s}$ function get_coeffs($s$, $\varepsilon$)
$\omega_{0}=1+\varepsilon/s^{2}$,
$\omega_{1}=T_{s}(\omega_{0})/T_{s}^{\prime}(\omega_{0})$
$\omega_{0}=1+\varepsilon/s^{2}$,
$\omega_{1}=T_{s}(\omega_{0})/T_{s}^{\prime}(\omega_{0})$

19: $\omega_{0}=1+\varepsilon/s^{2}$,
$\omega_{1}=T_{s}(\omega_{0})/T_{s}^{\prime}(\omega_{0})$
$b_{j}=1/T_{j}(\omega_{0})$, $j=0,\ldots,s$ $b_{j}=1/T_{j}(\omega_{0})$,
$j=0,\ldots,s$

20: $b_{j}=1/T_{j}(\omega_{0})$, $j=0,\ldots,s$
$\mu_{1}=\omega_{1}/\omega_{0}$ $\mu_{1}=\omega_{1}/\omega_{0}$

21: $\mu_{1}=\omega_{1}/\omega_{0}$ $c_{0}=0$, $c_{1}=\mu_{1}$ $c_{0}=0$,
$c_{1}=\mu_{1}$

22: $c_{0}=0$, $c_{1}=\mu_{1}$ $c_{0}=0$, $c_{1}=\mu_{1}$ $c_{0}=0$,
$c_{1}=\mu_{1}$ for $j=2,\ldots,s$ do $c_{0}=0$, $c_{1}=\mu_{1}$ $c_{0}=0$,
$c_{1}=\mu_{1}$ for $j=2,\ldots,s$ do

23: $c_{0}=0$, $c_{1}=\mu_{1}$ $c_{0}=0$, $c_{1}=\mu_{1}$ for $j=2,\ldots,s$
do $\mu_{j}=2\omega_{1}b_{j}/b_{j-1}$, $\nu_{j}=2\omega_{0}b_{j}/b_{j-1}$,
$\kappa_{j}=-b_{j}/b_{j-2}$ $\mu_{j}=2\omega_{1}b_{j}/b_{j-1}$,
$\nu_{j}=2\omega_{0}b_{j}/b_{j-1}$, $\kappa_{j}=-b_{j}/b_{j-2}$

24: $\mu_{j}=2\omega_{1}b_{j}/b_{j-1}$, $\nu_{j}=2\omega_{0}b_{j}/b_{j-1}$,
$\kappa_{j}=-b_{j}/b_{j-2}$ $c_{j}=\nu_{j}c_{j-1}+\kappa_{j}c_{j-2}+\mu_{j}$
$c_{j}=\nu_{j}c_{j-1}+\kappa_{j}c_{j-2}+\mu_{j}$
$c_{j}=\nu_{j}c_{j-1}+\kappa_{j}c_{j-2}+\mu_{j}$

25: $c_{j}=\nu_{j}c_{j-1}+\kappa_{j}c_{j-2}+\mu_{j}$ return
$\mu_{j},\nu_{j},\kappa_{j},c_{j}$ return $\mu_{j},\nu_{j},\kappa_{j},c_{j}$
return $\mu_{j},\nu_{j},\kappa_{j},c_{j}$

One step of the RKC method

At each time step, the function RKC_Step of Section 3.1 computes the
approximate solution $y_{n+1}$ at the new time $t_{n+1}$. Its input arguments
are the the solution $y_{n}$ at the current time $t_{n}$, the step size
$\Delta t$, the right-hand side $f$, and an approximation $\rho$ to the
spectral radius of $\partial f/\partial y(t_{n},y_{n})$. The spectral radius
$\rho$ is either known or estimated using any simple method, such as
Gershgorin’s theorem [31, pg. 89] or a nonlinear power iteration, as described
in Appendix A in Appendix A (see also [37, 38, 58, 57, 64]).

In Section 3.1, we note that the main iteration in 16 can be performed using
only three distinct vectors for $g_{j-2}$, $g_{j-1}$, and $g_{j}$. Moreover,
at 9 the number of stages is proportional to the square root of $\rho$,
instead of $\rho$, as in standard explicit (RK) methods. The damping parameter
$\varepsilon>0$ guarantees that the stability polynomial $R_{s}(z)$ satisfies
$|R_{s}(z)|\leq 1-\varepsilon$ along the negative real axis for $z\in[-\beta
s^{2},-\delta]$ and small $\delta>0$. It also slightly extends the stability
domain into the imaginary direction thereby ensuring a stable narrow strip
along the negative real axis; typically, we set $\varepsilon=0.05$. The
method’s coefficients $\mu,\nu,\kappa,c$ depend on Chebyshev polynomials of
the first kind, $T_{j}(x)$. Indeed, the stability polynomial $R_{s}(z)$ is a
shifted and scaled Chebyshev polynomial.

###  3.2. Modified equation and averaged force

Explicit stabilized methods, such as RKC, certainly alleviate the stringent
stability constraint on the time-step due to the stiff term $f_{F}$ in 8.
Nonetheless, the number $s$ of evaluations of the expensive term $f_{S}$ still
depends on the stiffness of $f_{F}$, as $\rho$, and hence $s$, are still
dictated by the stiffness of $f_{F}$. To ensure that the number of evaluations
of $f_{S}$ only depends on its own stiffness, and no longer on that of
$f_{F}$, the evaluation of $f_{F}$ and $f_{S}$ must be decoupled.

Following [3], we therefore introduce the modified equation

$y_{\eta}^{\prime}=f_{\eta}(y_{\eta}),\qquad y_{\eta}(0)=y_{0},$ (10)

with a free parameter $\eta>0$. Here, the averaged force
$f_{\eta}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is defined as

$f_{\eta}(y)=\frac{1}{\eta}(u(\eta)-y),$ (11)

where the auxiliary solution $u(s)$ is defined by the auxiliary equation

$u^{\prime}=f_{F}(u)+f_{S}(y)\quad s\in[0,\eta],\qquad u(0)=y.$ (12)

We now briefly summarize the main properties of the above modified equation
and averaged force; for further details, we refer to [3].

##### Accuracy.

The modified equation 10 is a first-order approximation to 8 with error
$\mathcal{O}(\eta)$, since it follows from 11 and 12 that

$f_{\eta}(y)=f_{S}(y)+\frac{1}{\eta}\int_{0}^{\eta}f_{F}(u(s))\operatorname{d\\!}s=f_{S}(y)+f_{F}(y)+\mathcal{O}(\eta).$
(13)

In fact, 13 also implies that $f_{\eta}$ is a particular average of
$f_{F}+f_{S}$ along the auxiliary solution $u$.

##### Stability.

If the multirate problem 8 is contractive [32, Chapter IV.12], then, under
some assumptions, 10 is also contractive and

$\|y(t)-y_{\eta}(t)\|\leq
C(\eta)\int_{0}^{t}e^{\mu(t-s)}\|f_{F}(y(s))+f_{S}(y(s))\|\operatorname{d\\!}s$
(14)

holds, where $y(t)$ is the exact solution to 8, $\mu<0$, and $C(\eta)\leq 1$
with $C(\eta)=\mathcal{O}(\eta)$ for $\eta\to 0$. Hence, the solutions of 8
and 10 remain close even during long-time integration.

##### Stiffness.

The spectral radius $\rho_{\eta}$ of the Jacobian of $f_{\eta}$ decreases with
increasing $\eta$. Hence, for $\eta$ sufficiently large, it holds

$\rho_{\eta}\leq\rho_{S},$ (15)

where $\rho_{S}$ if the spectral radius of the Jacobian of $f_{S}$. As the
stiffness of 10 then depends only on the slow term $f_{S}$, the numerical
solution of 10 with any explicit method will require fewer right-hand side
evaluations than that of 8.

Clearly, there is a trade-off between choosing $\eta$ sufficiently small to
preserve accuracy in 13 and 14 yet sufficiently large to satisfy the stiffness
condition 15. Fortunately, both conditions are typically satisfied
simultaneously, as 15 is already satisfied for quite small $\eta\leq\Delta t$
(often $\eta\ll\Delta t$), where $\Delta t$ is the step size used to solve the
modified equation 10. Then, 10 not only is a good approximation to 8 but also
much less stiff thanks to 15. Discretization of the modified equation 10
finally leads to the multirate RKC (mRKC) method, first proposed in [3].

###  3.3. Multirate explicit stabilized method algorithm

Instead of solving 8 directly, the multirate explicit stabilized (mRKC) method
[3] solves the modified problem 10 using an RKC method. In doing so, every
evaluation of $f_{\eta}$, $u(\eta)$ is approximated by solving 12 with yet
another RKC method. Special care is taken when choosing the number of stages
$s$ and $m$ needed for 10 and 12, respectively, and the value of $\eta$. For
an accuracy and stability analysis of the method, we refer to [3, 49].

In Section 3.3, we list the mRKC method for the numerical solution of 8. The
input parameters of mRKC_Step are $t_{n}$, $y_{n}$ and $\Delta t$ (as for
RKC_Step), the two right-hand sides $f_{F}$, $f_{S}$ and the spectral radii
$\rho_{F}$, $\rho_{S}$ of their respective Jacobians.

1: function mRKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f_{F}$, $f_{S}$,
$\rho_{F}$, $\rho_{S}$) function mRKC_Step($t_{n}$, $y_{n}$, $\Delta t$,
$f_{F}$, $f_{S}$, $\rho_{F}$, $\rho_{S}$)

2: function mRKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f_{F}$, $f_{S}$,
$\rho_{F}$, $\rho_{S}$) $s,\ \ell_{s}=\ $get_stages($\Delta t$, $\rho_{S}$,
$0.05$) $s,\ \ell_{s}=\ $get_stages($\Delta t$, $\rho_{S}$, $0.05$)

3: $s,\ \ell_{s}=\ $get_stages($\Delta t$, $\rho_{S}$, $0.05$) $\eta=2\Delta
t/\ell_{s}$ $\triangleright$ For RKC, $\ell_{s}=\beta s^{2}\approx 2s^{2}$
$\eta=2\Delta t/\ell_{s}$ $\triangleright$ For RKC, $\ell_{s}=\beta
s^{2}\approx 2s^{2}$

4: $\eta=2\Delta t/\ell_{s}$ $\triangleright$ For RKC, $\ell_{s}=\beta
s^{2}\approx 2s^{2}$ $m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$, $0.05$)
$m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$, $0.05$)

5: $m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$, $0.05$)
$\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$, $\eta$, $f_{F}$,
$f_{S}$, $m$)  $\triangleright$ $\overline{f}_{\eta}$ is the approximated
averaged force 11 $\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$,
$\eta$, $f_{F}$, $f_{S}$, $m$)  $\triangleright$ $\overline{f}_{\eta}$ is the
approximated averaged force 11

6: $\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$, $\eta$, $f_{F}$,
$f_{S}$, $m$)  $\triangleright$ $\overline{f}_{\eta}$ is the approximated
averaged force 11 $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$,
$\overline{f}_{\eta}$, $s$, $0.05$)  $\triangleright$ This call approximates
10 $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$,
$\overline{f}_{\eta}$, $s$, $0.05$)  $\triangleright$ This call approximates
10

7: $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$,
$\overline{f}_{\eta}$, $s$, $0.05$)  $\triangleright$ This call approximates
10 return $y_{n+1}$ return $y_{n+1}$ return $y_{n+1}$

8: return $y_{n+1}$ return $y_{n+1}$ function Averaged_Force($t$, $y$, $\eta$,
$f_{F}$, $f_{S}$, $m$) return $y_{n+1}$ return $y_{n+1}$ function
Averaged_Force($t$, $y$, $\eta$, $f_{F}$, $f_{S}$, $m$)

9: return $y_{n+1}$ return $y_{n+1}$ function Averaged_Force($t$, $y$, $\eta$,
$f_{F}$, $f_{S}$, $m$) $f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y)$ $\triangleright$
$f_{u}$ is the right-hand side of 12, $f_{S}$ is frozen
$f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y)$ $\triangleright$ $f_{u}$ is the right-hand
side of 12, $f_{S}$ is frozen

10: $f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y)$ $\triangleright$ $f_{u}$ is the right-
hand side of 12, $f_{S}$ is frozen $u_{\eta}=\ $RKC_Iteration($t$, $y$,
$\eta$, $f_{u}$, $m$, $0.05$)  $\triangleright$ This call approximates 12
$u_{\eta}=\ $RKC_Iteration($t$, $y$, $\eta$, $f_{u}$, $m$, $0.05$)
$\triangleright$ This call approximates 12

11: $u_{\eta}=\ $RKC_Iteration($t$, $y$, $\eta$, $f_{u}$, $m$, $0.05$)
$\triangleright$ This call approximates 12 return
$\tfrac{1}{\eta}(u_{\eta}-y)$ return $\tfrac{1}{\eta}(u_{\eta}-y)$ return
$\tfrac{1}{\eta}(u_{\eta}-y)$

One step of the mRKC method A single step of Section 3.3 mainly consists of an
outer RKC_Iteration (7), which solves 10 with the approximate averaged force
$\overline{f}_{\eta}$, and an inner RKC_Iteration (11), which solves 12
whenever $\overline{f}_{\eta}$ is needed.

Since the mRKC method is composed of two embedded RKC iterations, the
stabilization process is “multiplicative”, in the sense that the number of
outer RKC iterations or stages, $s$, suffices to stabilize the $f_{S}$ term,
whereas the product of outer and inner RKC iterations, $s\times m$, suffices
to stabilize the $f_{F}$ term. This property is reflected in the computational
cost.

##### Computational cost.

For each evaluation of $\overline{f}_{\eta}$, there is a call to the inner
loop, which involves one evaluation of $f_{S}$ and $m$ evaluations of $f_{F}$.
Since there are $s$ evaluations of $\overline{f}_{\eta}$ during each time
step, there are $s$ evaluations of $f_{S}$ and $s\times m$ evaluations of
$f_{F}$. According to Section 3.3, we compute

$\displaystyle\\#f_{S}\text{ eval.}$ $\displaystyle=s\approx\sqrt{\frac{\Delta
t\rho_{S}}{\beta}},$ (16) $\displaystyle\\#f_{F}\text{ eval.}$
$\displaystyle=s\times m\approx\sqrt{\frac{\Delta
t\rho_{S}}{\beta}}\sqrt{\frac{\eta\rho_{F}}{\beta}}=\sqrt{\frac{\Delta
t\rho_{S}}{\beta}}\sqrt{\frac{2\Delta
t\rho_{F}}{\beta^{2}s^{2}}}\approx\sqrt{\frac{2\Delta
t\rho_{F}}{\beta^{2}}}\approx\sqrt{\frac{\Delta t\rho_{F}}{\beta}}.$ (17)

Hence, the number of $f_{S}$ and $f_{F}$ evaluations depends only on their own
inherent stiffness.

The total number of evaluations of $f_{F}$ is approximately the same as with a
standard RKC method, for which the number of stages is given by $\sqrt{\Delta
t\rho/\beta}$ with $\rho\approx\rho_{F}$. In contrast, the number of
evaluations of $f_{S}$ are drastically reduced. Finally, the “multiplicative”
property of the mRKC method is clearly apparent in 17, as $s\times m$
satisfies precisely the stability condition $\Delta
t\rho_{F}\approx\beta(sm)^{2}$ for $f_{F}$.

###  3.4. The mRKC method for the monodomain model

To apply mRKC to the semi-discrete monodomain model 6 from Section 2.1, we
must first rewrite it as in 8. The monodomain model 6, however, fits the
general form 1, which involves the additional $f_{E}$ term from the ionic
model. In this section, we briefly describe a straightforward but somewhat
simplistic and at times inefficient approach to deal with $f_{E}$ and thus
apply mRKC directly to 6. Then, in Section 4, we shall adapt mRKC to
differential equations as in 1 by judiciously exploiting the particular
structure of $f_{E}$, which will yield a more efficient solver for 6.

To directly apply mRKC to 6, we first need to absorb the $f_{E}$ term from the
ionic model into $f_{F}$ and $f_{S}$. With $y,f_{F},f_{S},f_{E}$ as in 7a, we
define

$\tilde{y}^{\prime}=\tilde{f}_{F}(t,\tilde{y})+\tilde{f}_{S}(t,\tilde{y}),$
(18a) with
$\tilde{f}_{F}(t,y)=f_{F}(t,y)+f_{E}^{1}(t,y),\qquad\tilde{f}_{S}(t,y)=f_{S}(t,y)+f_{E}^{2}(t,y).$
(18b)

In doing so, we have split the ionic model as
$f_{E}^{1}(t,y)+f_{E}^{2}(t,y)=f_{E}(t,y)$, where $f_{E}^{1}$ contains only
the stiffest components while $f_{E}^{2}$ contains all the remaining ones.
Now, mRKC can be applied to the monodomain model 6 by solving 18.

Depending on the ionic model and the mesh size, the term $\tilde{f}_{F}$ in 18
is usually stiffer than $f_{F}$, but still cheap to evaluate since $f_{E}^{1}$
involves only a few components. In contrast, as $f_{E}^{2}$ is significantly
less stiff than $f_{E}^{1}$, the term $\tilde{f}_{S}$ is only mildly stiff but
more expensive to evaluate. Nevertheless, this approach not only remains
unsatisfactory, as it does not exploit the special structure of $f_{E}$, but
also cumbersome to implement because the stiffest gating variables in $f_{E}$
must be identified. In the next section, we shall introduce a different more
efficient mRKC approach, which exploits the particular structure of $f_{E}$
without the need for splitting the gating variables.

For the cardiac models considered here, it turns out that the stiffest
components $f_{E}^{1}$ from the ionic model always reduced to a single
variable for the entire simulation. To identify that stiffest gating variable,
we temporarily set $\tilde{f}_{F}(t,y)=f_{E}^{1}(t,y)$ thereby ignoring
diffusion. Next, selecting one gating variable at a time in $f_{E}^{1}$, we
run the simulation while evaluating the stiffness of $\tilde{f}_{F}$ and
$\tilde{f}_{S}$ at each time step with the nonlinear power iteration – see
Appendix A in the Appendix. The stiffest gating variable is that with the
highest spectral radius for $\tilde{f}_{F}$ over the entire simulation.

##  4\. Exponential Multirate Explicit Stabilized Methods

As shown in Section 3.4, the mRKC method can be directly applied to 1 by
splitting $f_{E}$ and absorbing the two terms into $f_{F}$ and $f_{S}$. Still,
whenever $f_{E}$ is stiffer than $f_{F}$, it will induce a higher number of
stages $m$ in Section 3.3. To improve on the efficiency, we shall now exploit
the special structure 2 of $f_{E}$, as in the Rush–Larsen method [53], by
using an exponential Euler step to integrate $f_{E}$. The combined mRKC /
exponential Euler method then yields the exponential mRKC (emRKC) method for
1, and hence for the monodomain equation 6.

###  4.1. Exponential Euler method

A very efficient method for the numerical solution of

$y^{\prime}=f_{E}(t,y),\qquad y(0)=y_{0},$ (19)

with $f_{E}(y)=\Lambda(y)(y-y_{\infty}(y))$ as in 2 and $\Lambda(y)$ a
diagonal matrix, is the exponential Euler method [33]:

$y_{n+1}=y_{n}+\Delta t\,\varphi(\Delta t\,\Lambda(y_{n}))f_{E}(y_{n}),$ (20)

where

$\varphi(z)=\frac{e^{z}-1}{z}$ (21)

is an exponential-like function. In contrast to the explicit Euler method,
$y_{n+1}=y_{n}+\Delta t\,f_{E}(y_{n})$, the exponential Euler method employs
the matrix exponential $\varphi(\Delta t\,\Lambda(y_{n}))$, which dampens all
large eigenvalues (in magnitude) and stabilizes the discrete dynamical system.
Hence, the exponential Euler method is not only unconditionally stable but
here also relatively cheap because $\Lambda(y)$ is diagonal. Higher-order
variants of this method exist, we refer to [33] for a review article.

###  4.2. Exponential Euler step for the auxiliary equation

Instead of splitting $f_{E}$, as in Section 3.4, we could also entirely
include it into the “fast” term $f_{F}$. Hence, the auxiliary equation 12
would read

$u^{\prime}=f_{F}(u)+f_{S}(y)+f_{E}(u)\quad t\in[0,\eta],\qquad u(0)=y.$ (22)

As in Section 3.4, however, the cost would again increase whenever $f_{E}$ is
stiffer than $f_{F}$.

To design an explicit multirate method for 1, which avoids any further
stability restriction due to $f_{E}$ while taking advantage of its particular
diagonal structure, we now modify the auxiliary problem 22 by taking
inspiration from splitting methods as follows: first, advance the $f_{E}$ term
with the exponential Euler method and then integrate the remaining terms. This
yields the split auxiliary equation:

$\displaystyle y_{E}$ $\displaystyle=y+\eta\,\varphi(\eta\Lambda(y))f_{E}(y),$
(23a) $\displaystyle u^{\prime}$ $\displaystyle=f_{F}(u)+f_{S}(y_{E})\quad
t\in[0,\eta],\qquad u(0)=y_{E}.$ (23b)

In contrast to 22, the stiffness of 23b is independent from the $f_{E}$ term;
hence, approximating $u$ with any explicit method is cheaper.

The emRKC method is derived similarly to mRKC though by considering the split
auxiliary equation 23 instead of 22. Hence, emRKC also solves the modified
equation 10 instead of the original problem 1 and the definition 11 of the
averaged force $f_{\eta}$ remains identical, except that $u(\eta)$ in 11 is
now computed from 23.

Alternatively, we may also split the auxiliary equation with $y_{E}$ as in 23a
yet with 23b replaced by
$u^{\prime}=f_{F}(u)+f_{S}(y_{E})+\varphi(\eta\Lambda(y))f_{E}(y)$ with
$u(0)=y$; hence, the contribution of $\varphi(\eta\Lambda(y))f_{E}(y)$ is
added progressively during the integration. In our numerical experiments, both
splittings led to a similar performance.

###  4.3. Pseudo-code for the exponential multirate explicit stabilized method

In Section 4.3, we list a pseudo-code for the emRKC method. Its main
differences with respect to Section 3.3 are: first, the additional input
argument, $f_{E}$, for emRKC_Step and Averaged_Force; second, the redefinition
of the $f_{u}$ term at 10 and 11; third, the different input argument in 12.
To avoid round-off errors due to small diagonal values in $\Lambda(y)$ or
small $\eta$, the computation of the exponential Euler step is done employing
the identity indicated in 14.

1: function emRKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f_{F}$, $f_{S}$,
$f_{E}$, $\rho_{F}$, $\rho_{S}$) function emRKC_Step($t_{n}$, $y_{n}$, $\Delta
t$, $f_{F}$, $f_{S}$, $f_{E}$, $\rho_{F}$, $\rho_{S}$)

2: function emRKC_Step($t_{n}$, $y_{n}$, $\Delta t$, $f_{F}$, $f_{S}$,
$f_{E}$, $\rho_{F}$, $\rho_{S}$) $s,\ \ell_{s}=\ $get_stages($\Delta t$,
$\rho_{S}$, $0.05$) $s,\ \ell_{s}=\ $get_stages($\Delta t$, $\rho_{S}$,
$0.05$)

3: $s,\ \ell_{s}=\ $get_stages($\Delta t$, $\rho_{S}$, $0.05$) $\eta=2\Delta
t/\ell_{s}$ $\eta=2\Delta t/\ell_{s}$

4: $\eta=2\Delta t/\ell_{s}$ $m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$,
$0.05$) $m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$, $0.05$)

5: $m,\ell_{m}=\ $get_stages($\eta$, $\rho_{F}$, $0.05$)
$\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$, $\eta$, $f_{F}$,
$f_{S}$, $f_{E}$, $m$) $\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$,
$\eta$, $f_{F}$, $f_{S}$, $f_{E}$, $m$)

6: $\overline{f}_{\eta}(t,y)=\ $Averaged_Force($t$, $y$, $\eta$, $f_{F}$,
$f_{S}$, $f_{E}$, $m$) $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$,
$\overline{f}_{\eta}$, $s$, $0.05$) $y_{n+1}=\ $RKC_Iteration($t_{n}$,
$y_{n}$, $\Delta t$, $\overline{f}_{\eta}$, $s$, $0.05$)

7: $y_{n+1}=\ $RKC_Iteration($t_{n}$, $y_{n}$, $\Delta t$,
$\overline{f}_{\eta}$, $s$, $0.05$) return $y_{n+1}$ return $y_{n+1}$ return
$y_{n+1}$

8: return $y_{n+1}$ return $y_{n+1}$ function Averaged_Force($t$, $y$, $\eta$,
$f_{F}$, $f_{S}$, $f_{E}$, $m$) return $y_{n+1}$ return $y_{n+1}$ function
Averaged_Force($t$, $y$, $\eta$, $f_{F}$, $f_{S}$, $f_{E}$, $m$)

9: return $y_{n+1}$ return $y_{n+1}$ function Averaged_Force($t$, $y$, $\eta$,
$f_{F}$, $f_{S}$, $f_{E}$, $m$) $y_{E}=\ $Exponential_Euler($y$, $\eta$,
$\Lambda(y)$, $y_{\infty}(y)$) $y_{E}=\ $Exponential_Euler($y$, $\eta$,
$\Lambda(y)$, $y_{\infty}(y)$)

10: $y_{E}=\ $Exponential_Euler($y$, $\eta$, $\Lambda(y)$, $y_{\infty}(y)$)
$f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y_{E})$ $f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y_{E})$

11: $f_{u}(r,u)=f_{F}(r,u)+f_{S}(t,y_{E})$ $u_{\eta}=\ $RKC_Iteration($t$,
$y_{E}$, $\eta$, $f_{u}$, $m$, $0.05$) $u_{\eta}=\ $RKC_Iteration($t$,
$y_{E}$, $\eta$, $f_{u}$, $m$, $0.05$)

12: $u_{\eta}=\ $RKC_Iteration($t$, $y_{E}$, $\eta$, $f_{u}$, $m$, $0.05$)
return $\tfrac{1}{\eta}(u_{\eta}-y)$ return $\tfrac{1}{\eta}(u_{\eta}-y)$
return $\tfrac{1}{\eta}(u_{\eta}-y)$

13: return $\tfrac{1}{\eta}(u_{\eta}-y)$ return $\tfrac{1}{\eta}(u_{\eta}-y)$
function Exponential_Euler($y$, $\eta$, $\Lambda$, $y_{\infty}$) return
$\tfrac{1}{\eta}(u_{\eta}-y)$ return $\tfrac{1}{\eta}(u_{\eta}-y)$ function
Exponential_Euler($y$, $\eta$, $\Lambda$, $y_{\infty}$)

14: return $\tfrac{1}{\eta}(u_{\eta}-y)$ return $\tfrac{1}{\eta}(u_{\eta}-y)$
function Exponential_Euler($y$, $\eta$, $\Lambda$, $y_{\infty}$) return
$y+(e^{\eta\Lambda}-I)(y-y_{\infty})$ $\triangleright$
$\eta\varphi(\eta\Lambda)\Lambda(y-y_{\infty})=(e^{\eta\Lambda}-I)(y-y_{\infty})$
return $y+(e^{\eta\Lambda}-I)(y-y_{\infty})$ $\triangleright$
$\eta\varphi(\eta\Lambda)\Lambda(y-y_{\infty})=(e^{\eta\Lambda}-I)(y-y_{\infty})$
return $y+(e^{\eta\Lambda}-I)(y-y_{\infty})$ $\triangleright$
$\eta\varphi(\eta\Lambda)\Lambda(y-y_{\infty})=(e^{\eta\Lambda}-I)(y-y_{\infty})$

One step of the emRKC method

###  4.4. Stability analysis

Here we analyse the stability of the emRKC method starting from previous
results for mRKC. In particular, we show that there are no constraints on the
step size $\Delta t$ for sufficiently many internal stages $s,m$. Hence, we
consider the standard linear scalar test equation,

$y^{\prime}=\lambda_{F}y+\lambda_{S}y+\lambda_{E}y,\qquad y(0)=y_{0},$ (24)

where $\lambda_{F},\lambda_{S},\lambda_{E}\leq 0$ correspond to
$f_{F},f_{S},f_{E}$, respectively. Before tackling the analysis for the emRKC
method, we briefly recall the stability analysis for the mRKC method [3, 49].

####  4.4.1. Stability analysis for the mRKC method

The stability analysis of the mRKC scheme [3, 49] also relies on 24 yet with
$\lambda_{E}=0$; hence, we consider the test equation

$y^{\prime}=\lambda_{F}y+\lambda_{S}y,\qquad y(0)=y_{0}.$ (25)

First, we recall a key result from [3], crucial for proving the stability of
the mRKC method when applied to the test equation 25.

######  Lemma 4.1.

Let $P_{m}(z)$ be the stability polynomial of the inner RKC iteration of mRKC
(11 of Section 3.3). Let $\Phi_{m}(z)=(P_{m}(z)-1)/z$. If

$-\ell_{s}\leq\Delta t\,\Phi_{m}(\eta\lambda_{F})(\lambda_{F}+\lambda_{S})\leq
0$ (26)

holds, the mRKC method applied to the test equation 25 is stable.

###### Sketch of proof..

Since mRKC is equivalent to an $s$-stage RKC method applied to the averaged
force $\overline{f}_{\eta}(y)$ (7 of Section 3.3), we conclude that mRKC is
stable if

$-\ell_{s}\leq\Delta t\frac{\partial\overline{f}_{\eta}}{\partial y}(y)\leq
0.$ (27)

To complete the proof, one shows that
$\overline{f}_{\eta}(y)=\Phi_{m}(\eta\lambda_{F})(\lambda_{F}+\lambda_{S})y$
when mRKC is applied to the test equation 25; thus, condition 26 is equivalent
to 27. ∎

In [49], the conditions on $s,m,\eta$ to satisfy 26 were studied in two
different scenarios. In the first case, $\lambda_{F},\lambda_{S}$ can take any
non-positive values and general stability conditions on $s,m,\eta$ are
derived. In the second case, it is assumed that $\lambda_{F}\ll\lambda_{S}\leq
0$, thus yielding milder stability conditions. We recall these conditions in
Lemma 4.2 below.

######  Lemma 4.2.

Let $\rho_{F}=|\lambda_{F}|$ and $\rho_{S}=|\lambda_{S}|$, $\ell_{s}=\beta
s^{2}$ and $\ell_{m}=\beta m^{2}$. Define the hypotheses:

  1. H1

: $\lambda_{F},\lambda_{S}\leq 0$ and $s,m,\eta$ such that

$\Delta
t\,\rho_{S}\leq\ell_{s},\qquad\eta\,\rho_{F}\leq\ell_{m},\qquad\eta=\frac{6\Delta
t}{\ell_{s}}\frac{m^{2}}{m^{2}-1}.$ (28)

  2. H2

: $\lambda_{F}\ll\lambda_{S}\leq 0$ and $s,m,\eta$ such that

$\Delta
t\,\rho_{S}\leq\ell_{s},\qquad\eta\,\rho_{F}\leq\ell_{m},\qquad\eta=\frac{2\Delta
t}{\ell_{s}}.$ (29)

If either 28 or 29 are true, then 26 is satisfied.

For the proof of this lemma and more details on the condition
$\lambda_{F}\ll\lambda_{S}\leq 0$111Essentially we need $\lambda_{F}\leq
4\lambda_{S}$. we refer to [49]. By combining Lemmas 4.1 and 4.2, we
eventually prove the stability of the mRKC method for the test equation.

######  Theorem 4.3.

If either hypothesis 28 or 29 in Lemma 4.2 is satisfied, the mRKC method
applied to the test equation 25 is stable.

###### Proof.

Follows from Lemmas 4.1 and 4.2. ∎

Note that 29 places stronger restrictions on $\lambda_{F},\lambda_{S}$, but
also yields a smaller $\eta$ than 28, and hence a smaller number of stages $m$
which results in fewer $f_{F}$ function evaluations and therefore a more
efficient algorithm. For the context of this work, it makes sense to assume
that 29 is satisfied, hence $\lambda_{F}\ll\lambda_{S}\leq 0$, which stands
for $f_{F}$ significantly stiffer than $f_{S}$. Indeed, it was found in [3],
that in general for systems of ODEs stemming from the spatial discretization
of parabolic PDEs, the stability conditions 29 yield stable solutions. Note
that 29 are the constraints encoded in Sections 3.3 and 4.3.

We also remark that there is no constraint on the step size $\Delta t$.
Indeed, given $\rho_{F},\rho_{S}$, one can always choose $\Delta t$ according
to the desired accuracy. Stability is then achieved by setting $s$
sufficiently large to ensure $\Delta t\rho_{S}\leq\ell_{s}$ and finally
computing $\eta$ and $m$.

####  4.4.2. Stability analysis for the emRKC method

The stability analysis of the emRKC method follows from that for the mRKC
method. In particular, we find that the same stability conditions hold and
therefore there are no constraints on the step size $\Delta t$.

######  Theorem 4.4.

If eithe hypothesis 28 or 29 in Lemma 4.2 is satisfied, the emRKC method
applied to the test equation 24 is stable.

###### Proof.

Following [3], we deduce for emRKC applied to 24 that

$u_{\eta}=(P_{m}(\eta\lambda_{F})+\Phi_{m}(\eta\lambda_{F})\eta\lambda_{S})y_{E}=(P_{m}(\eta\lambda_{F})+\Phi_{m}(\eta\lambda_{F})\eta\lambda_{S})e^{\eta\lambda_{E}}y.$
(30)

The averaged force $\overline{f}_{\eta}(y)$ is therefore given by

$\displaystyle\overline{f}_{\eta}(y)$
$\displaystyle=\frac{1}{\eta}(u_{\eta}-y)=\frac{1}{\eta}\left((P_{m}(\eta\lambda_{F})+\Phi_{m}(\eta\lambda_{F})\eta\lambda_{S})e^{\eta\lambda_{E}}-1\right)y$
(31)
$\displaystyle=\frac{1}{\eta}\left((P_{m}(\eta\lambda_{F})-1+\Phi_{m}(\eta\lambda_{F})\eta\lambda_{S})e^{\eta\lambda_{E}}+e^{\eta\lambda_{E}}-1\right)y$
$\displaystyle=\left(\Phi_{m}(\eta\lambda_{F})(\lambda_{F}+\lambda_{S})e^{\eta\lambda_{E}}+\varphi(\eta\lambda_{E})\lambda_{E}\right)y.$

Since emRKC is an $s$-stage RKC method applied to $\overline{f}_{\eta}(y)$ (7
of Section 4.3), it is stable if

$-\ell_{s}\leq\Delta t\frac{\partial\overline{f}_{\eta}}{\partial y}(y)\leq
0.$ (32)

By using $\varphi(z)\geq 0$, $\lambda_{E}\leq 0$ and 26, we obtain

$\displaystyle 0$ $\displaystyle\geq\Delta
t\frac{\partial\overline{f}_{\eta}}{\partial y}(y)=\Delta
t\left(\Phi_{m}(\eta\lambda_{F})(\lambda_{F}+\lambda_{S})e^{\eta\lambda_{E}}+\varphi(\eta\lambda_{E})\lambda_{E}\right)$
(33) $\displaystyle\geq-\ell_{s}e^{\eta\lambda_{E}}+\frac{\Delta
t}{\eta}(e^{\eta\lambda_{E}}-1)=-\ell_{s}+(e^{\eta\lambda_{E}}-1)\left(\frac{\Delta
t}{\eta}-\ell_{s}\right)\geq-\ell_{s}.$

The last inequality follows from $e^{\eta\lambda_{E}}-1\leq 0$ and the fact
that $\frac{\Delta t}{\eta}-\ell_{s}\leq 0$ holds for both choices of $\eta$
in 28 or 29. ∎

Note that $\lambda_{E}$ does not impact the choice of $s$,$m$, or $\eta$ due
to its exponential integration.

###  4.5. Accuracy analysis

To prove that the emRKC method is first-order accurate, we assume that
$f_{F},f_{S},f_{E}$ are sufficiently smooth. Furthermore, for simplicity of
presentation, we restrict ourselves to the autonomous case when
$f_{F},f_{S},f_{E}$ are time independent.

######  Theorem 4.5.

The emRKC method is first-order accurate.

###### Proof.

We start by computing a Taylor expansion of the exact auxiliary solution
$u(\eta)$ defined by the split auxiliary equation 23. We have $y_{E}=y+\eta
f_{E}(y)+\mathcal{O}(\eta^{2})$ and $u(s)=u(0)+\mathcal{O}(s)$. Since
$u(0)=y_{E}$, we thus obtain

$\displaystyle u(\eta)$
$\displaystyle=u(0)+\int_{0}^{\eta}u^{\prime}(s)\operatorname{d\\!}s=y_{E}+\int_{0}^{\eta}f_{F}(u(s))+f_{S}(y_{E})\operatorname{d\\!}s=y_{E}+\int_{0}^{\eta}f_{F}(y_{E})+f_{S}(y_{E})+\mathcal{O}(s)\operatorname{d\\!}s$
(34) $\displaystyle=y_{E}+\eta f_{F}(y_{E})+\eta
f_{S}(y_{E})+\mathcal{O}(\eta^{2})=y+\eta(f_{F}(y)+f_{S}(y)+f_{E}(y))+\mathcal{O}(\eta^{2}).$

Since the inner RKC method (12 of Section 4.3) used to approximate the
integral in 34 is first-order accurate, it introduces a local error of
$\mathcal{O}(\eta^{2})$; hence, the numerical auxiliary solution $u_{\eta}$
satisfies

$u_{\eta}=u(\eta)+\mathcal{O}(\eta^{2})=y+\eta(f_{F}(y)+f_{S}(y)+f_{E}(y))+\mathcal{O}(\eta^{2}).$
(35)

Therefore, the numerical averaged force $\overline{f}_{\eta}(y)$ computed in
Section 4.3 satisfies

$\overline{f}_{\eta}(y)=\frac{1}{\eta}(u_{\eta}-y)=f_{F}(y)+f_{S}(y)+f_{E}(y)+\mathcal{O}(\eta).$
(36)

Hence it is a first-order approximation in $\eta$ of the exact right-hand
side. At any time step given $y_{n}$, the outer RKC method (7 of Section 4.3)
computes $y_{n+1}$ as

$\displaystyle g_{0}$ $\displaystyle=y_{n},\quad g_{1}=g_{0}+\mu_{1}\Delta
t\overline{f}_{\eta}(g_{0}),$ (37) $\displaystyle g_{j}$
$\displaystyle=\nu_{j}g_{j-1}+\kappa_{j}g_{j-2}+\mu_{j}\Delta
t\overline{f}_{\eta}(g_{j-1}),\quad j=2,\ldots,s,$ $\displaystyle y_{n+1}$
$\displaystyle=g_{s}.$

By setting $d_{j}=g_{j}-y_{n}$ and using $\nu_{j}+\kappa_{j}=1$ we can rewrite
37 as

$\displaystyle d_{0}$ $\displaystyle=0,\quad d_{1}=\mu_{1}\Delta
t\overline{f}_{\eta}(t_{n},y_{n})$ (38) $\displaystyle d_{j}$
$\displaystyle=\nu_{j}d_{j-1}+\kappa_{j}d_{j-2}+\mu_{j}\Delta
t\overline{f}_{\eta}(y_{n}+d_{j-1}),\quad j=2,\ldots,s,$ $\displaystyle
y_{n+1}$ $\displaystyle=y_{n}+d_{s}.$

From [16, Lemma 3.1] with $A=0$ and $r_{j}=\mu_{j}\Delta
t\overline{f}_{\eta}(y_{n}+d_{j-1})$ we obtain

$y_{n+1}=y_{n}+\sum_{j=1}^{s}\frac{b_{s}}{b_{j}}U_{s-j}(\omega_{0})r_{j}=y_{n}+\Delta
t\sum_{j=1}^{s}\frac{b_{s}}{b_{j}}U_{s-j}(\omega_{0})\mu_{j}\overline{f}_{\eta}(y_{n}+d_{j-1}),$
(39)

where the $b_{j}$ are as in Section 3.1 and $U_{j}(x)$ denotes the Chebyshev
polynomial of the second kind of degree $j$. Since $d_{j}=\mathcal{O}(\Delta
t)$ (shown recursively) it follows

$y_{n+1}=y_{n}+\Delta
t\sum_{j=1}^{s}\frac{b_{s}}{b_{j}}U_{s-j}(\omega_{0})\mu_{j}\left(\overline{f}_{\eta}(y_{n})+\mathcal{O}(\Delta
t)\right)=y_{n}+\Delta t\overline{f}_{\eta}(y_{n})+\mathcal{O}(\Delta t^{2}),$
(40)

where the last equality follows from the identity
$\sum_{j=1}^{s}\frac{b_{s}}{b_{j}}U_{s-j}(\omega_{0})\mu_{j}=1$ (see [16]).
Finally, we use 40 in 36 to obtain

$y_{n+1}=y_{n}+\Delta
t(f_{F}(y_{n})+f_{S}(y_{n})+f_{E}(y_{n}))+\mathcal{O}(\Delta t^{2}+\Delta
t\eta).$ (41)

Since $\eta=\mathcal{O}(\Delta t)$, the local error is $\mathcal{O}(\Delta
t^{2})$ and the method is first-order accurate. ∎

##  5\. Numerical experiments

In this section, we apply the mRKC and emRKC methods to the monodomain model 3
and compare them with the popular implicit-explicit Rush–Larsen (IMEX-RL)
method described below. Their performance depends on a number of factors, such
as the mesh size and regularity, the nonlinearity, the stiffness, and the
number of ionic state variables; moreover, the performance of the IMEX-RL
method also depends on the type of solver used during each time step. To take
into account all these factors, we perform numerical experiments either with
structured or unstructured meshes for different mesh sizes, while considering
ionic models of varying size and stiffness, both in two and three space
dimensions thereby including both direct and iterative linear solvers into the
comparison.

More specifically, we first assess in Section 5.1 the effective stability
properties of the emRKC method when compared with the theoretical results of
Section 4.4. Then, in Section 5.2 we perform a two-dimensional experiment on a
structured mesh, where we compare the various methods for different ionic
models and mesh sizes. In Section 5.3, we repeat the experiment yet for a
three-dimensional structured mesh; hence, an iterative solver is used for the
linear systems needed in the IMEX-RL method. Next in Section 5.4, we consider
a more realistic three-dimensional geometry of the left atrium with an
unstructured mesh. Again we repeat the numerical experiment for the same
geometry but in the presence of unhealthy fibrotic tissue when the diffusion
tensor is discontinuous. We conclude in Section 5.6 with a scalability
experiment.

##### The implicit-explicit-Rush–Larsen method

The well-known IMEX-RL method for the monodomain model will be used as a
benchmark to assess the performance of mRKC and emRKC. It is based on an IMEX
scheme for 6a and the Rush–Larsen approach [53] for the ionic model 6b and 6c;
hence, $\bm{z}_{S}$ is integrated explicitly whereas the gating variables
$\bm{z}_{E}$ via an exponential integrator.

The IMEX-RL method proceeds as follows. First, it solves the ODE subsystem 6b
and 6c with the Rush–Larsen method [53]. Hence, at the beginning of each time
step, the method freezes the electric potential $\bm{V}_{m}$ and solves 6b
with the exponential Euler method. Then it solves 6c with the explicit Euler
method. Finally, it solves 6a by treating the linear term implicitly and the
remaining terms explicitly. A pseudo-code for the resulting IMEX-RL method is
provided in Section 5.

Since the stiffest terms $g_{E}(\bm{V}_{m},\bm{z}_{E})$ and
$\mathbf{A}\bm{V}_{m}$ are treated exponentially and implicitly, respectively,
the scheme only has a mild stability constraint. It is also quite efficient
because $\Lambda_{E}(\bm{V}_{m})$ is diagonal and the implicit step is linear
in the unknown.

1: function IMEX-RL_Step($t_{n}$, $\bm{V}_{m,n}$, $\bm{z}_{E,n}$,
$\bm{z}_{S,n}$, $\Delta t$) function IMEX-RL_Step($t_{n}$, $\bm{V}_{m,n}$,
$\bm{z}_{E,n}$, $\bm{z}_{S,n}$, $\Delta t$)

2: function IMEX-RL_Step($t_{n}$, $\bm{V}_{m,n}$, $\bm{z}_{E,n}$,
$\bm{z}_{S,n}$, $\Delta t$) $\bm{z}_{E,n+1}=\bm{z}_{E,n}+(e^{\Delta
t\Lambda_{E}(\bm{V}_{m,n})}-I)(\bm{z}_{E,n}-\bm{z}_{E,\infty}(\bm{V}_{m,n}))$
$\bm{z}_{E,n+1}=\bm{z}_{E,n}+(e^{\Delta
t\Lambda_{E}(\bm{V}_{m,n})}-I)(\bm{z}_{E,n}-\bm{z}_{E,\infty}(\bm{V}_{m,n}))$

3: $\bm{z}_{E,n+1}=\bm{z}_{E,n}+(e^{\Delta
t\Lambda_{E}(\bm{V}_{m,n})}-I)(\bm{z}_{E,n}-\bm{z}_{E,\infty}(\bm{V}_{m,n}))$
$\bm{z}_{S,n+1}=\bm{z}_{S,n}+\Delta
tg_{S}(\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n})$
$\bm{z}_{S,n+1}=\bm{z}_{S,n}+\Delta
tg_{S}(\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n})$

4: $\bm{z}_{S,n+1}=\bm{z}_{S,n}+\Delta
tg_{S}(\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n})$ Solve for $\bm{V}_{m,n+1}$:
$C_{m}\mathbf{M}\frac{\bm{V}_{m,n+1}-\bm{V}_{m,n}}{\Delta
t}=\chi^{-}1\mathbf{A}\bm{V}_{m,n+1}+\mathbf{M}\bm{I}_{stim}(t_{n})-\mathbf{M}I_{ion}(t_{n},\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n+1})$.
Solve for $\bm{V}_{m,n+1}$:
$C_{m}\mathbf{M}\frac{\bm{V}_{m,n+1}-\bm{V}_{m,n}}{\Delta
t}=\chi^{-}1\mathbf{A}\bm{V}_{m,n+1}+\mathbf{M}\bm{I}_{stim}(t_{n})-\mathbf{M}I_{ion}(t_{n},\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n+1})$.

5: Solve for $\bm{V}_{m,n+1}$:
$C_{m}\mathbf{M}\frac{\bm{V}_{m,n+1}-\bm{V}_{m,n}}{\Delta
t}=\chi^{-}1\mathbf{A}\bm{V}_{m,n+1}+\mathbf{M}\bm{I}_{stim}(t_{n})-\mathbf{M}I_{ion}(t_{n},\bm{V}_{m,n},\bm{z}_{E,n+1},\bm{z}_{S,n+1})$.
return $\bm{V}_{m,n+1}$, $\bm{z}_{E,n+1}$, $\bm{z}_{S,n+1}$ return
$\bm{V}_{m,n+1}$, $\bm{z}_{E,n+1}$, $\bm{z}_{S,n+1}$ return $\bm{V}_{m,n+1}$,
$\bm{z}_{E,n+1}$, $\bm{z}_{S,n+1}$

One step of the IMEX-RL method.

##### Computational setup

The monodomain model 3 is implemented in Python employing the FEniCSx [8, 55,
56] library, where, under the hood, all expensive computations are performed
in C/C++. In particular, all linear algebra operations are performed through
the PETSc [9, 17] library. As ionic models should not be implemented by hand,
due to their high complexity, their definitions were obtained from the CellML
[40] model repository; in doing so, we adapted the C code generated by the
Myokit [11] library and created wrappers to call it from Python. Numerical
experiments of Section 5.2 are performed on our workstation, while those of
Sections 5.3, 5.4, 5.5 and 5.6 on the Piz Daint supercomputer at the Swiss
National Supercomputing Centre (CSCS).

The mRKC method always incorporates the splitting of the gating variables
described in 18b—see also Remark 3.4. In contrast, the emRKC method employs
the same splitting as the IMEX-RL method, which is the usual one used in the
cardiac community.

Unless stated otherwise, the initial values for $\bm{V}_{m}$, $\bm{z}_{E}$,
$\bm{z}_{S}$ are provided by the ionic model’s definition and are uniformly
set throughout the computational domain. Here we shall consider the
Hodgkin–Huxley (HH) [34], the Courtemanche–Ramirez–Nattel (CRN) [15] and the
ten Tusscher–Panfilov (TTP) [62] model. Thus our numerical experiments span a
wide range of situations from the small nonstiff HH model, to the larger and
mildly stiff CRN model, and even to the large and very stiff TTP model, see
Table 2 and [60, 61]. In all numerical experiments, we use the parameter
values specified in Table 1. Finally, the vector fields
$\bm{a}_{l}(\bm{x}),\bm{a}_{t}(\bm{x}),\bm{a}_{n}(\bm{x})$ (see 4), which
define the fibers’ orientation in Section 5.4, are generated with lifex-fiber
[7, 47]. The applied stimulus $\bm{I}_{\text{stim}}(t)$ is specified
individually for each experiment.

| HH | CRN | TTP
---|---|---|---
$N$ | 3 | 20 | 18
$\rho_{\max}$ | 40 | 130 | 1000
Table 2. Ionic models’ size and stiffness, with $N$ the number of state
variables and $\rho_{\max}$ the maximal spectral radius of the right-hand
side’s Jacobian in 6.

The linear systems needed in the IMEX-RL method are solved differently in two
or three space dimensions: in 2–D we use the Cholesky factorization, and in
3–D the conjugate gradient method preconditioned using the algebraic multigrid
library hypre [19]. For mRKC and emRKC, the spectral radii $\rho_{F},\rho_{S}$
are computed only once at the beginning of the simulation employing Appendix A
in Appendix A.

During the numerical experiments, we always monitor the relative error in the
potential $\bm{V}$ at the final time with respect to the $L^{2}(\Omega)$-norm:

$\text{rel.
}L^{2}\text{-err.}=\frac{\|\bm{V^{\star}}-\bm{V}_{N}\|_{L^{2}(\Omega)}}{\|\bm{V^{\star}}\|_{L^{2}(\Omega)}},$
(42)

where $\bm{V}_{N}$ is the numerical solution obtained with either one of the
IMEX-RL, mRKC, emRKC methods and $\bm{V^{\star}}$ is a reference solution to 6
computed with a much smaller time-step $\Delta
t=${10}^{-4}\text{\,}\mathrm{ms}$$ using a method similar to IMEX-RL but fully
explicit, i.e., 5 in Section 5 is performed explicitly; we denote this method
EXEX-RL. The final time is chosen such that a propagating wave is located
inside the domain and the solution is still far from equilibrium.

###  5.1. Practical stability limit

In Section 4.4 we proved for a scalar linear test equation that the emRKC
method is stable for all step sizes $\Delta t>0$ if the numbers of stages
$s,m$ are sufficiently large and satisfy the stability condition 29. In this
first numerical experiment, we study the effective stability properties of
emRKC in a more realistic setting.

To assess the practical stability properties of the emRKC method, we solve the
monodomain equation 3 with increasing step size $\Delta t$ and verify that the
norm of the solution remains bounded, ignoring accuracy here. To push the
method to its stability limit, we choose a setting where the stiffness of
$f_{F}$ and $f_{E}$ is severe, i.e., we choose a relatively small mesh size
$\Delta x=$0.025\text{\,}\mathrm{mm}$$ in a one dimensional domain
$\Omega=[0,20]$\mathrm{mm}$$ and employ the TTP ionic model. Although the
stiffness of $f_{S}$ depends on the ionic model, too, it usually does only
mildly so. We integrate the monodomain equation up to the final time
$T=$100\text{\,}\mathrm{ms}$$ with an applied stimulus $\chi
I_{\text{stim}}(t,\bm{x})=$50\text{\,}\mathrm{\SIUnitSymbolMicro
A}\text{/}{\mathrm{mm}}^{3}$$ for $\bm{x}\in[0,1.5]$\mathrm{mm}$$ and
$t\leq$2\text{\,}\mathrm{ms}$$, and else $\chi I_{\text{stim}}(t,\bm{x})=0$.

In Fig. 1(a) we display the relative norm of the numerical solution
$\mathbf{V}_{N}$ with respect to a reference solution $\mathbf{V}^{*}$ as a
function of the step size $\Delta t$. We see that the method remains stable up
to $\Delta t\leq$4.5\text{\,}\mathrm{ms}$$. In contrast, a standard explicit
scheme, such as the EXEX-RL method described at the beginning of Section 5
(i.e. as in Section 5 but performing 5 explicitly), remains stable only up to
$\Delta t\leq$0.003\text{\,}\mathrm{ms}$$. In Fig. 1(b) we show the number of
stages $s,m$ taken by emRKC. We observe that $m>1$ already for small $\Delta
t$, indicating that the diffusion $f_{F}$ has to be stabilized very soon by
the inner iterations of emRKC. The outer stages are needed to stabilize the
non stiff term $f_{S}$, as indicated by $s>1$ only for relatively large
$\Delta t$. The “V” shaped behavior of $m$ is typical and further discussed in
the next Section 5.2.

In Section 4.4, we proved that for the test equation 24 the step size $\Delta
t$ can be taken arbitrarily large, as the emRKC method remains stable as long
as the number of stages $s,m$ are sufficiently large and satisfy the stability
condition 29. In contrast, as shown in the current experiment, due to
nonlinearity the maximal step size is bounded in practice. Nevertheless, this
practical bound is remarkably large, $\Delta t\leq$4.5\text{\,}\mathrm{ms}$$
when compared to the stability bound for a standard explicit method, where
$\Delta t\leq$0.003\text{\,}\mathrm{ms}$$, thus yielding a maximal step size
increase by a factor of 1500! In that sense, we claim that in practice emRKC
has no stability restriction on the step size, since we rarely use $\Delta
t>$0.1\text{\,}\mathrm{ms}$$ because of accuracy requirements.

(a) Relative norm of $\mathbf{V}_{N}$ as $\Delta t$ increases.

(b) Number of stages $s,m$ as $\Delta t$ increases.

Figure 1. Effective stability. Relative norm of the numerical solution
$\mathbf{V}_{N}$ and number of stages for increasing step size $\Delta t$.
Recall that $s$, $m$ are the number of stages used by the outer and inner RKC
iterations, respectively.

###  5.2. Two-dimensional rectangle

First, we consider the two-dimensional version of the benchmark problem
proposed in [42], hence 3 with $\Omega=[0,20]\times[0,7]\,${\mathrm{mm}}^{2}$$
and the final time $T=$25\text{\,}\mathrm{ms}$$. To initiate the propagation
of the action potential, the applied stimulus is defined as $\chi
I_{\text{stim}}(t,\bm{x})=$50\text{\,}\mathrm{\SIUnitSymbolMicro
A}\text{/}{\mathrm{mm}}^{3}$$ if $t\leq$2\text{\,}\mathrm{ms}$$ and
$\|\bm{x}\|_{\infty}\leq$1.5\text{\,}\mathrm{mm}$$, else $\chi
I_{\text{stim}}(t,\bm{x})=0$.

To begin we assess the accuracy of the various time integration methods, in
particular we verify that the first-order of convergence predicted by Theorem
4.5 also holds in practice. More precisely, for each mesh size $\Delta
x=$0.2\text{\,}\mathrm{mm}$$, $0.1\text{\,}\mathrm{mm}$,
$0.05\text{\,}\mathrm{mm}$ and every ionic model HH, CRN, TTP, we successively
reduce the step size as $\Delta t_{i}=2^{-i}$\mathrm{ms}$$ for $i=0,\ldots,9$,
i.e., from $1\text{\,}\mathrm{ms}$ down to approximately
$0.002\text{\,}\mathrm{ms}$. As shown in Fig. 2, all methods achieve the
expected first-order convergence. Moreover, in all but one case, the emRKC is
equally or even more accurate than the standard IMEX-RL method.

Nonetheless, in Figs. 2(a), 2(d) and 2(g) we also notice that some data points
for large $\Delta t$ are missing for the mRKC and emRKC methods, indicating
instability, yet only for the HH model and not for the more realistic CRN and
TTP models. We believe that instabilities in the HH model are due to the
overly large imaginary part of the eigenvalues of the right-hand side’s
Jacobian, compared to the more realistic CRN and TTP models whose eigenvalues
have relatively small imaginary parts [60, 61]. In fact, those instabilities
become even more pronounced as the mesh size decreases, probably because the
stability of mRKC and emRKC for non-real eigenvalues deterioriates as $f_{F}$
becomes increasingly stiffer; this effect, however, is not well understood
yet, since a stability analysis of mRKC and emRKC with complex eigenvalues is
still lacking.

(a) HH model, $\Delta x=0.2$\mathrm{mm}$$.

(b) CRN model, $\Delta x=0.2$\mathrm{mm}$$.

(c) TTP model, $\Delta x=0.2$\mathrm{mm}$$.

(d) HH model, $\Delta x=0.1$\mathrm{mm}$$.

(e) CRN model, $\Delta x=0.1$\mathrm{mm}$$.

(f) TTP model, $\Delta x=0.1$\mathrm{mm}$$.

(g) HH model, $\Delta x=0.05$\mathrm{mm}$$.

(h) CRN model, $\Delta x=0.05$\mathrm{mm}$$.

(i) TTP model, $\Delta x=0.05$\mathrm{mm}$$.

Figure 2. Two dimensional rectangle. Convergence experiments with different
mesh sizes $\Delta x$ and ionic models.

In Fig. 3, we display the number of stages $s,m$ used by the mRKC and emRKC
methods. For the CRN and TTP models, a closer look reveals that both tend to
lose some accuracy when the number of inner stages $m$ increases. As expected,
the number of outer stages $s$ typically remains the same for both methods,
since $s$ depends on $f_{S}$ for emRKC but on $f_{S}+f_{E}^{2}$ for mRKC, yet
with $f_{E}^{2}$ much less stiff than $f_{S}$. The number of inner stages $m$
for emRKC depends only on the diffusion term $f_{F}$, whereas for mRKC it
depends on $f_{F}+f_{E}^{1}$. In Fig. 3, we also observe that the number of
stages $m$ for the two methods often differ, thus indicating that $f_{E}^{1}$
is indeed stiffer than the diffusion $f_{F}$. In those instances, the
efficiency of mRKC deteriorates against emRKC. We also observe for mRKC that
$m$ depends not only on the ionic model but also on the mesh size, while for
emRKC it depends only on the mesh size. Finally, the “V” shape behavior in the
graph of $m$ apparent at larger step sizes stems from the “multiplicative”
stabilization procedure of the methods, as discussed in Section 3. Indeed, for
decreasing $\Delta t$, a constant number of outer stages $s$ also implies a
decreasing number of inner stages $m$. In contrast, for decreasing $s$, the
number of inner stages $m$ can slightly increase to compensate for the weaker
stabilization of the outer loop. Still, $s\times m$ always decreases with
$\Delta t$.

(a) HH model, $\Delta x=0.2$\mathrm{mm}$$.

(b) CRN model, $\Delta x=0.2$\mathrm{mm}$$.

(c) TTP model, $\Delta x=0.2$\mathrm{mm}$$.

(d) HH model, $\Delta x=0.1$\mathrm{mm}$$.

(e) CRN model, $\Delta x=0.1$\mathrm{mm}$$.

(f) TTP model, $\Delta x=0.1$\mathrm{mm}$$.

(g) HH model, $\Delta x=0.05$\mathrm{mm}$$.

(h) CRN model, $\Delta x=0.05$\mathrm{mm}$$.

(i) TTP model, $\Delta x=0.05$\mathrm{mm}$$.

Figure 3. Two dimensional rectangle. Number of stages taken by the mRKC and
emRKC methods, with different mesh sizes $\Delta x$ and ionic models. Recall
that $s$, $m$ are the number of stages taken by the outer and inner RKC
iterations, respectively.

Finally, in Fig. 4 we display work vs. accuracy by monitoring the $L^{2}$-norm
relative errors against the total CPU time. The emRKC method tends to be the
most efficient scheme except for one instance, while the mRKC method also
performs reasonably well. The high performance of emRKC compared to the other
methods stems from a number of factors. First, in comparison to IMEX-RL, the
explicit emRKC method is cheaper and tends to be more accurate because the
$f_{S},f_{E}$ terms take part in the stabilization procedure for the diffusion
term $f_{F}$ (12 of Section 4.3); thus, a strong coupling between the
different model’s components is guaranteed (not so in the splitting strategy
of IMEX-RL). Second, in comparison to mRKC, emRKC is cheaper due to the
smaller number of inner stages $m$ and has higher accuracy because it employs
the more accurate exponential integrator for the gating variables.

(a) HH model, $\Delta x=0.2$\mathrm{mm}$$.

(b) CRN model, $\Delta x=0.2$\mathrm{mm}$$.

(c) TTP model, $\Delta x=0.2$\mathrm{mm}$$.

(d) HH model, $\Delta x=0.1$\mathrm{mm}$$.

(e) CRN model, $\Delta x=0.1$\mathrm{mm}$$.

(f) TTP model, $\Delta x=0.1$\mathrm{mm}$$.

(g) HH model, $\Delta x=0.05$\mathrm{mm}$$.

(h) CRN model, $\Delta x=0.05$\mathrm{mm}$$.

(i) TTP model, $\Delta x=0.05$\mathrm{mm}$$.

Figure 4. Two dimensional rectangle. Efficiency experiments with different
mesh sizes $\Delta x$ and ionic models.

###  5.3. Three-dimensional cuboid

Next, we perform a series of numerical experiments very similar to those in
Section 5.2 but in three dimensions. Hence, we set
$\Omega=[0,20]\times[0,7]\times[0,3]\,${\mathrm{mm}}^{3}$$ while the applied
stimulus and final simulation time remain identical. Again all methods
achieved first-order convergence in time for all three ionic models – those
plots are omitted here for the sake of brevity. We also omit the graphs of the
number of stages taken by mRKC and emRKC, which are identical to the two-
dimensional case.

The work vs. accuracy diagrams, shown in Fig. 5, again confirm the high
efficiency of emRKC. Still, it tends to be slower than in 2-D probably because
the inner iterations are more expensive when the diffusion matrix stems from a
three-dimensional problem. Again we note that the mRKC and emRKC schemes are
stable for the CRN and TTP models, even at the largest step size $\Delta
t=$1\text{\,}\mathrm{ms}$$, but not for the HH model.

(a) HH model, $\Delta x=0.2$\mathrm{mm}$$.

(b) CRN model, $\Delta x=0.2$\mathrm{mm}$$.

(c) TTP model, $\Delta x=0.2$\mathrm{mm}$$.

(d) HH model, $\Delta x=0.1$\mathrm{mm}$$.

(e) CRN model, $\Delta x=0.1$\mathrm{mm}$$.

(f) TTP model, $\Delta x=0.1$\mathrm{mm}$$.

(g) HH model, $\Delta x=0.05$\mathrm{mm}$$.

(h) CRN model, $\Delta x=0.05$\mathrm{mm}$$.

(i) TTP model, $\Delta x=0.05$\mathrm{mm}$$.

Figure 5. Three dimensional cuboid. Efficiency experiments with different mesh
sizes $\Delta x$ and ionic models.

###  5.4. Left atrium geometry

Here we compare the performance of the mRKC and emRKC methods against the
IMEX-RL benchmark but now consider a realistic left atrium geometry. The
computational mesh was obtained by refining twice the realistic mesh used in
[20, 21], which led to an average mesh size of $0.26\text{\,}\mathrm{mm}$,
maximal mesh size of $0.69\text{\,}\mathrm{mm}$, minimal mesh of size
$0.026\text{\,}\mathrm{mm}$, and about 1.7 million degrees of freedom for the
FE-solution. The vector fields $\bm{a}_{l}(\bm{x})$ and $\bm{a}_{t}(\bm{x})$
which define the fiber and sheet orientation in 4, have been generated with
the lifex-fiber library [7, 47]. The initial value, shown in Fig. 6(a), is no
longer constant throughout $\Omega$ but instead generated by solving the
monodomain equation for $t\in[0,800]\,$\mathrm{ms}$$ and applying the same
stimulus-pacing protocol as in [23, 54] in a fixed region of the domain.

(a) Initial value.

(b) Solution at $t=$50\text{\,}\mathrm{ms}$$.

Figure 6. Left atrium geometry. Initial value and solution at
$t=$50\text{\,}\mathrm{ms}$$ for the CRN ionic model.

As in Sections 5.2 and 5.3, we perform convergence and efficiency experiments
for the three ionic models HH, CRN, and TTP. The final time is set to
$T=$50\text{\,}\mathrm{ms}$$, see Fig. 6(b) for an illustration of the
solution at $t=T$. Here, the applied stimulus is identically zero, as the
initial condition already contains a propagating wave. Again, all three
methods achieve first-order convergence in time for all three models, as shown
in Fig. 7. Unlike the previous cases, however, the mRKC method here is
unstable for the TTP model, whereas the emRKC method always remains stable.
The relative error for emRKC is always comparable if not better (HH case) to
the error for IMEX-RL, which is the standard method in the cardiac community.
Note that relative error of 1% to 5% is comparable to the spatial error, which
is of the same order for linear FE and mesh size of $\approx
0.25\text{\,}\mathrm{mm}$, see [45].

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 7. Left atrium geometry. Convergence experiments with different ionic
models.

In Fig. 8, we observe that the number of stages used by the mRKC and emRKC
methods is higher because of the smaller minimal mesh size. Therefore, the
mRKC method evaluates the stiff part $f_{E}^{1}$ of the ionic model many times
18b, which may be the source of internal instabilities. As the internal
stability analysis of mRKC is not available, this issue is not fully
understood yet. Finally, we note that both explicit multirate methods always
take a number of stages significantly larger than one, indicating that any
standard explicit method would be highly unstable for such large step sizes.
In Fig. 9, the results from the efficiency experiment again suggest that emRKC
is slightly faster than the IMEX-RL method, despite the high number of stages
(Fig. 8) due to the mesh induced increased stiffness.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 8. Left atrium geometry. Number of stages used by the mRKC and emRKC
methods with different ionic models. Recall that $s$, $m$ are the number of
stages taken by the outer and inner RKC iterations, respectively.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 9. Left atrium geometry. Efficiency experiments with different ionic
models.

Finally, we provide some code profiling results for the emRKC and IMEX-RL
methods. First, we measure how much time the emRKC method spends in evaluating
$f_{F}$, $f_{S}$, and in performing the exponential Euler step for $f_{E}$
(see Section 4.3). To do so, we apply emRKC with $\Delta
t=2^{-6}\,$\mathrm{ms}$=$0.015\,625\text{\,}\mathrm{ms}$$ and profile the
emRKC_Step function of Section 4.3. We choose this $\Delta t$ since it is the
largest one yielding a relative error of about 1% (Fig. 7). Moreover, in view
of the efficiency results depicted in Fig. 9, we deduce that the outcome of
the profiling is only weakly affected by the size of $\Delta t$. The results
for the three ionic models are displayed in Fig. 10, where for each graph we
show the average over five runs. The four wedges in the pie chart are as
follows: $F$ represents the time spent in evaluating $f_{F}$, $S$ the time
spent evaluating $f_{S}$, $E$ the time for computing the exponential step for
$f_{E}$, and “other” the time required for all remaining operations, such as
communications and vector additions. The size of each wedge is proportional to
the time spent in the corresponding task, with the total time (spent in the
emRKC_Step routine) at the center of the pie chart. At first sight, the time
spent in $F$ seems to depend on the ionic model, even though the number of
matrix-vector multiplications only depends on the mesh. Still, a closer look
at those timings reveals that this variability is due to large fluctuations in
the time spent computing matrix-vector multiplications. In contrast, the
execution timings for $S$ and $E$ are truly ionic model dependent.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 10. Left atrium geometry. Code profiling of the emRKC method. Wedges
$F$, $S$, $E$ indicate time spent evaluating $f_{F}$, $f_{S}$, and computing
the exponential step, respectively. Remaining operations are represented by
the “other” wedge.

Next, we perform the same profiling experiments for the IMEX-RL code while
measuring the time spent in IMEX-RL_Step of Section 5. The corresponding pie
charts are displayed in Fig. 11, where the total execution time is again
distributed between the four wedges $F$, $S$, $E$ and “other”. Here, $F$
represents the time spent in solving the linear system in Section 5 5, $S$ the
time spent in evaluating the slow components of the ionic model in 4 of
Section 5 , and $E$ the time spent in evaluating the exponential for the
gating variables in 3 of Section 5 .

As expected, both emRKC and IMEX-RL spend approximately the same amount of
time in $S$ and $E$, since the required operations are indeed identical. Due
to the larger number of vector additions, the emRKC method spends more time in
the “other” wedge. However, the total time required by IMEX-RL to solve the
linear systems is about 6-7 times larger than the time spent by emRKC in
evaluating $f_{F}$ (compare the two $F$ wedges). This significant reduction in
execution time implies that the stabilization procedure of emRKC for the
diffusion term $f_{F}$ is much cheaper than solving a linear system during
each implicit Euler step. Remarkably, we achieve this reduction in
computational cost without sacrificing accuracy, as illustrated in Fig. 7.
Overall, emRKC achieves a speed-up factor of about 2 over IMEX-RL for the two
most relevant ionic models CRN and TTP.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 11. Left atrium geometry. Code profiling of the IMEX-RL method. Wedges
$F$, $S$, $E$ indicate time spent solving the linear system for $f_{F}$,
evaluating $f_{S}$, and computing the exponential step, respectively.
Remaining operations are represented by the “other” wedge.

###  5.5. Left atrium with fibrosis

Finally, we repeat the numerical experiments from Section 5.4 using the same
computational mesh but in the presence of fibrosis in the cardiac tissue—a
common condition in patients with history of atrial fibrillation.
Mathematically, fibrosis can be modelled by locally reducing the electric
conductivity tensor and altering some ionic currents in the membrane model.
Here, we incorporate fibrosis only in the conductivity tensor by decoupling
cells in the cross-fiber direction [25, 26, 23]. The spatial heterogeneity is
obtained from a spatially-correlated random field $u\mathrel{\mathop{\mathchar
58\relax}}\Omega\rightarrow[0,1]$, where we set $\sigma_{t}=0$ if
$u(\bm{x})>q$, where $q$ is some quantile value; in our simulations, we use
the median. Following [46, Section 2], we generate the random field $u$, shown
in Fig. 12(a), by solving a stochastic differential equation. The
corresponding fibrosis pattern is shown in Fig. 12(b).

(a) Fibrosis random field $u(\bm{x})$.

(b) Fibrosis pattern, in red where $\sigma_{t}=0$.

(c) Initial value.

(d) Solution at $t=$50\text{\,}\mathrm{ms}$$.

Figure 12. Left atrium with fibrosis. The fibrosis pattern (panel B), has been
obtained by clipping at level 0.5 the random field (panel A). All elements
within the red region in panel (B) are considered “fibrotic” and assigned a
transverse conductivity of 0 (no transverse conduction). We also report the
initial value (panel C) and solution at $t=$50\text{\,}\mathrm{ms}$$ (panel D)
for the CRN ionic model.

The initial value is generated exactly as in Section 5.4 and is displayed in
Fig. 12(c), whereas the solution at the final time
$T=$50\text{\,}\mathrm{ms}$$ is shown in Fig. 12(d). Again in Fig. 13, we
verify first-order convergence for all three methods and models and also
observe instabilities in mRKC but not in emRKC. The number of stages taken by
the methods is omitted here since it essentially coincides with that in Fig.
8. In Fig. 14, we display the results from the efficiency experiments where we
again observe that emRKC is slightly faster than IMEX-RL.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 13. Left atrium with fibrosis. Convergence experiments with different
ionic models.

(a) HH model.

(b) CRN model.

(c) TTP model.

Figure 14. Left atrium with fibrosis. Efficiency experiments with different
ionic models.

###  5.6. Parallel performance and scalability

To illustrate the inherent parallelism of the explicit emRKC method, we now
conduct a strong scalability experiment where for a fixed problem size we
progressively increase the number of processors. Again we consider the left
atrium geometry from Section 5.4 and use the CRN ionic model with 21 state
variables (20 for the ionic model and 1 for the potential). Since the FE mesh
has about $1.7\times 10^{6}$ degrees of freedom (dof’s), the resulting total
problem size is about $36\times 10^{6}$ dof’s.

In Fig. 15, we display for both the emRKC and IMEX-RL methods the total
execution CPU time against the number of processors ranging from 16 to 1024;
hence, we start with about $2.2\times 10^{6}$ dof’s per processor and end up
with only about $35^{\prime}000$ dof’s per processor. Although both methods
yield an optimal speed-up up to 256 processors for this problem size, the
parallel efficiency of the implicit IMEX-RL then declines. Although those
preliminary computational results are limited to 1024 processors and a
relatively modest problem size, they already illustrate not only that the
explicit emRKC method is faster than the IMEX-RL methods, as expected from
Fig. 10 and Fig. 11, but also that the explicit emRKC method exhibits better
parallel scalability than the IMEX-RL method.

Figure 15. Strong scalability. CPU time vs. the number of processors.

##  6\. Conclusion

Explicit stabilized methods, such as RKC, strike a judicious balance between
explicit and implicit methods in particular for the time integration of large-
scale nonlinear parabolic PDE’s. The monodomain equation from cardiac
electrophysiology is in principle well-suited for explicit stabilized methods
due to its moderately stiff diffusive part. However, in the presence of
coupled stiff nonlinear ionic models, which are inherently multiscale in
nature, conventional RKC methods can become ineffective. By adapting the
number of stages locally to individual components depending on their
respective degree of stiffness, multirate RKC methods (mRKC) [3] significantly
mitigate this efficiency loss, yet without exploiting the particular diagonal
structure of ionic models as in the standard Rush–Larsen approach.

To further improve upon the mRKC approach, we have proposed the emRKC method
tailored specifically to the monodomain model. Drawing inspiration from the
Rush–Larsen approach, the emRKC method leverages the special diagonal
structure of ionic models for which it applies exponential integration. The
combination of mRKC with exponential time integration results in a
particularly efficient approach for solving the monodomain equations.
Moreover, since it is fully explicit, the emRKC is particularly well-suited
for implementation on GPU architectures in a parallel HPC environment. As
existing fully explicit solvers for the monodomain model are always limited to
structured grids [36], the emRKC thus removes this limitation and extends
explicit time integration to unstructured grids, which are crucial for
electro-mechanical [22, 48] and sub-cellular [51, 52, 35] models.

The stability and accuracy properties of the emRKC method have been analyzed
in Sections 4.4 and 4.5. In particular, in Theorem 4.4 we show that for a
linear scalar test equation the method is stable without constraint on the
step size for sufficiently many stages. We also prove in Theorem 4.5 that the
method is first-order accurate. Both stability and accuracy are also verified
numerically. In particular, we show in Section 5.1 through a realistic
nonlinear example that emRKC can take up to 1500 times larger time steps than
a standard explicit method.

In our numerical experiments, emRKC outperformed or was on par with a standard
baseline method, which combines IMEX integration for the parabolic component
with exponential integration for the ionic model, as such a fair benchmark for
comparison. All numerical experiments unequivocally showcased the good
performance of emRKC, even in the presence of severe stiffness and
nonlinearity.

The results of our study provide compelling evidence that explicit methods can
indeed be exploited effectively for the solution of large-scale, stiff,
nonlinear equations and can be a robust and efficient choice for simulating
multiscale dynamics. Furthermore, the outcomes of this study provide a
foundation for the design and implementation of explicit parallel-in-time
methods specifically tailored for addressing the monodomain equation.

##  Appendix A Nonlinear power iteration

Here we provide a pseudo-code for the nonlinear power iteration [37, 38, 58,
57, 64] used in the RKC, mRKC, and emRKC methods to estimate the spectral
radii of the right-hand sides $f$, $f_{F}$, $f_{S}$. The function Compute_Rho
below takes as input parameters the time $t$ and state $y$, for which the
spectral radius of $\tfrac{\partial f}{\partial y}(t,y)$ needs to be
estimated, the nonlinear function $f$, and an initial guess $v,\rho$ for the
dominant eigenvector of $\tfrac{\partial f}{\partial y}(t,y)$ and its spectral
radius. If no initial guess $v,\rho$ is available, we simply choose an
arbitrary random vector and set $\rho=0$. The computed values $v,\rho$
returned by Compute_Rho provide an initial guess for the next function call.

1: function compute_rho($t$, $y$, $f$, $v$, $\rho$) function compute_rho($t$,
$y$, $f$, $v$, $\rho$)

2: function compute_rho($t$, $y$, $f$, $v$, $\rho$) $n_{y}=\|y\|$,
$n_{v}=\|v\|$ $n_{y}=\|y\|$, $n_{v}=\|v\|$

3: $n_{y}=\|y\|$, $n_{v}=\|v\|$ $\epsilon=10^{-8}$, $\delta=n_{y}\epsilon$,
$q=\delta/n_{v}$ $\epsilon=10^{-8}$, $\delta=n_{y}\epsilon$, $q=\delta/n_{v}$

4: $\epsilon=10^{-8}$, $\delta=n_{y}\epsilon$, $q=\delta/n_{v}$ $tol=10^{-2}$,
$\rho_{\text{old}}=0$ $tol=10^{-2}$, $\rho_{\text{old}}=0$

5: $tol=10^{-2}$, $\rho_{\text{old}}=0$ $tol=10^{-2}$, $\rho_{\text{old}}=0$
$tol=10^{-2}$, $\rho_{\text{old}}=0$ while $|\rho-\rho_{\text{old}}|\geq
tol\cdot\rho$ do $tol=10^{-2}$, $\rho_{\text{old}}=0$ $tol=10^{-2}$,
$\rho_{\text{old}}=0$ while $|\rho-\rho_{\text{old}}|\geq tol\cdot\rho$ do

6: $tol=10^{-2}$, $\rho_{\text{old}}=0$ $tol=10^{-2}$, $\rho_{\text{old}}=0$
while $|\rho-\rho_{\text{old}}|\geq tol\cdot\rho$ do $\rho_{\text{old}}=\rho$,
$v_{\text{old}}=v$ $\rho_{\text{old}}=\rho$, $v_{\text{old}}=v$

7: $\rho_{\text{old}}=\rho$, $v_{\text{old}}=v$ $z=y+q\cdot v_{\text{old}}$
$\triangleright$ $\|z-y\|=\delta$ $z=y+q\cdot v_{\text{old}}$ $\triangleright$
$\|z-y\|=\delta$

8: $z=y+q\cdot v_{\text{old}}$ $\triangleright$ $\|z-y\|=\delta$
$v=f(t,z)-f(t,y)$ $\triangleright$ $v\approx\tfrac{\partial f}{\partial
y}(t,y)(z-y)=q\tfrac{\partial f}{\partial y}(t,y)v_{\text{old}}$
$v=f(t,z)-f(t,y)$ $\triangleright$ $v\approx\tfrac{\partial f}{\partial
y}(t,y)(z-y)=q\tfrac{\partial f}{\partial y}(t,y)v_{\text{old}}$

9: $v=f(t,z)-f(t,y)$ $\triangleright$ $v\approx\tfrac{\partial f}{\partial
y}(t,y)(z-y)=q\tfrac{\partial f}{\partial y}(t,y)v_{\text{old}}$
$n_{v}=\|v\|$, $q=\delta/n_{v}$ $n_{v}=\|v\|$, $q=\delta/n_{v}$

10: $n_{v}=\|v\|$, $q=\delta/n_{v}$ $\rho=n_{v}/\delta$ $\triangleright$
$n_{v}/\delta=\|f(t,z)-f(t,y)\|/\|z-y\|$ $\rho=n_{v}/\delta$ $\triangleright$
$n_{v}/\delta=\|f(t,z)-f(t,y)\|/\|z-y\|$ $\rho=n_{v}/\delta$ $\triangleright$
$n_{v}/\delta=\|f(t,z)-f(t,y)\|/\|z-y\|$

11: $\rho=n_{v}/\delta$ $\triangleright$
$n_{v}/\delta=\|f(t,z)-f(t,y)\|/\|z-y\|$ return $v,\rho$ return $v,\rho$
return $v,\rho$

Nonlinear power iteration

The function Compute_Rho indeed corresponds to a power iteration, as 9 of
Appendix A implies that

$v\approx\tfrac{\partial f}{\partial y}(t,y)(z-y)=q\tfrac{\partial f}{\partial
y}(t,y)v_{\text{old}},$ (43)

while 11

$\rho=\frac{n_{v}}{\delta}=\frac{\|f(t,z)-f(t,y)\|}{\|z-y\|}\geq\frac{\langle
f(t,z)-f(t,y),z-y\rangle}{\|z-y\|^{2}}\approx\frac{\langle\tfrac{\partial
f}{\partial y}(t,y)q\cdot v_{\text{old}},q\cdot
v_{\text{old}}\rangle}{\|q\cdot
v_{\text{old}}\|^{2}}=\frac{\langle\tfrac{\partial f}{\partial
y}(t,y)v_{\text{old}},v_{\text{old}}\rangle}{\|v_{\text{old}}\|^{2}}$ (44)

approximately bounds the Rayleigh quotient. In practice, $\rho$ is
subsequently multiplied by a safety factor $\alpha>1$, typically
$\alpha=1.05$.

The added computational cost from Appendix A during time integration scheme is
marginal. On the one hand, there is no need to update $v,\rho$ at every time
step, every so often, or even just once, is usually sufficient. On the other
hand, subsequent calls typically converge very fast ($1-2$ iterations) thanks
to the previous values of $v,\rho$ now available as initial guesses.

Acknowledgements The first author thanks the developers of FEniCSx and the
Swiss National Supercomputing Centre staff for their support.

## References

  * [1] A. Abdulle.  Fourth order Chebyshev methods with recurrence relation.  SIAM Journal on Scientific Computing, 23(6):2041–2054, 2002.
  * [2] A. Abdulle.  Explicit Stabilized Runge–Kutta Methods.  In Encyclopedia of Applied and Computational Mathematics, pages 460–468. Springer, 2015.
  * [3] A. Abdulle, M. J. Grote, and G. Rosilho de Souza.  Explicit stabilized multirate method for stiff differential equations.  Mathematics of Computation, 91(338):2681–2714, 2022.
  * [4] A. Abdulle and A. A. Medovikov.  Second order Chebyshev methods based on orthogonal polynomials.  Numerische Mathematik, 18:1–18, 2001.
  * [5] A. Abdulle and G. Vilmart.  PIROCK: A swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise.  Journal of Computational Physics, 242:869–888, 2013.
  * [6] P. Africa, M. Salvador, P. Gervasio, L. Dede’, and A. Quarteroni.  A matrix–free high–order solver for the numerical solution of cardiac electrophysiology.  Journal of Computational Physics, 478:111984, 2023.
  * [7] P. C. Africa, R. Piersanti, M. Fedele, L. Dede’, and A. Quarteroni.  Lifex-fiber: An open tool for myofibers generation in cardiac computational models.  BMC Bioinformatics, 24(1):143, Apr. 2023.
  * [8] M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells.  Unified Form Language: A domain-specific language for weak formulations of partial differential equations.  ACM Transactions on Mathematical Software, 40, 2014.
  * [9] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang.  PETSc Web page, 2023.
  * [10] R. H. Clayton, O. Bernus, E. M. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V. Panfilov, F. B. Sachse, G. Seemann, and H. Zhang.  Models of cardiac tissue electrophysiology: Progress, challenges and open questions.  Progress in Biophysics and Molecular Biology, 104(1-3):22–48, 2011\.
  * [11] M. Clerx, P. Collins, E. de Lange, and P. G. Volders.  Myokit: A simple interface to cardiac cellular electrophysiology.  Progress in Biophysics and Molecular Biology, 120(1-3):100–114, Jan. 2016.
  * [12] G. Cohen, P. Joly, J. E. Roberts, and N. Tordjman.  Higher order triangular finite elements with mass lumping for the wave equation.  SIAM J. Numer. Analysis, 38:2047–2078, 2001.
  * [13] P. Colli Franzone, L. F. Pavarino, and S. Scacchi.  Mathematical cardiac electrophysiology, volume 13.  Springer, 2014.
  * [14] Y. Coudière, C. D. Lontsi, and C. Pierre.  Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology.  Electronic Transactions on Numerical Analysis, 52:342–357, 2020\.
  * [15] M. Courtemanche, R. J. Ramirez, and S. Nattel.  Ionic mechanisms underlying human atrial action potential properties: Insights from a mathematical model.  American Journal of Physiology-Heart and Circulatory Physiology, 275(1):H301–H321, 1998.
  * [16] M. Croci and G. Rosilho de Souza.  Mixed-precision explicit stabilized Runge-Kutta methods for single- and multi-scale differential equations.  Journal of Computational Physics, 464, 2022.
  * [17] L. D. Dalcin, R. R. Paz, P. A. Kler, and A. Cosimo.  Parallel distributed computing using Python.  Advances in Water Resources, 34(9):1124–1139, 2011.
  * [18] T. Dumont, M. Duarte, S. Descombes, M. A. Dronne, M. Massot, and V. Louvet.  Simulation of human ischemic stroke in realistic 3D geometry.  Communications in Nonlinear Science and Numerical Simulation, 18(6):1539–1557, 2013.
  * [19] R. D. Falgout and U. M. Yang.  Hypre: A library of high performance preconditioners.  In International Conference on Computational Science, pages 632–641. Springer, 2002.
  * [20] T. E. Fastl, C. Tobon-Gomez, A. Crozier, J. Whitaker, R. Rajani, K. P. McCarthy, D. Sanchez-Quintana, S. Y. Ho, M. D. O’Neill, G. Plank, M. J. Bishop, and S. A. Niederer.  Personalized computational modeling of left atrial geometry and transmural myofiber architecture.  Medical Image Analysis, 47:180–190, July 2018.
  * [21] T. E. Fastl, C. Tobon-Gomez, A. Crozier, J. Whitaker, R. Rajani, K. P. McCarthy, D. Sanchez-Quintana, S. Y. Ho, M. D. O’Neill, G. Plank, M. J. Bishop, and S. A. Niederer.  Personalized computational finite element models for left atrial electromechanics, Sept. 2021.
  * [22] M. Favino, S. Pozzi, S. Pezzuto, F. W. Prinzen, A. Auricchio, and R. Krause.  Impact of mechanical deformation on pseudo-ECG: a simulation study.  EP Europace, 18(suppl_4):iv77–iv84, 2016.
  * [23] L. Gander, S. Pezzuto, A. Gharaviri, R. Krause, P. Perdikaris, and F. Sahli Costabal.  Fast Characterization of Inducible Regions of Atrial Fibrillation Models With Multi-Fidelity Gaussian Process Classification.  Frontiers in Physiology, 13:757159, Mar. 2022.
  * [24] S. Geevers, W. A. Mulder, and J. J. W. van der Vegt.  New higher-order mass-lumped tetrahedral elements for wave propagation modelling.  SIAM J. Sci. Comput., 40:A2830–A2857, 2018.
  * [25] A. Gharaviri, E. Bidar, M. Potse, S. Zeemering, S. Verheule, S. Pezzuto, R. Krause, J. G. Maessen, A. Auricchio, and U. Schotten.  Epicardial Fibrosis Explains Increased Endo–Epicardial Dissociation and Epicardial Breakthroughs in Human Atrial Fibrillation.  Frontiers in Physiology, 11:68, 2020.
  * [26] A. Gharaviri, S. Pezzuto, M. Potse, S. Verheule, G. Conte, R. Krause, U. Schotten, and A. Auricchio.  Left atrial appendage electrical isolation reduces atrial fibrillation recurrences: simulation study.  Circulation: Arrhythmia and Electrophysiology, 14(1):e009230, 2021\.
  * [27] M. J. Grote and J. Diaz.  Energy Conserving Explicit Local Time Stepping for Second-Order Wave Equations.  SIAM Journal on Scientific Computing, 31(3):1985–2014, 2009.
  * [28] M. J. Grote, M. Mehlin, and T. Mitkova.  Runge–Kutta-based explicit local time-stepping methods for wave propagation.  SIAM Journal on Scientific Computing, 37(2):A747–A775, 2015.
  * [29] M. J. Grote and T. Mitkova.  Explicit local time-stepping methods for Maxwell’s equations.  Journal of Computational and Applied Mathematics, 234(12):3283–3302, Oct. 2010.
  * [30] M. J. Grote and T. Mitkova.  High-order explicit local time-stepping methods for damped wave equations.  J. Comput. Appl. Math., 239:270–289, 2013.
  * [31] E. Hairer, S. P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations I, volume 8.  Springer-Verlag, Berlin, 2 edition, 2008.
  * [32] E. Hairer and G. Wanner.  Solving Ordinary Differential Equations II, volume 14.  Springer-Verlag, Berlin, 2002.
  * [33] M. Hochbruck and A. Ostermann.  Exponential integrators.  Acta Numerica, 19:209–286, 2010.
  * [34] A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of Physiology, 117(4):500–544, Aug. 1952.
  * [35] N. M. M. Huynh, F. Chegini, L. F. Pavarino, M. Weiser, and S. Scacchi.  Convergence analysis of BDDC preconditioners for composite DG discretizations of the cardiac cell-by-cell model.  SIAM Journal on Scientific Computing, 45(6):A2836–A2857, 2023.
  * [36] D. Krause, M. Potse, T. Dickopf, R. Krause, A. Auricchio, and F. Prinzen.  Hybrid Parallelization of a Large-Scale Heart Model.  In Facing the Multicore-Challenge II, pages 120–132. Springer, 2012.
  * [37] B. Lindberg.  IMPEX: A program package for solution of systems of stiff differential equations.  Technical report, Dept. of Information Processing, Royal Inst. of Tech., Stockholm, 1972.
  * [38] B. Lindberg.  IMPEX2: A procedure for solution of systems of stiff differential equations.  Technical Report TRITA-NA-7303, Department of Information Processing, Royal Institute of Technology, Stockholm, Sweden, 1973.
  * [39] L. P. Lindner, T. Gerach, T. Jahnke, A. Loewe, D. Weiss, and C. Wieners.  Efficient time splitting schemes for the monodomain equation in cardiac electrophysiology.  International Journal for Numerical Methods in Biomedical Engineering, 39(2), 2023.
  * [40] C. M. Lloyd, J. R. Lawson, P. J. Hunter, and P. F. Nielsen.  The CellML model repository.  Bioinformatics (Oxford, England), 24(18):2122–2123, 2008.
  * [41] C. D. Meyer, D. S. Balsara, and T. D. Aslam.  A stabilized Runge-Kutta-Legendre method for explicit super-time-stepping of parabolic and mixed equations.  Journal of Computational Physics, 257(PA):594–626, 2014.
  * [42] S. A. Niederer, E. Kerfoot, A. P. Benson, M. O. Bernabeu, O. Bernus, C. Bradley, E. M. Cherry, R. Clayton, F. H. Fenton, A. Garny, E. Heidenreich, S. Land, M. Maleckar, P. Pathmanathan, G. Plank, J. F. Rodríguez, I. Roy, F. B. Sachse, G. Seemann, O. Skavhaug, and N. P. Smith.  Verification of cardiac tissue electrophysiology simulators using an N-version benchmark.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369:4331–4351, 2011.
  * [43] D. Ogiermann, D. Balzani, and L. E. Perotti.  An explicit local space-time adaptive framework for monodomain models in cardiac electrophysiology.  Computer Methods in Applied Mechanics and Engineering, 422:116806, 2024.
  * [44] P. Pathmanathan, G. R. Mirams, J. Southern, and J. P. Whiteley.  The significant effect of the choice of ionic current integration method in cardiac electro-physiological simulations.  International Journal for Numerical Methods in Biomedical Engineering, 27(11):1751–1770, Nov. 2011.
  * [45] S. Pezzuto, J. Hake, and J. Sundnes.  Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology.  International Journal for Numerical Methods in Biomedical Engineering, 32(10):e02762, 2016.
  * [46] S. Pezzuto, A. Quaglino, and M. Potse.  On Sampling Spatially-Correlated Random Fields for Complex Geometries.  In Y. Coudière, V. Ozenne, E. Vigmond, and N. Zemzemi, editors, Functional Imaging and Modeling of the Heart, volume 11504, pages 103–111, Cham, 2019. Springer International Publishing.
  * [47] R. Piersanti, P. C. Africa, M. Fedele, C. Vergara, L. Dedè, A. F. Corno, and A. Quarteroni.  Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations.  Computer Methods in Applied Mechanics and Engineering, 373:113468, 2021.
  * [48] A. Quarteroni, L. Dedè, and F. Regazzoni.  Modeling the cardiac electromechanical function: A mathematical journey.  Bulletin of the American Mathematical Society, 59(3):371–403, 2022\.
  * [49] G. Rosilho De Souza.  Numerical Methods for Deterministic and Stochastic Differential Equations with Multiple Scales and High Contrasts.  PhD thesis, EPFL, 2020.
  * [50] G. Rosilho De Souza.  mRKC: A multirate Runge–Kutta–Chebyshev code, 2022.
  * [51] G. Rosilho de Souza, R. Krause, and S. Pezzuto.  Boundary integral formulation of the cell-by-cell model of cardiac electrophysiology.  Engineering Analysis with Boundary Elements, 158:239–251, 2024\.
  * [52] G. Rosilho de Souza, S. Pezzuto, and R. Krause.  Effect of gap junction distribution, size, and shape on the conduction velocity in a cell-by-cell model for electrophysiology.  In O. Bernard, P. Clarysse, N. Duchateau, J. Ohayon, and M. Viallon, editors, Functional Imaging and Modeling of the Heart, volume 13958 of Lecture Notes in Computer Science, pages 117–126, Cham, 2023. Springer.
  * [53] S. Rush and H. Larsen.  A Practical Algorithm for Solving Dynamic Membrane Equations.  IEEE Transactions on Biomedical Engineering, BME-25(4):389–392, 1978.
  * [54] F. Sahli Costabal, T. Banduc, L. Gander, and S. Pezzuto.  The fibrotic kernel signature: Simulation-free prediction of atrial fibrillation.  In O. Bernard, P. Clarysse, N. Duchateau, J. Ohayon, and M. Viallon, editors, Functional Imaging and Modeling of the Heart, volume 13958 of Lecture Notes in Computer Science, pages 87–96, Cham, 2023. Springer.
  * [55] M. W. Scroggs, I. A. Baratta, C. N. Richardson, and G. N. Wells.  Basix: A runtime finite element basis evaluation library.  Journal of Open Source Software, 7(73):3982, 2022.
  * [56] M. W. Scroggs, J. S. Dokken, C. N. Richardson, and G. N. Wells.  Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes.  ACM Transactions on Mathematical Software, 48(2):18:1–18:23, 2022\.
  * [57] L. F. Shampine.  Lipschitz constants and robust ode codes’, Computational Methods in Nonlinear Mechanics.  In Proceedings of the TICOM Second International Conference, JT Oden, Ed., North-Holland Publishing Company, New York, pages 427–449, 1980.
  * [58] L. F. Shampine.  Diagnosing Stiffness for Runge–Kutta Methods.  SIAM Journal on Scientific and Statistical Computing, 12(2):260–272, Mar. 1991.
  * [59] B. P. Sommeijer, L. Shampine, and J. G. Verwer.  RKC: An explicit solver for parabolic PDEs.  Journal of Computational and Applied Mathematics, 88(2):315–326, 1998.
  * [60] R. J. Spiteri and R. C. Dean.  Stiffness analysis of cardiac electrophysiological models.  Annals of Biomedical Engineering, 38(12):3592–3604, 2010.
  * [61] R. J. Spiteri and R. C. Dean.  Erratum: Stiffness analysis of cardiac electrophysiological models.  Annals of Biomedical Engineering, 40(7):1622–1625, 2012.
  * [62] K. H. W. J. ten Tusscher and A. V. Panfilov.  Alternans and spiral breakup in a human ventricular tissue model.  American Journal of Physiology. Heart and Circulatory Physiology, 291(3):H1088–1100, Sept. 2006.
  * [63] P. J. Van der Houwen and B. P. Sommeijer.  On the internal stability of explicit, $m$-stage Runge–Kutta methods for large $m$-values.  Zeitschrift für Angewandte Mathematik und Mechanik, 60(10):479–485, 1980.
  * [64] J. G. Verwer.  An implementation of a class of stabilized explicit methods for the time integration of parabolic equations.  ACM Transactions on Mathematical Software (TOMS), 6(2):188–205, 1980.
  * [65] J. G. Verwer.  Explicit Runge–Kutta methods for parabolic partial differential equations.  Applied Numerical Mathematics, 22(1-3):359–379, 1996.
  * [66] J. G. Verwer, W. Hundsdorfer, and B. P. Sommeijer.  Convergence properties of the Runge–Kutta–Chebyshev method.  Numerische Mathematik, 57(1):157–178, 1990.