11affiliationtext: Center for Computational Medicine in Cardiology, Euler
Institute, Università della Svizzera italiana, via G. Buffi 13, 6900 Lugano,
Switzerland22affiliationtext: Faculty of Mathematics and Informatics, FernUni,
Schinerstrasse 18, 3900, Brig, Switzerland33affiliationtext: Department of
Mathematics, Università di Trento, via Sommarive 14, 38123, Povo (Trento),
Italy

# Boundary Integral Formulation of the Cell-by-Cell Model of Cardiac
Electrophysiology111This work was supported by the European High-Performance
Computing Joint Undertaking EuroHPC under grant agreement No 955495
(MICROCARD) co-funded by the Horizon 2020 programme of the European Union (EU)
and the Swiss State Secretariat for Education, Research and Innovation.

Giacomo Rosilho de Souza<EMAIL_ADDRESS>Rolf Krause
<EMAIL_ADDRESS>Simone Pezzuto<EMAIL_ADDRESS>
###### Abstract

We propose a boundary element method for the accurate solution of the cell-by-
cell bidomain model of electrophysiology. The cell-by-cell model, also called
Extracellular-Membrane-Intracellular (EMI) model, is a system of reaction-
diffusion equations describing the evolution of the electric potential within
each domain: intra- and extra-cellular space and the cellular membrane. The
system is parabolic but degenerate because the time derivative is only in the
membrane domain. In this work, we adopt a boundary-integral formulation for
removing the degeneracy in the system and recast it to a parabolic equation on
the membrane. The formulation is also numerically advantageous since the
number of degrees of freedom is sensibly reduced compared to the original
model. Specifically, we prove that the boundary-element discretization of the
EMI model is equivalent to a system of ordinary differential equations, and we
consider a time discretization based on the multirate explicit stabilized
Runge–Kutta method. We numerically show that our scheme convergences
exponentially in space for the single-cell case. We finally provide several
numerical experiments of biological interest.

Keywords. Cell-by-cell model • EMI model • Boundary Element Method • Cardiac
Electrophysiology • Gap Junctions

##  1 Introduction

The human heart is composed of billions of electrically-active myocytes.
Altogether, myocytes form a syncytium of cells that enables electrical and
mechanical synchronization of the tissue [24]. Cardiac myocytes are excitable
cells that can react and transmit electric currents to communicate and
coordinate their action. Electrical propagation depends on the conductive
properties of the cytoplasm and the extracellular matrix. Cell-to-cell
conduction occurs via gap junctions, permeable channel mostly distributed in
the myocyte longitudinal direction. Myocyte excitability is due to hundreds of
thousands ion channels embedded in the cellular membrane. The overall
propagation of the cardiac action potential emerges from a balance of
diffusion and transmembrane currents.

Mathematically, cardiac electrophysiology models are systems of reaction-
diffusion equations. The reaction term results from transmembrane currents,
which are voltage-dependent and regulated through a gating mechanism. Ion
channel gating is typically modeled via Hodgkin–Huxley formalism, yielding a
possibly large set of ordinary differential equations. The diffusion term
captures the spatial coordination of the cardiac tissue. The state-of-the-art
model is the bidomain system. Here, the intracellular and extracellular spaces
are superimposed and homogenized [30]. Patient-specific organ-scale
simulations routinely employ the bidomain model and its monodomain
approximation. Despite being physiologically accurate, the bidomain model
fails to capture the sub-cellular tissue organization. The cell-by-cell
bidomain model accounts for the Extracellular-Membrane-Intracellular (EMI)
tissue components as separated (yet coupled) entities [37, 44]. The cell-by-
cell model enables a more accurate description of tissue heterogeneities, a
key aspect in heart failure and atrial fibrillation [36].

The cell-by-cell model presents several challenges. First, its unusual
mathematical formulation showing time dynamics at the boundaries, indeed it
presents an ordinary differential equation (ODE) on the transmembrane boundary
and a constraint on the gap junctions. Second, in addition to the natural
stiffness introduced by the Laplacian, the ionic model introduces stiff
nonlinear multiscale dynamics. Third, a full scale heart model would require
billions of cells leading to an incredibly large system of equations. Hence,
advanced tailored methods must be designed to solve cell-by-cell models.

In the literature, cell-by-cell models have already been solved by means of
the finite element or boundary element method. In the finite element community
the problem has been tackled by Stinstra and collaborators [41, 39, 42, 40]
and more recently by Tveito and collaborators [22, 43, 44, 12] as well. In
both cases the cell-by-cell model was employed to study the effects of the
cells microscopic structure on macroscopic values as conductive velocity or
effective tissue conductivity. Also, the cell-by-cell model was employed to
derive the effective parameters for the bidomain model under different
microstructural conditions [39, 37, 18]. Bécue, Potse and Coudière [4, 5, 6]
compared different gap junctions modelizations and studied existence of
solutions in [3]. In the context of the boundary element method the model was
solved only for very simple and structured geometries, for instance in [14,
21, 28] for a longitudinal array of non-touching cells, in [45] for two a two-
cells model, and more recently, from a theoretical point of view, in [19, 20]
it was analyzed for the case of isolated cells.

In this paper we propose a spatial discretization of the cell-by-cell model
based on the boundary element method (BEM) and reduce it to a single system of
ODEs living only on the transmembrane boundary. The great advantage brought by
the BEM is that only boundaries need to be discretized, leading to much
smaller systems of equations compared to more traditional methods as finite
elements or volumes. We stress that the approach presented here is easily
adapted for different gap-junction boundary conditions [44] or unbounded
extracellular domains (the “infinite bath” approximation). Also, any spatial
discretization method for which Dirichlet-to-Neumann maps can be computed
could be employed instead of the BEM. Compared to previous BEM approaches, our
methodology is independent from the cells structure and reduces any problem to
an ODE on the transmembrane boundary.

This paper is organized as follows. In Section 2 we treat the simple case
where only one myocyte cell is present, the purpose of this section is to
introduce the needed tools and our approach in a simplified setting. In
Section 3 instead we discretize in space the full problem with an arbitrary
number of cells, possibly in contact, and reduce it to a system of ODEs.
Finally, in Section 4 we present some numerical results.

##  2 The single-cell problem

The main purpose of this section is to introduce in a simplified setting the
boundary integral formulation and the boundary element method (BEM) employed
to discretize the full problem, done in Section 3, and as well the approach
used to reduce the space discrete problem into a system of ordinary
differential equations (ODE).

###  2.1 Problem formulation

images/text/one_cell

Figure 1: Geometrical setup of problem Eq. 1. The cell and the extra-cellular
space are respectively denoted by $\Omega_{1}$ and $\Omega_{0}$. The cellular
membrane is $\Gamma_{0}$.

Here we consider the EMI model for a single cell, denoted by the bounded
domain $\Omega_{1}\subset\mathbb{R}^{d}$ with $d=2$, embedded in the
extracellular space, denoted by $\Omega_{0}\subset\mathbb{R}^{d}$. See Fig. 1
for a schematic representation of the single-cell problem. Specifically, we
require that the intra- and extra-cellular domain do not overlap, that is
$\Omega_{0}\cap\Omega_{1}=\emptyset$, and that they share a common boundary
$\Gamma_{0}=\bar{\Omega}_{0}\cap\bar{\Omega}_{1}$. The boundary $\Gamma_{0}$
represents the cellular membrane. The membrane model and temporal dynamic of
the system, due to capacitative currents, is confined on $\Gamma_{0}$. We
finally assume that $\Omega_{0}$ is bounded with exterior boundary
$\Sigma=\partial\Omega_{0}\setminus\Gamma_{0}$. Also, we define
$\Gamma_{1}=\partial\Omega_{1}$ (note that for the single-cell problem
$\Gamma_{1}=\Gamma_{0}$). The single-cell problem reads as follows:

$\displaystyle-\sigma_{1}\Delta u_{1}$ $\displaystyle=0,$
$\displaystyle\text{in $\Omega_{1}$},$ (1a) $\displaystyle-\sigma_{0}\Delta
u_{0}$ $\displaystyle=0,$ $\displaystyle\text{in $\Omega_{0}$},$ (1b)
$\displaystyle\sigma_{1}\partial_{\bm{n}_{1}}u_{1}+\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$
$\displaystyle=0,$ $\displaystyle\text{on $\Gamma_{0}$},$ (1c)
$\displaystyle\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$
$\displaystyle=I_{\text{t}}(V_{\text{0}},z)=C_{\text{m}}\partial_{t}V_{\text{0}}+I_{\text{ion}}(V_{\text{0}},z),\qquad$
$\displaystyle\text{on $\Gamma_{0}$},$ (1d) $\displaystyle u_{1}-u_{0}$
$\displaystyle=V_{\text{0}},$ $\displaystyle\text{on $\Gamma_{0}$},$ (1e)
$\displaystyle\partial_{t}z$ $\displaystyle=g(V_{\text{0}},z),$
$\displaystyle\text{on $\Gamma_{0}$},$ (1f)
$\displaystyle\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$ $\displaystyle=0,\qquad$
$\displaystyle\text{on $\Sigma$},$ (1g)

where $u_{0}(x,t)$, $u_{1}(x,t)$, and $V_{\text{0}}(x,t)$ are respectively the
intra-cellular, extra-cellular, and transmembrane electric potential,
$\sigma_{0}$ and $\sigma_{1}$ are respectively the extra- and intra-cellular
electric conductivity, $C_{\text{m}}$ is the membrane capacitance, and
$\bm{n}_{i}$, $i=0,1$ is the outwards normal. The ionic model is defined by
$z(x,t)$, a vector of gating and concentration variables, its corresponding
dynamic in Eq. 1f, and the ionic currents $I_{\text{ion}}(V_{\text{0}},z)$.

The global well-posedness of the problem (1) in Sobolev spaces has been
studied by matano2011 [29]. The idea of the proof is similar to what we do
here, in the sense that the authors recast (1) to an ODE on the interface
$\Gamma_{0}$. The equation involves a pseudo-differential operator not
dissimilar to the Dirichlet-Neumann map, as done below in the discrete
settings with the operator $\mathcal{\psi}$. We are not aware of similar
results for the multi-cell problem.

###  2.2 Boundary integral formulation

Problem in Eq. 1 has already been tackled and carefully analyzed by HJA17
[20], where the BEM with a Galerkin approach was employed. Here, we derive a
boundary integral formulation of the unicellular problem Eq. 1 in terms of
trace operators and Poincaré–Steklov operators.

Let $\gamma^{1}_{t}$ be the trace operator and $\gamma^{1}_{n}$ the conormal
derivative on the boundary $\Gamma_{1}$. More specifically, we introduce the
operators as follows:

$\displaystyle\gamma^{1}_{t}$ $\displaystyle\colon H^{1}(\Omega_{1})\to
H^{1/2}(\Gamma_{1}),\quad$ $\displaystyle\gamma^{1}_{t}u_{1}(\bm{x})$
$\displaystyle=\lim_{\Omega_{1}\ni\bm{y}\to\bm{x}\in\Gamma_{1}}u_{1}(\bm{y}),$
(2) $\displaystyle\gamma^{1}_{n}$ $\displaystyle\colon H^{1}(\Omega_{1})\to
H^{-1/2}(\Gamma_{1}),$ $\displaystyle\gamma^{1}_{n}u_{1}(\bm{x})$
$\displaystyle=\lim_{\Omega_{1}\ni\bm{y}\to\bm{x}\in\Gamma_{1}}\langle\nabla
u_{1}(\bm{y}),\bm{n}_{1}\rangle,$

where the limits on the right must hold for smooth enough $u_{1}$. Let
$G(\bm{x},\bm{y})$ be the fundamental solution of the Laplacian in
$\mathbb{R}^{d}$, the Green representation formula for $u_{1}$ satisfying Eq.
1a implies that

$\displaystyle u_{1}(\bm{x})$
$\displaystyle=\int_{\Gamma_{1}}\gamma^{1}_{t,\bm{y}}G(\bm{x},\bm{y})\gamma^{1}_{n}u_{1}(\bm{y})\operatorname{d\\!}s_{\bm{y}}-\int_{\Gamma_{1}}\gamma^{1}_{n,\bm{y}}G(\bm{x},\bm{y})\gamma^{1}_{t}u_{1}(\bm{y})\operatorname{d\\!}s_{\bm{y}},$
$\displaystyle\bm{x}\in\Omega_{1}$ (3)

where $\bm{y}$ in $\gamma^{1}_{t,\bm{y}},\gamma^{1}_{n,\bm{y}}$ means that the
operators are applied to the second variable of $G(\bm{x},\bm{y})$. Taking the
trace $\gamma^{1}_{t}$ of Eq. 3 we obtain

$\gamma^{1}_{t}u_{1}=\mathcal{V}_{1}\gamma^{1}_{n}u_{1}-(\mathcal{K}_{1}-\tfrac{1}{2}I)\gamma^{1}_{t}u_{1},$
(4)

where $I$ is the identity operator, $\mathcal{V}_{1}$ and $\mathcal{K}_{1}$
are the single and double layer operators defined by

$\displaystyle\mathcal{V}_{1}$ $\displaystyle\colon H^{-1/2}(\Gamma_{1})\to
H^{1/2}(\Gamma_{1}),\quad$ $\displaystyle\mathcal{V}_{1}\rho(\bm{x})$
$\displaystyle=\int_{\Gamma_{1}}\gamma^{1}_{t,\bm{y}}G(\bm{x},\bm{y})\rho(\bm{y})\operatorname{d\\!}s_{\bm{y}},\quad$
$\displaystyle\bm{x}$ $\displaystyle\in\Gamma_{1},$ (5)
$\displaystyle\mathcal{K}_{1}$ $\displaystyle\colon H^{1/2}(\Gamma_{1})\to
H^{1/2}(\Gamma_{1}),$ $\displaystyle\mathcal{K}_{1}\rho(\bm{x})$
$\displaystyle=\int_{\Gamma_{1}}\gamma^{1}_{n,\bm{y}}G(\bm{x},\bm{y})\rho(\bm{y})\operatorname{d\\!}s_{\bm{y}},$
$\displaystyle\bm{x}$ $\displaystyle\in\Gamma_{1},$

and the term $\tfrac{1}{2}I$ in Eq. 4 comes from the jump of the double layer
potential as $\Omega_{1}\ni\bm{x}\to\bm{y}\in\Gamma_{1}$. Problem Eq. 4 can be
rewritten as

$\mathcal{V}_{1}\gamma^{1}_{n}u_{1}=(\mathcal{K}_{1}+\tfrac{1}{2}I)\gamma^{1}_{t}u_{1},$
(6)

or employing the Poincaré–Steklov operator (Dirichlet-to-Neumann map):

$\mathcal{P}_{1}\mathrel{\mathop{\ordinarycolon}}H^{1/2}(\Gamma_{1})\rightarrow
H^{-1/2}(\Gamma_{1}),\qquad\mathcal{P}_{1}\coloneqq\mathcal{V}_{1}^{-1}(\mathcal{K}_{1}+\tfrac{1}{2}I),$
(7)

then Eq. 6 becomes

$\gamma^{1}_{n}u_{1}=\mathcal{P}_{1}\gamma^{1}_{t}u_{1}.$ (8)

The Poincaré–Steklov operator $\mathcal{P}_{1}$ in Eq. 7 is known to be
symmetric [35, Section 3.7]. Now, let $\gamma^{0}_{t}u_{0}$,
$\gamma^{0}_{n}u_{0}$ be the trace and conormal derivative of $u_{0}$ on
$\Gamma_{0}=\Gamma_{1}$, respectively. In order to derive a Dirichlet-to-
Neumann map
$\mathcal{P}_{0}\mathrel{\mathop{\ordinarycolon}}H^{1/2}(\Gamma_{0})\rightarrow
H^{-1/2}(\Gamma_{0})$ in $\Omega_{0}$, hence

$\gamma^{0}_{n}u_{0}=\mathcal{P}_{0}\gamma^{0}_{t}u_{0},$ (9)

we need to take into account the boundary condition (1g) on the external
boundary $\Sigma$ of $\Omega_{0}$. In order to alleviate the presentation we
postpone the derivation of $\mathcal{P}_{0}$ to A. However, we would like to
note that $\mathcal{P}_{0}$ remains symmetric.

Due to the Green representation formula, Eqs. 1a and 1b can be dropped from
Eq. 1. Also, Eq. 1g is encoded into the definition of $\mathcal{P}_{0}$ (see
A). Finally, the boundary integral formulation of Eq. 1 is

$\displaystyle\sigma_{1}\mathcal{P}_{1}\gamma^{1}_{t}u_{1}+\sigma_{0}\mathcal{P}_{0}\gamma^{0}_{t}u_{0}$
$\displaystyle=0,$ (10a)
$\displaystyle\gamma^{1}_{t}u_{1}-\gamma^{0}_{t}u_{0}$ $\displaystyle=V_{0},$
(10b) $\displaystyle\sigma_{0}\mathcal{P}_{0}\gamma^{0}_{t}u_{0}$
$\displaystyle=I_{\text{t}}(V_{0},z),$ (10c) $\displaystyle\partial_{t}z$
$\displaystyle=g(V_{0},z).$ (10d)

###  2.3 Spatial discretization of the unicellular problem

We adopt the collocation BEM as spatial discretization scheme. Boundary
element methods have less degrees of freedom than other standard techniques,
while the collocation approach yields lower dimensional boundary integrals
than the variational method and hence faster computations. For extensive
presentations on the BEM we refer to [26, 35, 38].

We place $M$ collocation points $\bm{x}_{j}$, $j=1,\ldots,M$, on $\Gamma_{1}$
in a counterclockwise order. Then we compute a smooth parametrization
$\gamma_{\Gamma_{1}}\mathrel{\mathop{\ordinarycolon}}[0,1)\rightarrow\Gamma_{1}$
satisfying

$\gamma_{\Gamma_{1}}(t_{j})=\bm{x}_{j}\quad j=1,\ldots,M,$ (11)

where $\\{t_{j}\\}_{j=1}^{M}\subset[0,1)$ is an increasing sequence. (The
parametrization $\gamma_{\Gamma_{1}}$ is computed with Fourier interpolation.
For the unicellular problem, we could define $\gamma_{\Gamma_{1}}(t)$ first
and then set $\bm{x}_{j}$ as in Eq. 11. However, this is not possible for the
multi cell problems.) Finally, we represent $\gamma^{1}_{t}u_{1}$,
$\gamma^{1}_{n}u_{1}$ as

$\displaystyle\gamma^{1}_{t}u_{1}(\gamma_{\Gamma_{1}}(t))$
$\displaystyle=\sum_{j=1}^{M}u_{1}^{j}L_{j}(t),$
$\displaystyle\gamma^{1}_{n}u_{1}(\gamma_{\Gamma_{1}}(t))$
$\displaystyle=\sum_{j=1}^{M}\tilde{u}_{1}^{j}L_{j}(t),$ (12)

where $L_{j}(t)$ are trigonometric Lagrange polynomials satisfying
$L_{j}(t_{k})=\delta_{jk}$ for $j,k=1,\ldots,M$. Instead of Eq. 6 we solve the
weaker form

$\mathcal{V}_{1}\gamma^{1}_{n}u_{1}(\bm{x}_{k})=(\mathcal{K}_{1}+\tfrac{1}{2}I)\gamma^{1}_{t}u_{1}(\bm{x}_{k}),\qquad
k=1,\ldots,M,$ (13)

with $\bm{x}_{k}=\gamma_{\Gamma_{1}}(t_{k})$, which is equivalent to

$\sum_{j=1}^{M}\tilde{u}_{1}^{j}\mathcal{V}_{1}L_{j}(\gamma_{\Gamma_{1}}^{-1}(\bm{x}_{k}))=\sum_{j=1}^{M}u_{1}^{j}(\mathcal{K}_{1}+\tfrac{1}{2}I)L_{j}(\gamma_{\Gamma_{1}}^{-1}(\bm{x}_{k})),\qquad
k=1,\ldots,M,$ (14)

and hence the linear system

$V_{1}\bm{\tilde{u}}_{1}=(K_{1}+\tfrac{1}{2}I)\bm{u}_{1},$ (15)

with $\bm{u}_{1},\bm{\tilde{u}}_{1}$ the vectors of coefficients
$u_{1}^{j},\tilde{u}_{1}^{j}$, respectively, and

$\displaystyle(V_{1})_{kj}$
$\displaystyle\coloneqq\mathcal{V}_{1}(L_{j}\circ\gamma_{\Gamma_{1}}^{-1})(\bm{x}_{k})=\int_{\Gamma_{1}}\gamma^{1}_{t,\bm{y}}G(\bm{x}_{k},\bm{y})L_{j}(\gamma_{\Gamma_{1}}^{-1}(\bm{y}))\operatorname{d\\!}s_{\bm{y}}=\int_{0}^{1}\gamma^{1}_{t,\bm{y}}G(\bm{x}_{k},\gamma_{\Gamma_{1}}(t))L_{j}(t)\|\gamma_{\Gamma_{1}}^{\prime}(t)\|\operatorname{d\\!}t,$
$\displaystyle(K_{1})_{kj}$
$\displaystyle\coloneqq\mathcal{K}_{1}(L_{j}\circ\gamma_{\Gamma_{1}}^{-1})(\bm{x}_{k})=\int_{\Gamma_{1}}\gamma^{1}_{n,\bm{y}}G(\bm{x}_{k},\bm{y})L_{j}(\gamma_{\Gamma_{1}}^{-1}(\bm{y}))\operatorname{d\\!}s_{\bm{y}}=\int_{0}^{1}\gamma^{1}_{n,\bm{y}}G(\bm{x}_{k},\gamma_{\Gamma_{1}}(t))L_{j}(t)\|\gamma_{\Gamma_{1}}^{\prime}(t)\|\operatorname{d\\!}t.$

The matrix coefficients $(V_{1})_{kj}$, $(K_{1})_{kj}$ must be computed with
special care due to the singularities in the fundamental solution
$G(\bm{x},\bm{y})$ and its derivatives as
$\gamma_{\Gamma_{1}}(t)\to\bm{x}_{k}$, we refer to [16, 26] for the details.

Note that $\bm{u}_{1}$ and $\bm{\tilde{u}}_{1}$ are the vectors of coordinates
of $\gamma^{1}_{t}u_{1}$ and $\gamma^{1}_{n}u_{1}$, respectively, and that
from Eq. 15 follow the discrete version of Eq. 8

$\displaystyle\bm{\tilde{u}}_{1}$ $\displaystyle=P_{1}\bm{u}_{1},$
$\displaystyle P_{1}$
$\displaystyle\coloneqq(V_{1})^{-1}(K_{1}+\tfrac{1}{2}I),$ (16)

where $P_{1}$ is the discrete Poincaré–Steklov operator (Dirichlet-to-Neumann
map) in $\Omega_{1}$. Similarly, in A we derive the discrete version of Eq. 9
and obtain

$\displaystyle\bm{\tilde{u}}_{0}$ $\displaystyle=P_{0}\bm{u}_{0},$ (17)

with $\bm{u}_{0}$ and $\bm{\tilde{u}}_{0}$ the vectors of coordinates of
$\gamma^{0}_{t}u_{0}$ and $\gamma^{0}_{n}u_{0}$, respectively.

Finally, the space discretization of boundary integral formulation Eq. 10 is

$\displaystyle\sigma_{1}P_{1}\bm{u}_{1}+\sigma_{0}P_{0}\bm{u}_{0}$
$\displaystyle=0,$ (18a) $\displaystyle\bm{u}_{1}-\bm{u}_{0}$
$\displaystyle=\bm{V}_{\text{0}},$ (18b)
$\displaystyle\sigma_{0}P_{0}\bm{u}_{0}$
$\displaystyle=I_{\text{t}}(\bm{V}_{\text{0}},\bm{z}),$ (18c)
$\displaystyle\bm{z}^{\prime}$ $\displaystyle=g(\bm{V}_{\text{0}},\bm{z}),$
(18d)

where $\bm{V}_{\text{0}}(t)\in\mathbb{R}^{M}$ is the vector whose coefficients
represent $V_{\text{0}}(\bm{x}_{j},t)$ and analogously for $\bm{z}$. The
right-hand sides $I_{\text{t}}$ and $g$ are applied to
$\bm{V}_{\text{0}},\bm{z}$ component wise. If needed, the solution $u_{1}$
satisfying Eq. 1a is approximated via the Green identity Eq. 3 and Eq. 12, Eq.
16. We proceed similarly for $u_{0}$.

###  2.4 The Lagrange multipliers approach for the unicellular problem

Now we solve Eq. 18 and to do so we employ the Lagrange multiplier method. We
adopted this technique mainly for pedagogical reasons in regard of what will
be presented in Section 3, since for the unicellular problem Eq. 1 a more
direct approach could be used.

In the remaining of this section we construct the linear map

$\displaystyle\mathcal{\psi}$
$\displaystyle\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{M}\rightarrow\mathbb{R}^{M},\qquad$
$\displaystyle\mathcal{\psi}(\bm{V}_{\text{0}})=\sigma_{0}P_{0}\bm{u}_{0},$
(19)

where $\bm{u}_{0}$ satisfies Eqs. 18a and 18b (behind the scenes $\bm{u}_{1}$
is computed as well, but it is not needed as output of $\mathcal{\psi}$).
Inserting Eq. 19 and
$I_{\text{t}}(\bm{V}_{\text{0}},\bm{z})=C_{\text{m}}\bm{V}_{\text{0}}^{\prime}+I_{\text{ion}}(\bm{V}_{\text{0}},\bm{z})$
in Eqs. 18c and 18d the problem reduces to the ODE

$\displaystyle
C_{\text{m}}\bm{V}_{\text{0}}^{\prime}+I_{\text{ion}}(\bm{V}_{\text{0}},\bm{z})$
$\displaystyle=\mathcal{\psi}(\bm{V}_{\text{0}}),$ (20a)
$\displaystyle\bm{z}^{\prime}$ $\displaystyle=g(\bm{V}_{\text{0}},\bm{z}),$
(20b)

which can be integrated by any suitable time marching scheme. Eq. 20 has the
same structure as the one derived in [20], where a BEM for the cell-by-cell
model without gap junctions is derived.

Theorem 2.1 is the unicellular version of the more general Theorem 3.1 below,
which in turn takes inspiration from the work in [27]. We also remark that the
Theorem is independent on the spatial discretization. For instance, a finite
element discretization may be recast to Eq. 18 by static condensation, that it
by explicitly computing the discrete Poincaré–Steklov operator.

######  Theorem 2.1.

The linear map $\mathcal{\psi}$ from Eq. 19 satisfies
$\mathcal{\psi}(\bm{V}_{\text{0}})=\bm{\lambda}$, with
$\bm{\lambda}\in\mathbb{R}^{M}$ and $\beta_{1}\in\mathbb{R}$ solution to

$\begin{pmatrix}F&G\\\ G^{\top}&0\end{pmatrix}\begin{pmatrix}\bm{\lambda}\\\
\beta_{1}\end{pmatrix}=\begin{pmatrix}\bm{V}_{\text{0}}\\\ 0\end{pmatrix}.$
(21)

The matrices $F\in\mathbb{R}^{M\times M}$, $G\in\mathbb{R}^{M}$ are defined by

$F=-(\sigma_{1}^{-1}(P_{1}^{+})^{-1}+\sigma_{0}^{-1}(P_{0}^{+})^{-1}),\qquad
G=\bm{e},$ (22)

$\bm{e}\in\mathbb{R}^{M}$ is the vector of ones and

$P_{0}^{+}=P_{0}+\alpha_{0}\bm{e}\bm{e}^{\top},\qquad
P_{1}^{+}=P_{1}+\alpha_{1}\bm{e}\bm{e}^{\top},$ (23)

with $\alpha_{0},\alpha_{1}>0$. If needed, $\bm{u}_{1},\bm{u}_{0}$ are
computed with

$\bm{u}_{0}=\sigma_{0}^{-1}(P_{0}^{+})^{-1}\bm{\lambda},\qquad\bm{u}_{1}=-\sigma_{1}^{-1}(P_{1}^{+})^{-1}\bm{\lambda}+\beta_{1}\bm{e}.$
(24)

Hence, when solving the ODE system Eq. 20 with a time integration scheme,
every time that $\mathcal{\psi}(\bm{V}_{\text{0}})$ needs to be evaluated
system Eq. 21 is solved and $\mathcal{\psi}(\bm{V}_{\text{0}})=\bm{\lambda}$
is inserted in Eq. 20.

###### Proof of Theorem 2.1..

Vectors $\bm{u}_{1}$, $\bm{u}_{0}$ are solutions to the smaller system in Eqs.
18a and 18b. Since $P_{0}$, $P_{1}$ are symmetric, $\bm{u}_{1}$, $\bm{u}_{0}$
are also solution to the constrained minimization problem

$\min_{\bm{u}_{1},\bm{u}_{0}}\frac{\sigma_{0}}{2}\langle
P_{0}\bm{u}_{0},\bm{u}_{0}\rangle+\frac{\sigma_{1}}{2}\langle
P_{1}\bm{u}_{1},\bm{u}_{1}\rangle\quad\text{with}\quad\bm{u}_{1}-\bm{u}_{0}=\bm{V}_{\text{0}},$
(25)

where $\langle\cdot,\cdot\rangle$ is the Euclidean inner product in
$\mathbb{R}^{M}$. Let

$\mathcal{L}(\bm{u}_{0},\bm{u}_{1},\bm{\lambda})=\frac{\sigma_{0}}{2}\langle
P_{0}\bm{u}_{0},\bm{u}_{0}\rangle+\frac{\sigma_{1}}{2}\langle
P_{1}\bm{u}_{1},\bm{u}_{1}\rangle+\langle\bm{u}_{1}-\bm{u}_{0}-\bm{V}_{\text{0}},\bm{\lambda}\rangle$
(26)

be the Lagrangian function, imposing
$\nabla\mathcal{L}(\bm{u}_{0},\bm{u}_{1},\bm{\lambda})=\bm{0}$ yields

$\sigma_{0}P_{0}\bm{u}_{0}-\bm{\lambda}=\bm{0},\qquad\sigma_{1}P_{1}\bm{u}_{1}+\bm{\lambda}=\bm{0},\qquad\bm{u}_{1}-\bm{u}_{0}=\bm{V}_{\text{0}}.$
(27)

Note that $P_{0}\bm{e}=\bm{0}$, with $\bm{e}\in\mathbb{R}^{M}$ a vector of
ones, thus $P_{0}$ is singular and the first equation of (27)
$\sigma_{0}P_{0}\bm{u}_{0}=\bm{\lambda}$ has a solution only if

$\langle\bm{\lambda},\bm{e}\rangle=0.$ (28)

Let

$P_{0}^{+}=P_{0}+\alpha_{0}\bm{e}\bm{e}^{\top}$ (29)

with $\alpha_{0}>0$, then $P_{0}^{+}$ is invertible and it can be verified
that if Eq. 28 holds then

$\bm{u}_{0}=\sigma_{0}^{-1}(P_{0}^{+})^{-1}\bm{\lambda}+\beta_{0}\bm{e},\qquad\beta_{0}\in\mathbb{R},$
(30)

is solution to $\sigma_{0}P_{0}\bm{u}_{0}=\bm{\lambda}$. Similarly,

$\bm{u}_{1}=-\sigma_{1}^{-1}(P_{1}^{+})^{-1}\bm{\lambda}+\beta_{1}\bm{e},\qquad\beta_{1}\in\mathbb{R},$
(31)

is solution to $\sigma_{1}P_{1}\bm{u}_{1}=-\bm{\lambda}$, with $P_{1}^{+}$
defined analogously to $P_{0}^{+}$. Note that if $\bm{u}_{0},\bm{u}_{1}$ are
solutions to Eqs. 18a and 18b then $\bm{u}_{0}+C\bm{e}$, $\bm{u}_{1}+C\bm{e}$
are solutions for all $C\in\mathbb{R}$. We choose to fix such free constant by
setting $\beta_{0}=0$, which implies $\sum_{j=1}^{M}u_{0}^{j}=0$. The last
equation of Eq. 27 yields

$\displaystyle\bm{V}_{\text{0}}$
$\displaystyle=\bm{u}_{1}-\bm{u}_{0}=-(\sigma_{1}^{-1}(P_{1}^{+})^{-1}+\sigma_{0}^{-1}(P_{0}^{+})^{-1})\bm{\lambda}+\beta_{1}\bm{e}=F\bm{\lambda}+G\beta_{1},$
(32)

with matrices $F$, $G$ as in Eq. 22. Together with Eq. 28 it yields system Eq.
21. ∎

######  Remark 2.2.

Note that the content of this section is readily adapted to a problem Eq. 1
with unbounded domain $\Omega_{0}$, hence without boundary condition (1g). In
that case, $P_{0}$ is derived analogously to $P_{1}$. However, $P_{0}$ would
be non singular hence in Theorem 2.1 we consider Eq. 29 with $\alpha_{0}=0$.
Condition Eq. 28 is still required for the existence of a solution to
$\sigma_{1}P_{1}\bm{u}_{1}=-\bm{\lambda}$.

##  3 Discretization of the full cell-by-cell model

We introduce here the general cell-by-cell model. We consider an extracellular
domain $\Omega_{0}\subset\mathbb{R}^{d}$, $d\geq 2$, an intracellular domain
$\Omega_{\mathrm{I}}\subset\mathbb{R}^{d}$, and an interface domain
$\Gamma_{0}=\bar{\Omega}_{0}\cap\bar{\Omega}_{\mathrm{I}}$. (See Figure 2 for
a graphical illustration of the model.) We suppose that $\Omega_{\mathrm{I}}$
and $\Omega_{0}$ are disjoint and we denote by $\Omega$ the whole tissue,
$\Omega=\Omega_{\mathrm{I}}\cup\Omega_{0}\cup\Gamma_{0}$. The domain $\Omega$
is always assumed connected and bounded, with $\Sigma=\partial\Omega$. For the
sake of simplicity, $\partial\Omega_{0}\setminus\Gamma_{0}=\Sigma$ and
$\partial\Omega_{\mathrm{I}}=\Gamma_{0}$, that is the exterior boundary of
$\Omega$ always corresponds to the extracellular matrix. Next, the
intracellular space is described by the union of disjoint cells, denoted by
$\Omega_{i}$, $i=1,\ldots,N$. (Conveniently, $\Omega_{i}$ for $i=0$
corresponds to the extracellular space.) Thus,
$\Omega_{\mathrm{I}}=\bigcup_{i=1}^{N}\Omega_{i}$. We denote
$\Gamma_{i}=\partial\Omega_{i}$, $i=1,\ldots,N$. The cell-to-cell
interconnections are denoted by $\Gamma_{ij}=\Gamma_{i}\cap\Gamma_{j}$, $0\leq
j<i\leq N$. Note that the boundary of each cell is either in contact with
another cell or with the extracellular space. The cell-by-cell model reads as
follows:

$\displaystyle-\sigma_{i}\Delta u_{i}$ $\displaystyle=0,$
$\displaystyle\text{in $\Omega_{i}$, $i=0,\ldots,N$},$ (33a) $\displaystyle
u_{i}-u_{0}$ $\displaystyle=V_{i0},$ $\displaystyle\text{on $\Gamma_{i0}$ for
$1\leq i\leq N$},$ (33b) $\displaystyle-\sigma_{i}\partial_{\bm{n}_{i}}u_{i}$
$\displaystyle=C_{\text{m}}\partial_{t}V_{i0}+I_{\text{ion}}(V_{i0},z_{i}),$
$\displaystyle\text{on $\Gamma_{i0}$ for $1\leq i\leq N$},$ (33c)
$\displaystyle-\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$
$\displaystyle=-C_{\text{m}}\partial_{t}V_{i0}-I_{\text{ion}}(V_{i0},z_{i}),\qquad$
$\displaystyle\text{on $\Gamma_{i0}$ for $1\leq i\leq N$},$ (33d)
$\displaystyle\partial_{t}z_{i}$ $\displaystyle=g(V_{i0},z_{i}),$
$\displaystyle\text{on $\Gamma_{0}$},$ (33e)
$\displaystyle-\sigma_{i}\partial_{\bm{n}_{i}}u_{i}$
$\displaystyle=\kappa(u_{i}-u_{j}),$ $\displaystyle\text{on $\Gamma_{ij}$ for
$1\leq j<i\leq N$},$ (33f)
$\displaystyle-\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$ $\displaystyle=0,\qquad$
$\displaystyle\text{on $\Sigma$}.$ (33g)

The constant conductivities are $\sigma_{i}$, $i=0,\ldots,N$. The gap
junctions (intercellular connections) are represented by $\Gamma_{ij}$ for
$1\leq j<i\leq N$, with permeability $\kappa$. The normals $\bm{n}_{i}$ point
outwards to $\Omega_{i}$.

Figure 2: Illustration of problem Eq. 33.

The intracellular potentials are $u_{i}$ for $i=1,\ldots,N$, the extracellular
potential is $u_{0}$ and $V_{i0}$ is the transmembrane potential on
$\Gamma_{i0}$. The membrane electric capacitance is $C_{\text{m}}$ and
$I_{\text{ion}}$ represents the sum of ionic currents. The transmembrane
potential $V_{i0}$ is regulated by the ionic currents, which in turn depend on
ionic concentrations and their transmembrane fluxes through ion channels,
which are governed by gating variables. Ion concentrations and gating
variables are represented by $z_{i}$ and the pair $I_{\text{ion}},g$ describe
the membrane ionic model. Several ionic models exist and they typically
consist of few to hundreds of equations. We remark that there is no
restriction in the system (33) for having different ionic models on each cell.

Model Eq. 33 is a slight simplification of a more detailed model by TMR21
[44], where the dynamics at the gap junctions is time dependent and nonlinear
in $u_{i}-u_{j}$. The simplification adopted here follows from linearization
and an equilibrium assumption. This procedure leads to a less computationally
intensive model. Solving the complete model and compare the results is subject
of a future work.

In this section we adapt the techniques used in Section 2 to the full problem
Eq. 33. First, in Section 3.1, we perform the spatial discretization of the
cell-by-cell model, obtaining a differential algebraic equation. Then, in
Section 3.2, we reduce the problem to a system of ordinary differential
equations on the transmembrane boundary.

We start rewriting Eq. 33 as follows. Let $V_{ij}=u_{i}-u_{j}$ be the
difference of potential defined on the gap junctions $\Gamma_{ij}$ for $1\leq
j<i\leq N$. Please note that we consider $\Gamma_{ij}$ with $j<i$ only, this
is to avoid any confusion regarding the sign of $u_{i}$ and $u_{j}$ in the
definition of $V_{ij}$. Let $\Gamma_{g}=\cup_{1\leq j<i\leq N}\Gamma_{ij}$ be
the union of all gap junctions, $\Gamma_{0}$ is the transmembrane boundary and
$\Gamma=\Gamma_{0}\cup\Gamma_{g}$ the union of all internal boundaries. We
define $V$ on $\Gamma$ by $V|_{\Gamma_{ij}}=V_{ij}$. This yields Eq. 34c
instead of Eq. 33b. Condition (33d) yields Eq. 34e, while summing Eqs. 33d and
33c we obtain Eq. 34b for $j=0$. Summing Eq. 33f inverting the roles of $i,j$
yields Eq. 34b for $j\geq 1$, while taking the difference gives Eq. 34d.

$\displaystyle-\sigma_{i}\Delta u_{i}$ $\displaystyle=0,$
$\displaystyle\text{in $\Omega_{i}$ for $i=0,\ldots,N$},$ (34a)
$\displaystyle\sigma_{i}\partial_{\bm{n}_{i}}u_{i}+\sigma_{j}\partial_{\bm{n}_{j}}u_{j}$
$\displaystyle=0,$ $\displaystyle\text{on $\Gamma_{ij}$ for $0\leq j<i\leq
N$},$ (34b) $\displaystyle u_{i}-u_{j}$ $\displaystyle=V,$
$\displaystyle\text{on $\Gamma_{ij}$ for $0\leq j<i\leq N$},$ (34c)
$\displaystyle\sigma_{j}\partial_{\bm{n}_{j}}u_{j}-\sigma_{i}\partial_{\bm{n}_{i}}u_{i}$
$\displaystyle=2\kappa V,$ $\displaystyle\text{on $\Gamma_{ij}$ for $1\leq
j<i\leq N$},$ (34d) $\displaystyle\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$
$\displaystyle=I_{\text{t}}(V,z),\qquad$ $\displaystyle\text{on
$\Gamma_{0}$},$ (34e) $\displaystyle\partial_{t}z$ $\displaystyle=g(V,z),$
$\displaystyle\text{on $\Gamma_{0}$}.$ (34f)
$\displaystyle\sigma_{0}\partial_{\bm{n}_{0}}u_{0}$ $\displaystyle=0,\qquad$
$\displaystyle\text{on $\Sigma$}.$ (34g)

Model Eq. 34 is equivalent to Eq. 33, however it is written in a more
“symmetric” manner.

###  3.1 Spatial discretization of the cell-by-cell model

We discretize all boundary segments $\Gamma_{ij}$ with $M_{ij}$ collocation
points $\bm{x}_{ij}^{k}\in\mathring{\Gamma}_{ij}$, for $k=1,\ldots,M_{ij}$ and
$0\leq j<i\leq N$. Let $M_{i}$ be the number of discretization points lying on
boundary $\Gamma_{i}$, $i=0,\ldots,N$. The total number of collocation points
on $\Gamma=\cup_{i=0}^{N}\Gamma_{i}$ is $M=\sum_{0\leq j<i\leq
N}M_{ij}=\tfrac{1}{2}\sum_{i=0}^{N}M_{i}$. We denote $\bm{x}^{l}$,
$l=1,\ldots,M$, the global collocations points on $\Gamma$ and by
$\bm{x}_{i}^{k}$, $k=1,\ldots,M_{i}$, the local collocation points on
$\Gamma_{i}$. Note that every $\bm{x}^{l}$ lies on some $\Gamma_{ij}$, hence
there are $\bm{x}_{ij}^{k_{1}}$, $\bm{x}_{i}^{k_{2}}$, $\bm{x}_{j}^{k_{3}}$
satisfying
$\bm{x}^{l}=\bm{x}_{ij}^{k_{1}}=\bm{x}_{i}^{k_{2}}=\bm{x}_{j}^{k_{3}}$.

Let $A_{i}\in\mathbb{R}^{M_{i}\times M}$ be the boolean connectivity matrix
mapping a vector $\bm{v}\in\mathbb{R}^{M}$ of global nodal values on $\Gamma$
to the vector $\bm{v}_{i}\in\mathbb{R}^{M_{i}}$ of local nodal values on
$\Gamma_{i}$. Every line of $A_{i}$ has exactly one non zero element:
$(A_{i})_{kl}=1$ for $k,l$ such that $\bm{x}^{l}=\bm{x}_{i}^{k}$. Note that
$A_{i}^{\top}$ maps local to global degrees of freedom. We also define
$B_{i}\in\mathbb{R}^{M_{i}\times M}$ having the same sparsity pattern as
$A_{i}$. Let $(B_{i})_{kl}$ be the only non zero element in the $k$-th line,
hence $\bm{x}^{l}=\bm{x}_{i}^{k}$. If $\bm{x}^{l}\in\Gamma_{ij}$ with $j<i$
then $(B_{i})_{kl}=1$, else $(B_{i})_{kl}=-1$.

Let $P_{i}\in\mathbb{R}^{M_{i}\times M_{i}}$ be the discrete Poincaré–Steklov
operator on each domain $\Omega_{i}$ and $\bm{u}_{i}\in\mathbb{R}^{M_{i}}$ the
vector of coordinates representing $\gamma^{\text{i}}_{0}u_{i}$. The vector of
coordinates $\bm{V}\in\mathbb{R}^{M}$ represents $V$ and
$\bm{V}_{\text{0}}=A_{0}\bm{V}$ represents $V|_{\Gamma_{0}}$. The spatial
discretization of Eqs. 34c and 34b is given by

$\sum_{i=0}^{N}\sigma_{i}A_{i}^{\top}P_{i}\bm{u}_{i}=\bm{0},\qquad\qquad\sum_{i=0}^{N}B_{i}^{\top}\bm{u}_{i}=\bm{V}.$
(35)

Recall that $A_{0}$ is the connectivity matrix mapping a global vector
$\bm{v}\in\mathbb{R}^{M}$ to a local vector $\bm{v}_{0}\in\mathbb{R}^{M_{0}}$
on the transmembrane boundary $\Gamma_{0}$. Let $M_{g}=M-M_{0}$ be the number
of points on the gap junctions $\Gamma_{g}$ and
$A_{g}\in\mathbb{R}^{M_{g}\times M}$ the matrix mapping a global vector to a
local vector $\bm{v}_{g}\in\mathbb{R}^{M_{g}}$ on $\Gamma_{g}$. The spatial
discretization of Eqs. 34e and 34d is

$\sigma_{0}P_{0}\bm{u}_{0}=I_{\text{t}}(A_{0}\bm{V},\bm{z}),\qquad\qquad\sum_{i=1}^{N}\sigma_{i}A_{g}B_{i}^{\top}P_{i}\bm{u}_{i}=-2\kappa
A_{g}\bm{V}.$ (36)

As in Section 2, conditions Eqs. 34a and 34g are automatically satisfied by
the Green representation formula (3) and the definition of the
Poincaré–Steklov operator $P_{0}$ on $\Omega_{0}$. Finally, the spatial
discretization of Eq. 34f is

$\bm{z}^{\prime}=g(A_{0}\bm{V},\bm{z}).$ (37)

Hence, the spatial discretization of Eq. 34 is given by Eqs. 35, 36 and 37.

###  3.2 Reduction to an ordinary differential equation

In this section we transform the space discretization Eqs. 35, 36 and 37 into
an ordinary differential equation. First, similarly to Section 2.4, we search
for linear maps

$\mathcal{\psi}_{i}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{M}\rightarrow\mathbb{R}^{M_{i}},\qquad\qquad\mathcal{\psi}_{i}(\bm{V})=\sigma_{i}P_{i}\bm{u}_{i},\qquad\qquad
i=0,\ldots,N,$ (38)

where the $\bm{u}_{i}$ satisfy Eq. 35. With the help of these maps we can
dispose of Eq. 35 by inserting Eq. 38 into Eq. 36 and obtain the system of
equations

$\mathcal{\psi}_{0}(\bm{V})=I_{\text{t}}(A_{0}\bm{V},\bm{z})=C_{\text{m}}A_{0}\bm{V}^{\prime}+I_{\text{ion}}(A_{0}\bm{V},\bm{z}),\qquad\qquad\sum_{i=1}^{N}A_{g}B_{i}^{\top}\mathcal{\psi}_{i}(\bm{V})=-2\kappa
A_{g}\bm{V}.$ (39)

However, Eq. 39 is a differential algebraic equation (DAE), which requires
more involved time marching schemes than a simple ODE. Therefore, departing
from the definition of the maps $\mathcal{\psi}_{i}$ given in Theorem 3.1, in
Theorem 3.2 we derive a new map which takes into account also the algebraic
condition (second equality in Eq. 39). This new map will allow us to derive an
ODE instead of a DAE.

We start with the theorem below, where we compute the maps
$\mathcal{\psi}_{i}$ of Eq. 38. The procedure adopted here is inspired from
[27], where a domain decomposition technique for the BEM is presented.

######  Theorem 3.1.

The linear maps $\mathcal{\psi}_{i}$ from Eq. 38 satisfy

$\mathcal{\psi}_{i}(\bm{V})=-B_{i}\bm{\lambda},$ (40)

with $\bm{\lambda}\in\mathbb{R}^{M}$ and $\bm{\beta}\in\mathbb{R}^{N}$
solution to

$\begin{pmatrix}F&G\\\ G^{\top}&0\end{pmatrix}\begin{pmatrix}\bm{\lambda}\\\
\bm{\beta}\end{pmatrix}=\begin{pmatrix}\bm{V}\\\ \bm{0}\end{pmatrix}.$ (41)

The matrices $F\in\mathbb{R}^{M\times M}$, $G\in\mathbb{R}^{M\times N}$ are
defined by

$F=-\sum_{i=0}^{N}\sigma_{i}^{-1}B_{i}^{\top}(P_{i}^{+})^{-1}B_{i},\qquad
G=(B_{1}^{\top}\bm{e}_{1},\ldots,B_{N}^{\top}\bm{e}_{N}),$ (42)

$\bm{e}_{i}\in\mathbb{R}^{M_{i}}$ is the vector of ones and

$P_{i}^{+}=P_{i}+\alpha_{i}\bm{e}_{i}\bm{e}_{i}^{\top},$ (43)

with $\alpha_{i}>0$, $i=0,\ldots,N$. If needed, $\bm{u}_{i}$ for
$i=0,\ldots,N$ is computed with

$\bm{u}_{i}=-\sigma_{i}^{-1}(P_{i}^{+})^{-1}B_{i}\bm{\lambda}+\beta_{i}\bm{e}_{i},$
(44)

where $\bm{\beta}=(\beta_{1},\ldots,\beta_{N})^{\top}$ and $\beta_{0}=0$.

###### Proof.

As in Section 2.4 we notice that instead of solving Eq. 35 we can solve a
constrained minimization problem with Lagrangian function

$\mathcal{L}(\bm{u}_{0},\ldots,\bm{u}_{N},\bm{\lambda})=\sum_{i=0}^{N}\frac{\sigma_{i}}{2}\langle
P_{i}\bm{u}_{i},\bm{u}_{i}\rangle+\sum_{i=0}^{N}\langle
B_{i}^{\top}\bm{u}_{i},\bm{\lambda}\rangle-\langle\bm{V},\bm{\lambda}\rangle.$
(45)

Indeed,

$\nabla_{\bm{\lambda}}\mathcal{L}(\bm{u}_{0},\ldots,\bm{u}_{N},\bm{\lambda})=\sum_{i=0}^{N}B_{i}^{\top}\bm{u}_{i}-\bm{V}=\bm{0}$
(46)

is equivalent to the second equality in Eq. 35. The first equality in Eq. 35
follows from

$\nabla_{\bm{u}_{i}}\mathcal{L}(\bm{u}_{0},\ldots,\bm{u}_{N},\bm{\lambda})=\sigma_{i}P_{i}\bm{u}_{i}+B_{i}\bm{\lambda}=\bm{0}$
(47)

and

$\bm{0}=\sum_{i=0}^{N}A_{i}^{\top}\nabla_{\bm{u}_{i}}\mathcal{L}(\bm{u}_{0},\ldots,\bm{u}_{N},\bm{\lambda})=\sum_{i=0}^{N}\sigma_{i}A_{i}^{\top}P_{i}\bm{u}_{i}+A_{i}^{\top}B_{i}\bm{\lambda}=\sum_{i=0}^{N}\sigma_{i}A_{i}^{\top}P_{i}\bm{u}_{i},$
(48)

where we used $\sum_{i=0}^{N}A_{i}^{\top}B_{i}=0$ (for every $1$ there is a
$-1$).

Let us solve Eq. 47. We denote $\bm{e}_{i}\in\mathbb{R}^{M_{i}}$ the vector of
ones, imposing

$\langle B_{i}\bm{\lambda},\bm{e}_{i}\rangle=0,\qquad i=0,\ldots,N,$ (49)

a solution $\bm{u}_{i}$ to Eq. 47 exists and is given by

$\bm{u}_{i}=-\sigma_{i}^{-1}(P_{i}^{+})^{-1}B_{i}\bm{\lambda}+\beta_{i}\bm{e}_{i},\quad\text{with}\quad
P_{i}^{+}=P_{i}+\alpha_{i}\bm{e}_{i}\bm{e}_{i}^{\top},$ (50)

$\alpha_{i}>0$ and $\beta_{i}\in\mathbb{R}$. Since if $\bm{u}_{i}$,
$i=0,\ldots,N$, are solutions to Eq. 35 then also $\bm{u}_{i}+C\bm{e}_{i}$ are
solutions, we choose to set $\beta_{0}=0$ and remove this degree of freedom.
Note as well that

$\langle B_{0}\bm{\lambda},\bm{e}_{0}\rangle=-\sum_{i=1}^{N}\langle
B_{i}\bm{\lambda},\bm{e}_{i}\rangle\qquad\forall\bm{\lambda}\in\mathbb{R}^{M},$
(51)

hence Eq. 49 is replaced by the sufficient one

$\langle B_{i}\bm{\lambda},\bm{e}_{i}\rangle=0,\qquad i=1,\ldots,N.$ (52)

Inserting Eq. 50 into Eq. 46 yields

$\bm{V}=\sum_{i=0}^{N}B_{i}^{\top}(-\sigma_{i}^{-1}(P_{i}^{+})^{-1}B_{i}\bm{\lambda}+\beta_{i}\bm{e}_{i})=-\sum_{i=0}^{N}\sigma_{i}^{-1}B_{i}^{\top}(P_{i}^{+})^{-1}B_{i}\bm{\lambda}+\sum_{i=1}^{N}\beta_{i}B_{i}^{\top}\bm{e}_{i},$
(53)

recall that $\beta_{0}=0$. From Eqs. 53 and 52 follow Eq. 41. ∎

Now we use the result of Theorem 3.1 and the second equality of Eq. 39 in
order to derive a standard ODE problem. We recall that
$\bm{V}_{\text{0}}=A_{0}\bm{V}$.

######  Theorem 3.2.

The space discretization Eqs. 35, 36 and 37 of Eq. 33 is equivalent to the
ordinary differential equations system

$\displaystyle
C_{\text{m}}\bm{V}_{\text{0}}^{\prime}+I_{\text{ion}}(\bm{V}_{\text{0}},\bm{z})$
$\displaystyle=\mathcal{\psi}(\bm{V}_{\text{0}}),$ (54a)
$\displaystyle\bm{z}^{\prime}$ $\displaystyle=g(\bm{V}_{\text{0}},\bm{z}),$
(54b)

where $\mathcal{\psi}(\bm{V}_{\text{0}})=\bm{\lambda}_{0}$ and
$\bm{\lambda}_{0}\in\mathbb{R}^{M_{0}}$,
$\bm{\lambda}_{g}\in\mathbb{R}^{M_{g}}$, $\bm{\beta}\in\mathbb{R}^{N}$ are
solutions to

$\begin{pmatrix}F_{00}&F_{0g}&A_{0}G\\\ F_{g0}&F_{gg}-\kappa^{-1}I&A_{g}G\\\
G^{\top}A_{0}^{\top}&G^{\top}A_{g}^{\top}&0\end{pmatrix}\begin{pmatrix}\bm{\lambda}_{0}\\\
\bm{\lambda}_{g}\\\
\bm{\beta}\end{pmatrix}=\begin{pmatrix}\bm{V}_{\text{0}}\\\ \bm{0}\\\
\bm{0}\end{pmatrix},$ (55)

with

$F_{00}=A_{0}FA^{\top}_{0},\quad F_{0g}=A_{0}FA^{\top}_{g},\quad
F_{g0}=A_{g}FA^{\top}_{0},\quad F_{gg}=A_{g}FA^{\top}_{g}.$ (56)

###### Proof.

We denote $\bm{V}_{\text{g}}=A_{g}\bm{V}$, $\bm{V}_{\text{0}}=A_{0}\bm{V}$,
$\bm{\lambda}_{g}=A_{g}\bm{\lambda}$ and $\bm{\lambda}_{0}=A_{0}\bm{\lambda}$.
Also, note that $B_{i}^{\top}B_{i}=A_{i}^{\top}A_{i}$ since
$B_{i}^{\top}B_{i}$ projects a global vector forth and back from $\Gamma_{i}$
and if a sign change happens it occurs twice. Therefore

$\left(\sum_{i=0}^{N}B_{i}^{\top}B_{i}\right)\bm{v}=\left(\sum_{i=0}^{N}A_{i}^{\top}A_{i}\right)\bm{v}=2\bm{v}\qquad\forall\,\bm{v}\in\mathbb{R}^{M},$
(57)

indeed every segment $\Gamma_{ij}$ will receive the contribution from exactly
two neighbouring domains. From Theorem 3.1 we have that
$\mathcal{\psi}_{i}(\bm{V})=-B_{i}\bm{\lambda}$, which inserted into the
second equality of Eq. 39 yields

$\kappa\bm{V}_{\text{g}}=\kappa
A_{g}\bm{V}=-\frac{1}{2}\sum_{i=1}^{N}A_{g}B_{i}^{\top}\mathcal{\psi}_{i}(\bm{V})=\frac{1}{2}A_{g}\left(\sum_{i=0}^{N}B_{i}^{\top}B_{i}\right)\bm{\lambda}=A_{g}\bm{\lambda}=\bm{\lambda}_{g}.$
(58)

Note as well that $A_{0}^{\top}A_{0}+A_{g}^{\top}A_{g}$ is the identity matrix
in $\mathbb{R}^{M}$, hence multiplying the first line
$F\bm{\lambda}+G\bm{\beta}=\bm{V}$ of Eq. 41 with $A_{0}$ yields

$\displaystyle\bm{V}_{\text{0}}$
$\displaystyle=A_{0}\bm{V}=A_{0}F(A_{0}^{\top}A_{0}\bm{\lambda}+A_{g}^{\top}A_{g}\bm{\lambda})+A_{0}G\bm{\beta}=F_{00}\bm{\lambda}_{0}+F_{0g}\bm{\lambda}_{g}+A_{0}G\bm{\beta}.$
(59)

Similarly, multiplication by $A_{g}$ yields
$\bm{V}_{\text{g}}=F_{g0}\bm{\lambda}_{0}+F_{gg}\bm{\lambda}_{g}+A_{g}G\bm{\beta}$
and thus

$\bm{0}=F_{g0}\bm{\lambda}_{0}+(F_{gg}-\kappa^{-1}I)\bm{\lambda}_{g}+A_{g}G\bm{\beta}.$
(60)

For the second line of Eq. 41 we have

$\bm{0}=G^{\top}\bm{\lambda}=G^{\top}A_{0}^{\top}\bm{\lambda}_{0}+G^{\top}A_{g}^{\top}\bm{\lambda}_{g}.$
(61)

Relations Eqs. 59, 60 and 61 yield Eq. 55. The identity $-B_{0}=A_{0}$ and
$\mathcal{\psi}_{0}(\bm{V})=-B_{0}\bm{\lambda}=A_{0}\bm{\lambda}=\bm{\lambda}_{0}$
implies $\mathcal{\psi}_{0}(\bm{V})=\mathcal{\psi}(\bm{V}_{0})$ and hence Eq.
54. ∎

From the proof of Theorem 3.2 we see that $\bm{V}$, $\bm{\lambda}$ of Theorem
3.1 are given by
$\bm{V}=A_{0}^{\top}\bm{V}_{0}+A_{g}^{\top}\bm{V}_{\text{g}}$,
$\bm{\lambda}=A_{0}^{\top}\bm{\lambda}_{0}+A_{g}^{\top}\bm{\lambda}_{g}$,
$\kappa\bm{V}_{\text{g}}=\bm{\lambda}_{g}$ and moreover $\bm{\beta}$ is the
same as in Theorem 3.1; hence, if needed, $\bm{u}_{i}$ for $i=0,\ldots,N$ can
be computed as in Theorem 3.1. Note as well that in Eq. 55 we have chosen to
use $\bm{\lambda}_{g}$ as unknown, instead of the alternative
$\bm{V}_{\text{g}}$. If we used $\bm{V}_{\text{g}}$ we would obtain the same
matrix as in Eq. 55 but with the second column multiplied by $\kappa$ and
therefore break the symmetry.

###  3.3 Time integration

For the time integration of Eq. 54, we use the multirate explicit stabilized
method mRKC [1] for problems

$\bm{y}^{\prime}=f_{F}(t,\bm{y})+f_{S}(t,\bm{y}),\qquad\bm{y}(0)=\bm{y}_{0},$
(62)

where $f_{F}$ is a stiff term and $f_{S}$ is a mildly stiff but more expensive
term. The mRKC scheme is fully explicit and does not have any step size
restriction. Its stability properties are inherited from the RKC methods [13],
which use an increased number of stages, with respect to classical methods, to
increase stability. Since stability grows quadratically with the work load,
the methods are particularly efficient. For the integration of Eq. 54 with
mRKC, we rewrite Eq. 54 as Eq. 62, with

$\bm{y}=\begin{pmatrix}\bm{V}_{\text{0}}\\\ \bm{z}\end{pmatrix},\qquad
f_{F}(t,\bm{y})=\begin{pmatrix}\mathcal{\psi}(\bm{V}_{\text{0}})/C_{\text{m}}\\\
\bm{0}\end{pmatrix},\qquad
f_{S}(t,\bm{y})=\begin{pmatrix}-(I_{\text{ion}}(\bm{V}_{\text{0}},\bm{z})+I_{\text{stim}}(t))/C_{\text{m}}\\\
g(\bm{V}_{\text{0}},\bm{z})\end{pmatrix},$ (63)

where $I_{\text{stim}}(t)$ is used to stimulate some cells and initiate an
action potential propagation.

##  4 Numerical experiments

In this section we perform some numerical experiments in order to asses the
accuracy of the space-time discretization of the cell-by-cell model (33) but
also investigate the regularity properties of the model itself.

We start with two experiments, in Sections 4.1 and 4.2, where we investigate
the convergence rates of the maps $\mathcal{\psi}_{i}$ from Theorem 3.1 and
then the impact of the mesh and step size on the accuracy of the conduction
velocity (CV). These experiments are crucial to understand which
discretization parameters yield solutions within a certain error tolerance.

In the subsequent experiments the goal is to study the model itself. For
instance, in Section 4.3 we investigate the effect of the gap junction’s
permeability $\kappa$ on the CV and in Section 4.4 we study how the contact
area between cells affects CV. Before presenting the results, we resume here
below our computational setting.

#### Computational setup

The following numerical experiments have been performed with our C++ code,
where for the dense linear algebra routines we employ the Eigen library [17].
The ionic model, is taken from CellML [11] and the relative C code is produced
with the Myokit library [9]. Concerning the model Eq. 33, the number of cells
$N$ and the domains $\Omega_{i}$ vary from one experiment to another and are
specified later. If not stated otherwise, in the next experiments we use the
coefficients $C_{\text{m}}$, $\sigma_{i}$, $\kappa$, given in Table 1. The
values for $C_{\text{m}}$, $\sigma_{i}$, $\kappa$ are taken from [40], where
for $\kappa$ we consider $\kappa=1/R_{m}$ with $R_{m}=0.00145\
$\mathrm{k\SIUnitSymbolOhm}\text{\,}{\mathrm{cm}}^{2}$$. If not specified, we
consider the ionic model from CRN98 [10]. The initial values for $V$ for and
the ionic model’s state variables are uniform on the transmembrane boundary
and are taken from the Myokit’s code. For instance, for the Courtemanche-
Nattel-Ramiréz model the initial value for $V$ is $V_{0}=-81.18\
$\mathrm{mV}$$.

$C_{m}$ | $\sigma_{0}$ | $\sigma_{1},\ldots,\sigma_{N}$ | $\kappa$
---|---|---|---
$1\text{\,}\mathrm{\SIUnitSymbolMicro F}\text{\,}{\mathrm{cm}}^{-2}$ | $20\text{\,}\mathrm{mS}\text{\,}{\mathrm{cm}}^{-1}$ | $3\text{\,}\mathrm{mS}\text{\,}{\mathrm{cm}}^{-1}$ | $690\text{\,}\mathrm{mS}\text{\,}{\mathrm{cm}}^{-2}$
Table 1: Model’s coefficients employed in numerical experiments.

###  4.1 Approximation properties of the $\mathcal{\psi}$ operators

In this experiment we investigate the convergence rate of the
$\mathcal{\psi}_{i}$ operators defined in Theorem 3.1. For this purpose, we
conduct four convergence experiments, one for every geometry displayed in Fig.
3, and we display the errors on the trace and normal derivative against the
number of degrees of freedom $M$ in Fig. 4.

../images/conv_psi/illustrations/circle

(a) Single cell

../images/conv_psi/illustrations/splitted_circle

(b) Two cells

../images/conv_psi/illustrations/two_semicircles

(c) Isolated, non-smooth cells

../images/conv_psi/illustrations/smoothed_semicircles

(d) Isolated, smooth cells

Figure 3: Illustration of the geometrical settings employed in Section 4.1.

Let us describe the geometries of Fig. 3. In Fig. 3(a) we have a model
composed of one cell ($N=1$), defined by
$\Omega_{0}=\\{\bm{x}\in\mathbb{R}^{2}\mathrel{\mathop{\ordinarycolon}}2<\|\bm{x}\|<4\\}$
and
$\Omega_{1}=\\{\bm{x}\in\mathbb{R}^{2}\mathrel{\mathop{\ordinarycolon}}\|\bm{x}\|<2\\}$.
In Fig. 3(b) we split the cell by introducing a vertical gap junction, hence
we have the same $\Omega_{0}$ but
$\Omega_{1}=\\{\bm{x}\in\mathbb{R}^{2}\mathrel{\mathop{\ordinarycolon}}\|\bm{x}\|<2,x_{1}<0\\}$
and
$\Omega_{2}=\\{\bm{x}\in\mathbb{R}^{2}\mathrel{\mathop{\ordinarycolon}}\|\bm{x}\|<2,x_{1}>0\\}$.
In Fig. 3(c) we keep the same cells but remove the gap junction by introducing
an horizontal gap of size $0.4$ between $\Omega_{1}$ and $\Omega_{2}$.
Finally, in Fig. 3(d) we keep the separation of $\Omega_{1}$, $\Omega_{2}$ but
smooth out the corners by introducing quarter of circles of radius
$0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

For the setting of Figs. 3(a) and 3(b) an exact solution to Eqs. 34a, 34c, 34b
and 34d, with $V$ defined by Eq. 34c, is given by

$\displaystyle u_{0}(\bm{x})$
$\displaystyle=\frac{\sigma_{1}}{\sigma_{0}}\frac{16+\|\bm{x}\|^{2}}{6\|\bm{x}\|^{2}}x_{2},\qquad$
$\displaystyle u_{1}(\bm{x})=u_{2}(\bm{x})$ $\displaystyle=-\frac{1}{2}x_{2}.$
(64)

Therefore, for different values of $M$ (i.e. number of collocation points), we
can compute the vector of coefficients $\bm{u}_{i}$, define $\bm{V}$ as in Eq.
35, solve Eq. 41 and compute the errors

$\displaystyle e_{1}$
$\displaystyle=\max_{i=0,\ldots,N}\|\mathcal{\psi}_{i}(\bm{V})-\sigma_{i}\partial_{\bm{n}_{i}}u_{i}\|_{L^{2}(\Gamma_{i})},\qquad$
$\displaystyle e_{0}$
$\displaystyle=\max_{i=0,\ldots,N}\|\widetilde{\mathcal{\psi}}_{i}(\bm{V})-u_{i}\|_{L^{2}(\Gamma_{i})/\mathbb{R}},$
(65)

where
$\widetilde{\mathcal{\psi}}_{i}(\bm{V})=-\sigma_{i}^{-1}(P_{i}^{+})^{-1}B_{i}\bm{\lambda}+\beta_{i}\bm{e}_{i}$
and thus $\widetilde{\mathcal{\psi}}_{i}(\bm{V}_{m})$ approximates $u_{i}$
(cf. Eq. 44), up to a constant.

For the geometries of Figs. 3(c) and 3(d) we do not possess an exact solution.
Hence, we set $V(\bm{x})=\cos(\pi x_{1})\sin(\pi x_{2})$ and errors $e_{1}$,
$e_{0}$ are now computed as

$\displaystyle e_{1}$
$\displaystyle=\max_{i=0,\ldots,N}\|\mathcal{\psi}_{i}(\bm{V})-\mathcal{\psi}_{i}(\bm{V}^{*})\|_{L^{2}(\Gamma_{i})},\qquad$
$\displaystyle e_{0}$
$\displaystyle=\max_{i=0,\ldots,N}\|\widetilde{\mathcal{\psi}}_{i}(\bm{V})-\widetilde{\mathcal{\psi}}_{i}(\bm{V}^{*})\|_{L^{2}(\Gamma_{i})/\mathbb{R}},$
(66)

where $\mathcal{\psi}_{i}(\bm{V}^{*})$,
$\widetilde{\mathcal{\psi}}_{i}(\bm{V}^{*})$ are reference solutions
calculated on a finer mesh.

../images/conv_psi/plots/res_all

Figure 4: Convergence rates of the $\mathcal{\psi}$ operator defined in
Theorem 3.1 for the problems depicted in Fig. 3.

We display the errors $e_{0},e_{1}$ with respect to $M$, for the geometries of
Fig. 3, in Fig. 4. Due to the smoothness of the solutions and the boundaries,
we remark in LABEL:fig:plot_conv_psi_a that for the first problem we obtain
exponential convergence thanks to the trigonometric Lagrange basis functions;
with very few degrees of freedom machine precision is achieved. This result is
in line with the theory and experiments performed in [19]. In the second
problem, the boundary is Lipschitz continuous only, which prevents exponential
convergence. Indeed, in LABEL:fig:plot_conv_psi_b the convergence rates for
the trace and normal derivative are $1.5$ and $0.5$, respectively.

In LABEL:fig:plot_conv_psi_a and LABEL:fig:plot_conv_psi_b we see how the
convergence rates decrease dramatically when the circle is divided into two
half-circles. The purpose of two last experiments is to demonstrate
numerically that this phenomenon is due to the non smooth boundaries, rather
than the introduction of a gap junction. Indeed, in LABEL:fig:plot_conv_psi_c
we observe the same convergence rates as in LABEL:fig:plot_conv_psi_b, while
in LABEL:fig:plot_conv_psi_d we obtain higher convergence rates.

###  4.2 Impact of discretization parameters on conduction velocity

It is already known that discretization methods, mesh size and step size
affect conduction velocity (CV) in the monodomain and bidomain models for
cardiac electrophysiology [8, 32, 33]. In this experiment we investigate how
mesh and step size affect the CV for the cell-by-cell model discretized with
the BEM in space and the mRKC method [1] in time. In order to be able to
employ relatively uniform mesh sizes in this experiment we consider
rectangular cells.

To measure the CV we design the following experiment. We consider an array of
$2\times 30$ connected rectangular cells of width
$c_{w}=$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and length
$c_{l}=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, cells are positioned so
that their bottom left vertex has coordinates $(i\cdot c_{l},j\cdot c_{w})$
for $i=0,\ldots,29$, $j=0,1$; yielding a block of cells of width $2\cdot
c_{w}=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and length $30\cdot
c_{l}=$3000\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The outer domain is a
rectangle of size $$440\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$5000\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ centred on the array of
cells. To initiate an action potential traversing the cell’s array a stimulus
of $300\text{\,}\mathrm{\SIUnitSymbolMicro A}\text{\,}{\mathrm{cm}}^{-2}$ is
applied for a duration of $1\text{\,}\mathrm{ms}$ at the transmembrane
boundary of the two leftmost cells. CV is computed as the average over
$CV_{k}=\|\bm{p}_{k}-\bm{q}_{k}\|/(t_{\bm{p}_{k}}-t_{\bm{q}_{k}})$ for
$k=1,\ldots,5$, where: $\bm{p}_{k}=((7.5+k)\cdot
c_{l})\,$\mathrm{\SIUnitSymbolMicro m}$$, $\bm{q}_{k}=((17.5+k)\cdot
c_{l},0)\,$\mathrm{\SIUnitSymbolMicro m}$$ and $t_{\bm{p}_{k}}$,
$t_{\bm{q}_{k}}$ are the time instants in which $V$ exceeds the threshold of
$V_{\mathrm{th}}=$-20\text{\,}\mathrm{mV}$$ in $\bm{p}_{k}$, $\bm{q}_{k}$,
respectively. The choice of $\bm{p}_{k}$, $\bm{q}_{k}$ is such that measures
are taken sufficiently far from the stimulated point and to avoid boundary
effects as well. See Fig. 5 for an illustration of the solution at time
$t=$2\text{\,}\mathrm{ms}$$.

Figure 5: Action potential propagation in an array of $2\times 20$ cells
($t=2$\mathrm{ms}$$).

We solve Eq. 54 with different step sizes $\Delta t$ and mesh size $\Delta x$
and compute the signed relative error on CV:
$E_{\mathrm{CV}}=(\mathrm{CV}-\mathrm{CV}^{*})/\mathrm{CV}^{*}$, with
$\mathrm{CV}^{*}$ a reference solution. We display $E_{\mathrm{CV}}$ as
function of $\Delta t$, $\Delta x$ in Fig. 6. The reference value of CV is
$\mathrm{CV}^{*}\approx 1.27153$.

../images/dep_cv_on_dtdx/CV_vs_dt_vs_dx

Figure 6: Conduction velocity accuracy with respect to timestep $\Delta t$ and
mesh size $\Delta x$.

First, we notice that for coarse space grids the true $\mathrm{CV}^{*}$ tends
to be overestimated, whereas for large time steps it is underestimated. Then,
we remark that even with relatively large mesh sizes $\Delta
x=$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ the estimated CV remains within
a $2\%$ error. Figure 6 (middle and right panels) also shows the same results,
but for the unsigned relative error $|E_{\mathrm{CV}}|$ and fixing either
$\Delta x$ or $\Delta t$. We observe that the local minimal appearing in the
curves is due to the cancellation of the positive spatial discretization error
with the negative time discretization error.

Based on the results of this section, in the forthcoming experiments we
consider $\Delta t\leq$0.02\text{\,}\mathrm{ms}$$ and $\Delta
x\leq$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, which, for this experiment,
yield a relative error of less than $5\%$.

###  4.3 Dependence of conduction velocity on gap junctions’ permeability and
cells inner conductivity

In this experiment we study how CV depends on the gap junctions’ permeability
$\kappa$ and the inner conductivity $\sigma_{i}$, $i=1,\ldots,N$. For that
purpose, we consider again an array of $2\times 30$ cells and cells of size
$c_{w}\times c_{l}$, with fixed $c_{w}=$10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$ and either $c_{l}=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ or
$c_{l}=$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. First, we measure CV for
varying $\kappa$ but keeping the other coefficients fixed, a stimulus is
initiated applying a stimulus of $200\text{\,}\mathrm{\SIUnitSymbolMicro
A}\text{\,}{\mathrm{cm}}^{-2}$ to the transmembrane boundary of the two
leftmost cells. Results for $c_{l}=$100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$ and $c_{l}=$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ are displayed in
LABEL:fig:cv_vs_kappa. We observe as CV decreases with $\kappa$ and also that
the physiological value $\kappa=690$ is in the range where CV is maximal. We
note that for values of $\kappa\leq 2\cdot 10^{-4}$ the action potential does
not propagate. Then, we measure CV for varying $\sigma_{i}$, $i=1,\ldots,N$,
and fixed $\kappa$, results are displayed in LABEL:fig:cv_vs_si for
$c_{l}=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$c_{l}=$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. Conduction velocity
increases with $\sigma_{i}$, specially for the shorter cells.

../images/dep_cv_on_kappa_si/CV_all

Figure 7: Dependence of CV on permeability $\kappa$ and intracellular
conductivity $\sigma_{i}$ (red: reference value).

###  4.4 Dependence of conduction velocity on gap junctions’ surface area

In general, gap junctions perpendicular to the fiber direction are not flat
surfaces and are better modelled by intercalated discs [21]. In two dimensions
we model these gap junctions with a sinus wave (see Fig. 8(a) for an
illustration) of amplitude $a$ and frequency $k$.

In this experiment we consider and array of $2\times 30$ cells of size
$$10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and measure CV as in
Section 4.2. In Fig. 8(b) we show the conductive velocity as a function of the
frequency $k$ for a fixed amplitude of
$a=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. We note that for moderate
frequency $k$ the CV increases due to an increase of contact surface area.
However, for larger frequencies CV decreases, probably because of a flux
saturation at the narrower junctions. In Fig. 8(c) we show the conductive
velocity as a function of the amplitude $a$ for a fixed frequency of $k=3$.
Again, for larger amplitude $a$ the conduction velocity decreases.

../images/dep_cv_on_area/intercalated_discs

(a) Two cells with an intercalated discs at gap junctions.

../images/dep_cv_on_area/CV_vs_freq

(b) Impact of intercalated discs frequency on CV for
$a=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.

../images/dep_cv_on_area/CV_vs_amp

(c) Impact of intercalated discs amplitude on CV for $k=3$.

Figure 8: Effect of gap junctions surface area on CV.

###  4.5 Dependence of conduction velocity on cells size and aspect ratio

Finally, we investigate how the cell’s size and aspect ratio impact the
conduction velocity. We consider an array of $2\times 30$ cells of size
$c_{w}\times c_{l}$. First, we fix
$c_{w}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and vary $c_{l}$, results
are reported in LABEL:fig:cv_vs_cl. We observe as CV decreases as $c_{l}$
increases. In LABEL:fig:cv_vs_cw we display the results for fixed
$c_{l}=$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and varying $c_{w}$, here
CV increases with $c_{w}$. In the last figure LABEL:fig:cv_vs_clcw we vary
both $c_{l}$ and $c_{w}$ while keeping a constant aspect ratio $c_{l}=10\cdot
c_{w}$, more precisely they vary from
$(c_{w},c_{l})=($2.5\text{\,}\mathrm{\SIUnitSymbolMicro
m}$,$25\text{\,}\mathrm{\SIUnitSymbolMicro m}$)$ to
$(c_{w},c_{l})=($14\text{\,}\mathrm{\SIUnitSymbolMicro
m}$,$140\text{\,}\mathrm{\SIUnitSymbolMicro m}$)$. We see as CV increases with
the cells area $A=c_{l}\cdot c_{w}$.

../images/dep_cv_on_cells_size/CV_vs_all

Figure 9: Impact of cell length $c_{l}$, cell width $c_{w}$, and cell area
with fixed aspect ratio on the CV.

##  5 Conclusion

In this paper we solve the cell-by-cell or EMI model for cardiac
electrophysiology via the boundary element method, with no geometrical
restrictions. The cell-by-cell model consists in Laplace equations inside and
outside the cellular domains coupled with an ordinary differential equation on
the transmembrane boundary and an algebraic condition on the gap junctions.
Due to the boundary integral formulation, Laplace equations are cleared away,
yielding a differential algebraic equation living on the cell’s boundaries
only. In a subsequent step, the differential algebraic equation is reduced to
an ordinary differential equation lying on the transmembrane boundary only.
Finally, we provide numerical results where: first, we study the accuracy of
the numerical method and then we investigate the model properties and
sensibility with respect to its parameters.

The convergence rate shows that the solution of the general problem is non-
smooth, due to the presence of multiple cell contact (or 2 cells and the
extracellular domain.) The single cell problem is instead smooth. To the best
of our knowledge, there are no regularity results for the single cell problem,
except for those provided in [29] on asymptotic solutions. For the one-cell
problem (1) with smooth interface, but $\sigma_{1}\neq\sigma_{0}$, the
solution is probably regular. Intuitively, the interface problem with
piecewise smooth coefficients and fixed transmembrane potential (that is, at
the equilibrium) can be recast, via lifting [7], to a more classical interface
problem already studied by babuska1970 [2, 23], who showed that the solution
is at least $H^{2}$ on each subdomain.

Well-posedness results for the general EMI problem are found in [3, 15],
however we are not aware of regularity results for the general EMI problem
(33). For two or more cells in contact, subdomains must be polygonal, which
limits the regularity. The singularities introduced by the contact have been
analyzed by nicaise1994 [31]. Moreover, the transmembrane voltage cannot be
smooth on $\Gamma_{0}$, because it has multiple branches. These results should
be taken in consideration in the development of higher order numerical schemes
for the solution of the EMI problem.

We also show that the discretization parameters are not too restrictive, when
compared to the more standard (homogenized) bidomain model. A typical time
step is $0.02\text{\,}\mathrm{ms}$ or lower for IMEX solvers [25]. In space,
we observed here that a mesh resolution for the membrane of
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is sufficiently accurate for the
cell-by-cell model. On the other hand, for the standard bidomain the mesh size
depends on the front thickness, in turn depending on tissue excitability and
conductivity. An accepted value is $100\text{\,}\mathrm{\SIUnitSymbolMicro
m}200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ [32]. However, this is only true
in the fiber direction and for healthy tissue, since in fibrotic tissue or in
the cross-fiber direction the front thickness is generally lower [33, 34].
Finally, we observed here that propagation failure can occur in the EMI model,
in contrast to the bidomain model. This aspect is very important in the study
of pathological situations.

This work paves the way in two directions. First, for designing another method
where a more realistic cell-by-cell model is solved, i.e., where the linear
algebraic condition on the gap junctions is replaced by a stiff nonlinear
ordinary differential equation. Second, it provides the mathematical framework
for solving the EMI model in three dimensions with a boundary integral
formulation. Once assembled, the solution of the problem is very fast, opening
interesting opportunities for long simulations.

## Acknowledgement

We are very grateful to Michael Multerer for kindly sharing the BEM code and
the fruitful discussion on the problem. We also thank the MICROCARD consortium
for the suggestions for improving the manuscript.

## References

  * [1] Assyr Abdulle, Marcus J. Grote, and Giacomo Rosilho de Souza.  Explicit stabilized multirate method for stiff differential equations.  Mathematics of Computation, 91:2681–2714, 2022.
  * [2] Ivo Babuška.  The finite element method for elliptic equations with discontinuous coefficients.  Computing, 5(3):207–213, 1970.
  * [3] Pierre-Elliott Becue.  Modélisation et simulation de l’électrophysiologie cardiaque à l’échelle microscopique. Theses, Université de Bordeaux, December 2018.
  * [4] Pierre-Elliott Bécue, Florian Caro, Mark Potse, and Yves Coudière.  Theoretical and Numerical Study of Cardiac Electrophysiology Problems at the Microscopic Scale.  SIAM Conference on the Life Sciences (LS16), July 2016.  Poster.
  * [5] Pierre-Elliott Bécue, Mark Potse, and Yves Coudière.  Microscopic simulation of the cardiac electrophysiology: A study of the influence of different gap junctions models.  In 2018 Computing in Cardiology Conference (CinC), volume 45, pages 1–4. IEEE, 2018.
  * [6] Pierre-Elliott Bécue, Mark Potse, and Yves Coudière.  A three-dimensional computational model of action potential propagation through a network of individual cells.  In 2017 Computing in Cardiology (CinC), pages 1–4, 2017.
  * [7] Zhiming Chen and Jun Zou.  Finite element methods and their convergence for elliptic and parabolic interface problems.  Numerische Mathematik, 79(2):175–202, 1998.
  * [8] 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:22–48, 2011.
  * [9] Michael Clerx, Pieter Collins, Enno de Lange, and Paul G.A. Volders.  Myokit: A simple interface to cardiac cellular electrophysiology.  Progress in Biophysics and Molecular Biology, 120:100–114, 1 2016\.
  * [10] Marc Courtemanche, Rafael J Ramirez, and Stanley Nattel.  Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model.  American Journal of Physiology-Heart and Circulatory Physiology, 275:H301–H321, 1998.
  * [11] Autumn A. Cuellar, Catherine M. Lloyd, Poul F. Nielsen, David P. Bullivant, David P. Nickerson, and Peter J. Hunter.  An overview of cellml 1.1, a biological model description language.  SIMULATION, 79:740–747, 12 2003.
  * [12] Cécile Daversin-Catty, Chris N Richardson, Ada J Ellingsrud, and Marie E Rognes.  Abstractions and automated algorithms for mixed domain finite element methods.  ACM Transactions on Mathematical Software (TOMS), 47(4):1–36, 2021\.
  * [13] 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:479–485, 1980.
  * [14] K. R. Foster and A. E. Sowers.  Dielectrophoretic forces and potentials induced on pairs of cells in an electric field.  Biophysical Journal, 69:777–784, 1995.
  * [15] Piero Colli Franzone and Giuseppe Savaré.  Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level.  Evolution Equations, Semigroups and Functional Analysis: in memory of Brunello Terreni, pages 49–78, 2002.
  * [16] Lia Gander, Rolf Krause, Michael Multerer, and Simone Pezzuto.  Space–time shape uncertainties in the forward and inverse problem of electrocardiography.  International Journal for Numerical Methods in Biomedical Engineering, 37:1–23, 2021.
  * [17] Gaël Guennebaud, Benoît Jacob, et al.  Eigen v3.  http://eigen.tuxfamily.org, 2010.
  * [18] Paul E Hand, Boyce E Griffith, and Charles S Peskin.  Deriving macroscopic myocardial conductivities by homogenization of microscopic models.  Bulletin of Mathematical Biology, 71:1707–1726, 2009.
  * [19] Fernando Henríquez and Carlos Jerez-Hanckes.  Multiple traces formulation and semi-implicit scheme for modelling biological cells under electrical stimulation.  ESAIM: Mathematical Modelling and Numerical Analysis, 52:659–702, 2018.
  * [20] Fernando Henríquez, Carlos Jerez-Hanckes, and Fernando Altermatt.  Boundary integral formulation and semi-implicit scheme coupling for modeling cells under electrical stimulation.  Numerische Mathematik, 136:101–145, 5 2017.
  * [21] Hervé Hogues, L. Joshua Leon, and Fernand A. Roberge.  A model study of electric field interactions between cardiac myocytes.  IEEE Transactions on Biomedical Engineering, 39:1232–1243, 1992\.
  * [22] Karoline Horgmo Jæger, Andrew G. Edwards, Wayne R. Giles, and Aslak Tveito.  Arrhythmogenic influence of mutations in a myocyte-based computational model of the pulmonary vein sleeve.  Scientific Reports, 12, 12 2022.
  * [23] Bruce R Kellogg.  On the poisson equation with intersecting interfaces.  Applicable Analysis, 4(2):101–129, 1974.
  * [24] André G Kléber and Yoram Rudy.  Basic mechanisms of cardiac impulse propagation and associated arrhythmias.  Physiological reviews, 2004.
  * [25] Dorian Krause, Mark Potse, Thomas Dickopf, Rolf Krause, Angelo Auricchio, and Frits W. Prinzen.  Hybrid parallelization of a large-scale heart model.  In Rainer Keller, David Kramer, and Jan-Philipp Weiss, editors, Facing the Multicore-Challenge II, volume 7174 of Lecture Notes in Computer Science, pages 120–132, Berlin, 2012. Springer.
  * [26] Rainer Kress.  Linear Integral Equations.  Springer, 1989.
  * [27] U. Langer and O. Steinbach.  Boundary element tearing and interconnecting methods.  Computing, 71:205–228, 2003.
  * [28] L. J. Leon and F. A. Roberge.  A model study of extracellular stimulation of cardiac cells.  IEEE Transactions on Biomedical Engineering, 40:1307–1319, 1993\.
  * [29] Hiroshi Matano and Yoichiro Mori.  Global existence and uniqueness of a three-dimensional model of cellular electrophysiology.  Discrete Contin. Dyn. Syst, 29(4):1573–1636, 2011.
  * [30] JC Neu and W Krassowska.  Homogenization of syncytial tissues.  Critical reviews in biomedical engineering, 21(2):137–199, 1993\.
  * [31] Serge Nicaise and Anna-Margarete Sändig.  General interface problems—i.  Mathematical Methods in the Applied Sciences, 17(6):395–429, 1994\.
  * [32] Steven A. Niederer, Eric Kerfoot, Alan P. Benson, Miguel O. Bernabeu, Olivier Bernus, Chris Bradley, Elizabeth M. Cherry, Richard Clayton, Flavio H. Fenton, Alan Garny, Elvio Heidenreich, Sander Land, Mary Maleckar, Pras Pathmanathan, Gernot Plank, José F. Rodríguez, Ishani Roy, Frank B. Sachse, Gunnar Seemann, Ola Skavhaug, and Nic 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.
  * [33] Simone Pezzuto, J. Hake, and Joakim Sundnes.  Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology.  International Journal for Numerical Methods in Biomedical Engineering, 32:e02762, 2016.  Referenze 12,13,14.
  * [34] Carlos Ruiz Herrera, Thomas Grandits, Gernot Plank, Paris Perdikaris, Francisco Sahli Costabal, and Simone Pezzuto.  Physics-informed neural networks to learn cardiac fiber orientation from multiple electroanatomical maps.  Engineering with Computers, 38(5):3957–3973, 2022.
  * [35] Stefan A. Sauter and Cristoph Schwab.  Boundary element methods.  Springer, 2011.
  * [36] Ulrich Schotten, Sander Verheule, Paulus Kirchhof, and Andreas Goette.  Pathophysiological mechanisms of atrial fibrillation: a translational appraisal.  Physiological Reviews, 91(1):265–325, 2011.
  * [37] Madison S Spach and J Francis Heidlage.  The stochastic nature of cardiac propagation at a microscopic level: electrical description of myocardial architecture and its application to conduction.  Circulation research, 76(3):366–380, 1995.
  * [38] Olaf Steinbach.  Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements.  Springer, 2007.
  * [39] Jeroen Stinstra, Rob MacLeod, and Craig Henriquez.  Incorporating histology into a 3d microscopic computer model of myocardium to study propagation at a cellular level.  Annals of Biomedical Engineering, 38:1399–1414, 4 2010.
  * [40] Jeroen G. Stinstra, Bruce Hopenfeld, and Rob S. MacLeod.  On the passive cardiac conductivity.  Annals of Biomedical Engineering, 33:1743–1751, 2005.
  * [41] JG Stinstra, CS Henriquez, and RS MacLeod.  Comparison of microscopic and bidomain models of anisotropic conduction.  Computers in Cardiology, 36:657–660, 2009.
  * [42] JG Stinstra, S Poelzing, RS MacLeod, and CS Henriquez.  A model for estimating the anisotropy of the conduction velocity in cardiac tissue based on the tissue morphology.  Computers in Cardiology, 34:129–132, 2007.
  * [43] Aslak Tveito, Karoline H. Jæger, Miroslav Kuchta, Kent Andre Mardal, and Marie E. Rognes.  A cell-based framework for numerical modeling of electrical conduction in cardiac tissue.  Frontiers in Physics, 5, 10 2017.
  * [44] Aslak Tveito, Kent-Andre Mardal, and Marie E Rognes.  Modeling Excitable Tissue: The EMI Framework.  Springer Nature, 2021.
  * [45] Edward J Vigmond, Jose L Perez Velazquez, Taufik A Valiante, Berj L Bardakjian, and Peter L Carlen.  Mechanisms of electrical coupling between pyramidal cells.  Journal of Neurophysiology, 78(6):3107–3116, 1997.

##  Appendix A Poincaré–Steklov map on $\Omega_{0}$

Here we briefly describe how to derive the Poincaré–Steklov operators
$\mathcal{P}_{0},P_{0}$ on $\Omega_{0}$ used in Eqs. 9 and 17, respectively.

Let ${\Phi}\in\\{\Gamma_{0},\Sigma\\}$, we introduce the restricted trace
operators

$\displaystyle\gamma^{\text{0}}_{t,{\Phi}}$ $\displaystyle\colon
H^{1}(\Omega_{0})\to H^{1/2}({\Phi}),\quad$
$\displaystyle\gamma^{\text{0}}_{t,{\Phi}}u_{0}(\bm{x})$
$\displaystyle=\lim_{\Omega_{0}\ni\bm{y}\to\bm{x}\in{\Phi}}u_{0}(\bm{y}),$
(67) $\displaystyle\gamma^{\text{0}}_{n,{\Phi}}$ $\displaystyle\colon
H^{1}(\Omega_{0})\to H^{-1/2}({\Phi}),$
$\displaystyle\gamma^{\text{0}}_{n,{\Phi}}u_{0}(\bm{x})$
$\displaystyle=\lim_{\Omega_{0}\ni\bm{y}\to\bm{x}\in{\Phi}}\langle\nabla
u_{0}(\bm{y}),\bm{n}_{0}\rangle.$

From the Green’s representation formula we have

$\displaystyle u_{0}(\bm{x})$
$\displaystyle=\int_{\Gamma_{0}}\gamma^{\text{0}}_{t,\bm{y}}G(\bm{x},\bm{y})\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}(\bm{y})\operatorname{d\\!}s_{\bm{y}}+\int_{\Sigma}\gamma^{\text{0}}_{t,\bm{y}}G(\bm{x},\bm{y})\gamma^{\text{0}}_{n,\Sigma}u_{0}(\bm{y})\operatorname{d\\!}s_{\bm{y}}$
(68)
$\displaystyle\quad-\int_{\Gamma_{0}}\gamma^{\text{0}}_{n,\bm{y}}G(\bm{x},\bm{y})\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}(\bm{y})\operatorname{d\\!}s_{\bm{y}}-\int_{\Sigma}\gamma^{\text{0}}_{n,\bm{y}}G(\bm{x},\bm{y})\gamma^{\text{0}}_{t,\Sigma}u_{0}(\bm{y})\operatorname{d\\!}s_{\bm{y}},$
$\displaystyle\bm{x}\in\Omega_{0}.$

For $\Phi,\Psi\in\\{\Gamma_{0},\Sigma\\}$ , we define

$\displaystyle\mathcal{V}_{0}^{\Psi,\Phi}$ $\displaystyle\colon
H^{-1/2}(\Phi)\to H^{1/2}(\Psi),\quad$
$\displaystyle\mathcal{V}_{0}^{\Psi,\Phi}\rho(\bm{x})$
$\displaystyle=\int_{\Phi}\gamma^{\text{0}}_{0,\bm{y}}G(\bm{x},\bm{y})\rho(\bm{y})\operatorname{d\\!}s_{\bm{y}},\quad$
$\displaystyle\bm{x}$ $\displaystyle\in\Psi,$ (69)
$\displaystyle\mathcal{K}_{0}^{\Psi,\Phi}$ $\displaystyle\colon
H^{1/2}(\Phi)\to H^{1/2}(\Psi),$
$\displaystyle\mathcal{K}_{0}^{\Psi,\Phi}\rho(\bm{x})$
$\displaystyle=\int_{\Phi}\gamma^{\text{0}}_{1,\bm{y}}G(\bm{x},\bm{y})\rho(\bm{y})\operatorname{d\\!}s_{\bm{y}},$
$\displaystyle\bm{x}$ $\displaystyle\in\Psi,$

applying the trace operators $\gamma^{\text{0}}_{t,\Gamma_{0}}$ and
$\gamma^{\text{0}}_{n,\Sigma}$ to Eq. 68 yields

$\displaystyle\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}$
$\displaystyle=\mathcal{V}_{0}^{\Gamma_{0},\Gamma_{0}}\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}+\mathcal{V}_{0}^{\Gamma_{0},\Sigma}\gamma^{\text{0}}_{n,\Sigma}u_{0}-(\mathcal{K}_{0}^{\Gamma_{0},\Gamma_{0}}-\tfrac{1}{2}I)\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}-\mathcal{K}_{0}^{\Gamma_{0},\Sigma}\gamma^{\text{0}}_{t,\Sigma}u_{0},$
(70) $\displaystyle\gamma^{\text{0}}_{t,\Sigma}u_{0}$
$\displaystyle=\mathcal{V}_{0}^{\Sigma,\Gamma_{0}}\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}+\mathcal{V}_{0}^{\Sigma,\Sigma}\gamma^{\text{0}}_{n,\Sigma}u_{0}-\mathcal{K}_{0}^{\Sigma,\Gamma_{0}}\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}-(\mathcal{K}_{0}^{\Sigma,\Sigma}-\tfrac{1}{2}I)\gamma^{\text{0}}_{t,\Sigma}u_{0},$

which, after manipulation and setting $\gamma^{\text{0}}_{n,\Sigma}u_{0}=0$
(cf. Eq. 33g), result in

$\displaystyle(\mathcal{K}_{0}^{\Gamma_{0},\Gamma_{0}}+\tfrac{1}{2}I)\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}$
$\displaystyle=\mathcal{V}_{0}^{\Gamma_{0},\Gamma_{0}}\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}-\mathcal{K}_{0}^{\Gamma_{0},\Sigma}\gamma^{\text{0}}_{t,\Sigma}u_{0},$
(71)
$\displaystyle(\mathcal{K}_{0}^{\Sigma,\Sigma}+\tfrac{1}{2}I)\gamma^{\text{0}}_{t,\Sigma}u_{0}$
$\displaystyle=\mathcal{V}_{0}^{\Sigma,\Gamma_{0}}\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}-\mathcal{K}_{0}^{\Sigma,\Gamma_{0}}\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}.$

Solving for $\gamma_{n,\Gamma_{0}}u_{0},\gamma_{t,\Sigma}u_{0}$ with respect
to $\gamma_{t,\Gamma_{0}}u_{0}$ yields the linear relation Eq. 9 (dropping
$\Gamma_{0}$ from the notation).

We discretize $\Gamma_{0}=\Gamma_{1}$ as in Section 2.2 (same collocation
points) and place $\overline{M}$ collocation points $\overline{\bm{x}}_{j}$ on
$\Sigma$. We compute a smooth parametrization
$\gamma_{\Sigma}\mathrel{\mathop{\ordinarycolon}}[0,1)\rightarrow\mathbb{R}^{d}$
of $\Sigma$, $\gamma_{\Sigma}(s_{j})=\overline{\bm{x}}_{j}$, by Fourier
interpolation (as for $\Gamma_{0}$) and represent

$\displaystyle\gamma^{\text{0}}_{t,\Gamma_{0}}u_{0}(\gamma_{\Gamma_{0}}(t))$
$\displaystyle=\sum_{j=1}^{M}u_{0,\Gamma_{0}}^{j}L_{j}(t),$
$\displaystyle\gamma^{\text{0}}_{n,\Gamma_{0}}u_{0}(\gamma_{\Gamma_{0}}(t))$
$\displaystyle=\sum_{j=1}^{M}\tilde{u}_{0,\Gamma_{0}}^{j}L_{j}(t),$ (72)
$\displaystyle\gamma^{\text{0}}_{t,\Sigma}u_{0}(\gamma_{\Sigma}(t))$
$\displaystyle=\sum_{j=1}^{\overline{M}}u_{0,\Sigma}^{j}\overline{L}_{j}(t),$

with $\overline{L}_{j}(s)$ the trigonometric Lagrange basis functions
satisfying $\overline{L}_{j}(s_{i})=\delta_{ij}$, $i,j=1,\ldots,\overline{M}$.
Inserting Eq. 72 into Eq. 71 yields

$\displaystyle(K_{0}^{\Gamma_{0},\Gamma_{0}}+\tfrac{1}{2}I)\bm{u}_{0,\Gamma_{0}}$
$\displaystyle=V_{0}^{\Gamma_{0},\Gamma_{0}}\bm{\tilde{u}}_{0,\Gamma_{0}}-K_{0}^{\Gamma_{0},\Sigma}\bm{u}_{0,\Sigma},$
(73)
$\displaystyle(K_{0}^{\Sigma,\Sigma}+\tfrac{1}{2}\overline{I})\bm{u}_{0,\Sigma}$
$\displaystyle=V_{0}^{\Sigma,\Gamma_{0}}\bm{\tilde{u}}_{0,\Gamma_{0}}-K_{0}^{\Sigma,\Gamma_{0}}\bm{u}_{0,\Gamma_{0}},$

with
$(\bm{u}_{0,\Gamma_{0}})_{j}=u_{0,\Gamma_{0}}^{j},(\bm{\tilde{u}}_{0,\Gamma_{0}})_{j}=\tilde{u}_{0,\Gamma_{0}}^{j},(\bm{u}_{0,\Sigma})_{j}=u_{0,\Sigma}^{j}$
and

$\displaystyle(K_{0}^{\Gamma_{0},\Gamma_{0}})_{kj}$
$\displaystyle=\mathcal{K}_{0}^{\Gamma_{0},\Gamma_{0}}(L_{j}\circ\gamma_{\Gamma_{0}}^{-1})(\bm{x}_{k}),$
$\displaystyle(K_{0}^{\Sigma,\Gamma_{0}})_{kj}$
$\displaystyle=\mathcal{K}_{0}^{\Sigma,\Gamma_{0}}(L_{j}\circ\gamma_{\Gamma_{0}}^{-1})(\overline{\bm{x}}_{k}),$
(74) $\displaystyle(K_{0}^{\Sigma,\Sigma})_{kj}$
$\displaystyle=\mathcal{K}_{0}^{\Sigma,\Sigma}(\overline{L}_{j}\circ\gamma_{\Sigma}^{-1})(\overline{\bm{x}}_{k}),$
$\displaystyle(K_{0}^{\Gamma_{0},\Sigma})_{kj}$
$\displaystyle=\mathcal{K}_{0}^{\Gamma_{0},\Sigma}(\overline{L}_{j}\circ\gamma_{\Sigma}^{-1})(\bm{x}_{k}),$
(75)

and similarly for $V_{0}^{\Sigma,\Gamma_{0}},V_{0}^{\Gamma_{0},\Gamma_{0}}$.
Solving for $\bm{\tilde{u}}_{0,\Gamma_{0}},\bm{u}_{0,\Sigma}$ with respect to
$\bm{u}_{0,\Gamma_{0}}$ yields

$\bm{\tilde{u}}_{0,\Gamma_{0}}=P_{0}\bm{u}_{0,\Gamma_{0}},$ (76)

which is employed in Eqs. 17 and 18 (dropping $\Gamma_{0}$ from the notation).