# The Two-Point Finite Volume Scheme for the Microscopic Bidomain Model of Electrocardiology

Zeina Chehade, Yves Coudière

<figure><img src="image_1.png" /></figure>

► To cite this version:

Zeina Chehade, Yves Coudière. The Two-Point Finite Volume Scheme for the Microscopic Bidomain Model of Electrocardiology. FVCA 10 - Finite Volume for Complex applications 10, Oct 2023, Strasbourg, France. hal-04403881

HAL Id: hal-04403881

https://inria.hal.science/hal-04403881v1

Submitted on 18 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The Two-Point Finite Volume Scheme for the Microscopic Bidomain Model of Electrocardiology

Zeina Chehade¹ and Yves Coudière¹

Univ. Bordeaux, CNRS, Inria, Bordeaux INP, IMB, UMR 5251, IHU Liryc, F-33400 Talence, France,
zeina.chehade@u-bordeaux.fr, yves.coudiere@u-bordeaux.fr

**Abstract.** We are interested in the approximation of a cardiac microscopic model. It is a set of Laplace equations with non-standard, time-dependent, transmission conditions, for which finite volume methods are of real interest. The transmission conditions state as ordinary differential equations on the jump of the potential, namely the transmembrane voltage, so that we keep this voltage as an unknown in our scheme. Here we extend the two-point flux approximation to the discretization of this model, show that it converges, and compute error estimates.

**Keywords:** finite volumes, error estimate, cardiac EMI model

## 1 Introduction

The standard mathematical model for describing the spread of excitation of a cardiac tissue is a system of reaction-diffusion equations, obtained by homogenization of a microscopic model. It fails to capture effects of small scale tissue organization, relevant to cardiac fibrillation [3]. For these reasons, we aim at finding numerical approximations of the microscopic cardiac model, called the extracellular-membrane-intracellular (EMI) model, as defined in [1, 5].

The general EMI model is a set of $N_C + 1$ Laplace equations. For sake of simplicity, we will work only with $N_C = 1$ cell, denoted $\Omega_1$ within a tissue sample $\Omega (\bar{\Omega}_1 \subset \Omega)$. The extracellular matrix (ECM) is denoted by $\Omega_0 := \Omega \setminus \bar{\Omega}_1$ (see Fig. 1). Scalar electrical conductivity coefficients $\sigma_0 > 0$ and $\sigma_1 > 0$ are given in the two subdomains. The main physical unknowns are the electrical potentials $u_0$ and $u_1$. The normal flux of current at the interface between subdomains is continuous, whereas the electrical potential jumps between subdomains. This jump defines the voltage $v = u_1 - u_0$ on the cell membrane $\Sigma := \bar{\Omega}_0 \cap \bar{\Omega}_1$. The system is closed by ordinary differential equations (ODEs) that link the flux of current to the voltage on the cell membrane. The equations read

$$ -\sigma_k \Delta u_k = 0, \quad \text{in } \Omega_k, \quad k = 0, 1, \quad (1) $$

$$ -\sigma_0 \nabla u_0 \cdot n_0 = \sigma_1 \nabla u_1 \cdot n_1 = -(c_m \partial_t v + f(v)), \quad \text{on } \Sigma, \quad (2) $$

where the unit vectors $n_0, n_1$ are normal to $\Sigma$ outward of $\Omega_0$ and $\Omega_1$, as depicted on Fig. 1. The function $f$ defines the ionic current, and models membrane
electrophysiology, while $c_m > 0$ is the membrane surface capacitance. In general, the function $f$ also depends on a set of state variables $w$ defined on $\Sigma$ (i.e. $f = f(v, w)$), and the model is closed by a system of ODE, $\partial_t w = g(v, w)$ on $\Sigma$. Here we consider a simplified model without state variables $w$. The equations are completed by mixed boundary conditions $-\sigma_0 \nabla u_0 \cdot n_0 = g^N$ on $\Gamma^N$, and $u_0 = g^D$ on $\Gamma^D$, where $\Gamma^D \cup \Gamma^N = \partial\Omega$ is the boundary of the tissue sample. The equations are supplemented with initial data on $\Sigma$, $v(0, \cdot) = v^0$, and we look for a solution $u_0, u_1$ for $t \in [0, T]$ for any $T > 0$. From previous work [1], we know that there exists a weak solution in $L^2(0, T; H^1(\Omega_0) \times H^1(\Omega_1))$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Geometrical setup of the EMI model with one cell</figcaption></figure>

There has been some attempts at solving this system of equations by finite elements or boundary element methods [1, 4, 5], but not with finite volumes methods (FVM). Since the dynamics of the phenomena is written on the interface between subdomains (through jump conditions), and the model relies on the continuity of the flux of current, we believe that FVM are a relevant choice for the EMI model. We start here by studying the simple two-point flux approximation (TPFA) for Eqs (1) and (2). It is essential to carefully write the discrete flux on the cell membrane $\Sigma$. Following the usual track, and assuming enough regularity (see section 4), we are able to compute some error estimates for the jump $v$ in $L^\infty(0, T; L^2(\Sigma))$, as well as in a discrete $L^2(0, T, H^1)$-like norm for $u_0$ and $u_1$. For this purpose, we assume that $\Omega_0$ and $\Omega_1$ are polygonal domains.

Section 2, 3, and 4 describe the scheme, its well-posedness, and the error estimates; section 5 discusses the current result, and on-going work.

## 2 The Numerical Scheme

We consider a standard FV-admissible mesh $\mathcal{T}$, as defined in e.g. [2], with control volumes $K \in \mathcal{T}$, and cell centers $(x_K)_{K \in \mathcal{T}}$. We assume that the mesh $\mathcal{T}$ is consistent with the subdomains $\Omega_0$ and $\Omega_1$: any $K \in \mathcal{T}$ is such that, either $K \subset \Omega_0$, or $K \subset \Omega_1$, so that we define $\mathcal{T}_i = \{K \in \mathcal{T}, K \subset \Omega_i\}$ ($i = 0, 1$). We also split the set of the mesh interfaces into: the set $\mathcal{E}^*$ of interfaces $e = K|L$ (using notations from [2]) such that $(K, L) \in \mathcal{T}_i \times \mathcal{T}_i$ for some $i \in \{0, 1\}$; the set $\mathcal{E}^\Sigma$ of interfaces $e = K|L$ such that $(K, L) \in \mathcal{T}_0 \times \mathcal{T}_1$; and the set $\mathcal{E}^D$ (resp. $\mathcal{E}^N$) of the boundary faces $e = K$ for which $K \subset \Omega_0$ and $e \in \Gamma^D$ (resp. $\Gamma^N$). Let
$\mathcal{E} = \mathcal{E}^* \cup \mathcal{E}^\Sigma \cup \mathcal{E}^D \cup \mathcal{E}^N$, and denote by $(x_e)_{e \in \mathcal{E}}$ the points $x_e$ at the intersection of the interface $e$ and the line $(x_K, x_L)$ if $e \in \mathcal{E}^* \cup \mathcal{E}^\Sigma$, and the perpendicular projection of $x_K$ on $e$ if $e \in \mathcal{E}^D \cup \mathcal{E}^N$. Regarding the time discretization, we set $\Delta t = \frac{T}{N}$, for any $N > 0$ and $t^n = n\Delta t$, for $n = 0, \dots, N$.

An edge $e = K|L \in \mathcal{E}^\Sigma$ may be such that $(K, L) \in \mathcal{T}_0 \times \mathcal{T}_1$ or $(K, L) \in \mathcal{T}_1 \times \mathcal{T}_0$. Below, we always use the convention that $K \in \mathcal{T}_0$ and $L \in \mathcal{T}_1$.

The integral form of Eq. (1) on any cell $K \in \mathcal{T}$ at time $t^n$ reads

$$ - \sum_{e \in \mathcal{E}_K} \int_e \sigma_k \nabla u_k(t^n, \cdot) \cdot n_{Ke} = 0, \quad (3) $$

where $\mathcal{E}_K$ is the set of interfaces $e \in \mathcal{E}$ that form the boundary of $K \in \mathcal{T}$, and the vector $n_{Ke}$ is the unit normal to $e$ outward of $K$. The integral of Eq. (2) on any interface $e = K|L \in \mathcal{E}^\Sigma$ (with the above convention) at time $t^n$ reads

$$ \begin{aligned} - \int_e \sigma_0 \nabla u_0(t^n, \cdot) \cdot n_{Ke} &= \int_e \sigma_1 \nabla u_1(t^n, \cdot) \cdot n_{Le} \\ &= - \left( c_m \int_e \partial_t v(t^n, \cdot) + \int_e f(v(t^n, \cdot)) \right). \end{aligned} \quad (4) $$

There are $N_\mathcal{T} = \text{card}(\mathcal{T})$ equations (3), and $N_\Sigma = \text{card}(\mathcal{E}^\Sigma)$ equations (4), so that we take $N_\mathcal{T} + N_\Sigma$ unknowns at each time $t^n$ ($n = 1 \dots N$), which are the vectors $u_\mathcal{T}^n := (u_K^n)_{K \in \mathcal{T}} \in \mathbb{R}^{N_\mathcal{T}}$, approximating the $u(t^n, x_K)$, and $v_\mathcal{E}^n := (v_e^n)_{e \in \mathcal{E}^\Sigma} \in \mathbb{R}^{N_\Sigma}$, approximating the $v(t^n, x_e)$. Using a semi-implicit Euler time-stepping method, our numerical scheme states, for $n = 1, 2 \dots N$, as the $N_\mathcal{T}+N_\Sigma$ equations

$$ - \sum_{e \in \mathcal{E}_K} F_{Ke}^n = 0, \quad \text{for all } K \in \mathcal{T}, \quad (5) $$

$$ -F_{Ke}^n = F_{Le}^n = - \left( \frac{c_m}{\Delta t} (v_e^n - v_e^{n-1}) + f(v_e^{n-1}) \right) |e|, \quad \text{for all } e = K|L \in \mathcal{E}^\Sigma, \quad (6) $$

with some given initial values $v_\mathcal{E}^0 = (v_e^0)_{e \in \mathcal{E}^\Sigma} \in \mathbb{R}^{N_\Sigma}$. The numerical flux formula $F_{Ke}^n$ for $e \in \mathcal{E}_K$ and $K \in \mathcal{T}$ approximate the exact flux:

$$ F_{Ke}^n = \tau_e (u_L^n - u_K^n) \quad \text{if } e = K|L \in \mathcal{E}_K \cap \mathcal{E}^*, \quad K \in \mathcal{T}, \quad (7) $$

$$ F_{Ke}^n = \tau_e (u_L^n - u_K^n - v_e^n) \quad \text{if } e = K|L \in \mathcal{E}_K \cap \mathcal{E}^\Sigma, \quad K \in \mathcal{T}_0, L \in \mathcal{T}_1, \quad (8) $$

$$ F_{Ke}^n = \tau_e (g_e^{D,n} - u_K^n) \quad \text{if } e = K| \in \mathcal{E}_K \cap \mathcal{E}^D, \quad K \in \mathcal{T}_0, \quad (9) $$

$$ F_{Ke}^n = -g_e^{N,n} |e| \quad \text{if } e = K| \in \mathcal{E}_K \cap \mathcal{E}^N, \quad K \in \mathcal{T}_0. \quad (10) $$

In these expressions, the coefficient $\tau_e$ has the usual value $\tau_e = \frac{\tau_{Ke} \tau_{Le}}{\tau_{Ke} + \tau_{Le}}$ if $e = K|L \in \mathcal{E}^* \cup \mathcal{E}^\Sigma$, and $\tau_e = \tau_{Ke}$ if $e = K| \in \mathcal{E}^D$, where $\tau_{Ke} = \frac{\sigma_i |e|}{d_{Ke}}$ for all $e \in \mathcal{E}_K$ and $K \in \mathcal{T}_i$ ($i = 0, 1$). The quantities $|e|$ and $d_{Ke}$ are the measure of the interface $e$ and the Euclidean distance $d_{Ke} = d(x_K, e)$. At last, we take $g_e^{D,n} = g^D(t^n, x_e)$ and $g_e^{N,n} |e| = \int_e g^N(t^n, \cdot)$.
Expressions (7), (9), and (10) are standard, but expression (8) is obtained on an interface $e = K|L$ between $K \in \mathcal{T}_0$ and $L \in \mathcal{T}_1$ after introducing two auxiliary unknowns, named $u_{e,K}$ and $u_{e,L}$, approximating $u_0(t^n, x_e)$ and $u_1(t^n, x_e)$, and flux expressions $F_{Ke}^n = \tau_{Ke}(u_{K,e} - u_K)$ and $F_{Le}^n = \tau_{Le}(u_{L,e} - u_L)$. The auxiliary unknowns are eliminated with the conservation and jump conditions:

$$ \tau_{Ke}(u_{K,e} - u_K^n) + \tau_{Le}(u_{L,e} - u_L^n) = 0, \quad u_{L,e} - u_{K,e} = v_e^n. \quad (11) $$

Instead of the usual solution, we find that $u_{K,e} = \frac{\tau_{Ke}u_K^n + \tau_{Le}u_L^n - \tau_{Le}v_e^n}{\tau_{Ke} + \tau_{Le}}$, and $u_{L,e} = \frac{\tau_{Ke}u_K^n + \tau_{Le}u_L^n + \tau_{Ke}v_e^n}{\tau_{Ke} + \tau_{Le}}$, and then Formula (8) for the flux.

### 3 Discrete norms, discrete problem, coercivity

In this section, we establish existence and uniqueness of the discrete solution.

First, from Eqs. (5) and (6), and the definitions (7) - (10), we observe that the scheme writes as the following block linear system, at each time step

$$ \begin{pmatrix} A & B \\ B^\top & C + \frac{c_m}{\Delta t} D \end{pmatrix} \begin{pmatrix} u_{\mathcal{T}}^n \\ v_{\mathcal{E}\Sigma}^n \end{pmatrix} = \left( - \left( -\frac{c_m}{\Delta t} v_{\mathcal{E}\Sigma}^{n-1} + f(v_{\mathcal{E}\Sigma}^{n-1}) \right) \middle|_{e} \right), \quad (12) $$

where $A \in \mathbb{R}^{N_{\mathcal{T}} \times N_{\mathcal{T}}}$ is the usual TPFA matrix, with nonzero entries $a_{KL} = -\tau_e$ if $e = K|L$, and $a_{KK} = \sum_{e=K|L \in \mathcal{E}_K} \tau_e$. The matrix $B \in \mathbb{R}^{N_{\mathcal{T}} \times N_{\mathcal{E}\Sigma}}$ has nonzero entries $b_{Ke} = \tau_e$, and $b_{Le} = -\tau_e$ (for $e \in \mathcal{E}^{\Sigma}$), and the matrices $C, D \in \mathbb{R}^{N_{\mathcal{E}\Sigma} \times N_{\mathcal{E}\Sigma}}$ are diagonal, with $c_{ee} = \tau_e$ and $d_{ee} = |e|$ (for $e \in \mathcal{E}^{\Sigma}$). The vector $G^n \in \mathbb{R}^{N_{\mathcal{T}}}$ gathers the contributions of the boundary data.

We multiply scalarly Eq. (12) by the unknown vector $(u_{\mathcal{T}}^n, v_{\mathcal{E}\Sigma}^n)^\top$, or equivalently, we multiply Eq. (5) by $u_K^n$ and sum over $K \in \mathcal{T}$, multiply Eq. (6) by $v_e^n$ and sum over $e \in \mathcal{E}^{\Sigma}$. After reordering the summation over the set of edges, it yields

$$ \begin{aligned} & \sum_{e \in \mathcal{E}^*} \tau_e |u_L^n - u_K^n|^2 + \sum_{e \in \mathcal{E}^{\Sigma}} \tau_e |u_L^n - u_K^n - v_e^n|^2 + \sum_{e \in \mathcal{E}^D} \tau_e |u_K^n|^2 + \sum_{e \in \mathcal{E}^{\Sigma}} \frac{c_m}{\Delta t} |v_e^n|^2 |e| \\ & = \sum_{e \in \mathcal{E}^D} \tau_e g_e^{D,n} u_K^n - \sum_{e \in \mathcal{E}^N} g_e^{N,n} u_K^n |e| + \sum_{e \in \mathcal{E}^{\Sigma}} \frac{c_m}{\Delta t} (v_e^{n-1} - f(v_e^{n-1})) v_e^n |e|. \end{aligned} \quad (13) $$

Equation (13) shows that the linear system (12) is symmetric and positive-definite. In addition, we can introduce the semi-norm of a vector $(u_{\mathcal{T}}, v_{\mathcal{E}\Sigma}) = ((u_K)_{K \in \mathcal{T}}, (v_e)_{e \in \mathcal{E}^{\Sigma}})$ by

$$ \begin{aligned} & |(u_{\mathcal{T}}, v_{\mathcal{E}\Sigma})|_{1,\mathcal{T}} \\ & := \left( \sum_{e \in \mathcal{E}^*} \tau_e |u_L - u_K|^2 + \sum_{e \in \mathcal{E}^{\Sigma}} \tau_e |u_L - u_K - v_e|^2 + \sum_{e \in \mathcal{E}^D} \tau_e |u_K|^2 \right)^{1/2}, \end{aligned} \quad (14) $$

and the norm of $v_{\mathcal{E}\Sigma} = (v_e)_{e \in \mathcal{E}^{\Sigma}} \in \mathbb{R}^{N_{\mathcal{E}\Sigma}}$ by $\|v_{\mathcal{E}\Sigma}\|_{0,\mathcal{E}\Sigma}^2 = \sum_{e \in \mathcal{E}^{\Sigma}} |v_e|^2 |e|$. This latter formula defines a norm, because, if $|(u_{\mathcal{T}}, v_{\mathcal{E}\Sigma})|_{1,\mathcal{T}} = 0$, then one observes
that $u_K = 0$ for all $K \in \mathcal{T}_0$, and there exists $u \in \mathbb{R}$ such that $u_L = u$ for all $L \in \mathcal{T}_1$, and at last $v_e = u$ on all $e \in \mathcal{E}^{\Sigma}$, so that if $|(u_{\mathcal{T}}, v_{\mathcal{E}\Sigma})|_{1,\mathcal{T}} = 0$, then there exists $u \in \mathbb{R}$ such that $u_K = v_e = u$ for all $K \in \mathcal{T}$ and $e \in \mathcal{E}^{\Sigma}$.

Using these norms, Eq. (13) rewrites as:

$$|(u_{\mathcal{T}}^n, v_{\mathcal{E}\Sigma}^n)|_{1,\mathcal{T}}^2 + \frac{c_m}{\Delta t} \|v_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 = (u_{\mathcal{T}}^{n\top} v_{\mathcal{E}\Sigma}^{n\top}) \left( \frac{A}{B^T} + \frac{c_m}{\Delta t} D \right) \left( \begin{pmatrix} u_{\mathcal{T}}^n \\ v_{\mathcal{E}\Sigma}^n \end{pmatrix} \right). \quad (15)$$

Finally, the FV formulation allows to write the problem as an evolution problem on $\Sigma$ only, computing the Schur complement of (12) and iterating with

$$(C - B^{\top} A^{-1} B + \frac{c_m}{\Delta t} D) v_{\mathcal{E}\Sigma}^n = \left( \frac{c_m}{\Delta t} v_{\mathcal{E}\Sigma}^{n-1} - f(v_{\mathcal{E}\Sigma}^{n-1}) \right) |e| - B^{\top} A^{-1} G^n. \quad (16)$$

## 4 Convergence analysis

In this section, a convergence analysis is carried out by proving the following error estimates theorem.

**Theorem 1.** Assume that $f$ is Lipschitz continuous in $\mathbb{R}$. Given $T > 0$, and the discrete setup from above, assume that $u_k \in C^2([0, T] \times \bar{\Omega}_k)$, and consider $\bar{u}_{\mathcal{T}}^n$ and $\bar{v}_{\mathcal{E}\Sigma}^n$ defined by $\bar{u}_K^n = u_i(t^n, x_K)$ for all $K \in \mathcal{T}_i$ ($i = 0, 1$), and $\bar{v}_e^n = u_1(t^n, x_e) - u_0(t^n, x_e)$ for all $e \in \mathcal{E}^{\Sigma}$, and $n = 0...N$. The discrete errors $(\epsilon_{\mathcal{T}}^n, \eta_{\mathcal{E}\Sigma}^n)$ are defined by $\epsilon_{\mathcal{T}}^n = \bar{u}_{\mathcal{T}}^n - u_{\mathcal{T}}^n$ and $\eta_{\mathcal{E}\Sigma}^n = \bar{v}_{\mathcal{E}\Sigma}^n - v_{\mathcal{E}\Sigma}^n$. We associate to $(\eta_{\mathcal{E}\Sigma}^n)_n$ the $L^\infty(0, T; L^2(\Sigma))$ function $\eta_{\mathcal{E}\Sigma}(t, x) = \eta_e^n$ for $t \in ]t^{n-1}, t^n[$ ($n = 1...N$), and $x \in e$ ($e \in \mathcal{E}^{\Sigma}$). If the initial approximation $v_{\mathcal{E}\Sigma}^0$ is such that $\|v^0 - v_{\mathcal{E}\Sigma}^0\|_{0,\Sigma} \le Ch$, then there exists $C > 0$, depending only on the data, such that

$$\|\eta_{\mathcal{E}\Sigma}\|_{L^{\infty}(0,T;L^2(\Sigma))} \le C(h + \Delta t), \quad \left( \sum_{n=1}^{N} \Delta t |(\epsilon_{\mathcal{T}}^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\mathcal{T}} \right)^{1/2} \le C(h + \Delta t), \quad (17)$$

where we have defined $h = \max_{K \in \mathcal{T}} \text{diam}(K)$.

*Proof.* We denote by $\hat{F}_{Ke}^n := \int_e \sigma_i \nabla u_i \cdot n_{Ke}$ the exact flux out of cell $K \in \mathcal{T}_i$ ($i = 0, 1$) through the edge $e \in \mathcal{E}_K$, and by $\bar{F}_{Ke}^n$ the flux associated to $(\bar{u}_{\mathcal{T}}^n, \bar{v}_{\mathcal{E}\Sigma}^n)$ with Formulas (7) - (10). By subtracting Eqs. (5) and (6) defining the approximation solution, from Eqs (3) and (4) verified by the exact solution, the equations on the error vectors $(\epsilon_{\mathcal{T}}^n, \eta_{\mathcal{E}\Sigma}^n)$ for $n = 1...N$, write

$$- \sum_{e \in \mathcal{E}_K} (\bar{F}_{Ke}^n - F_{Ke}^n) = - \sum_{e \in \mathcal{E}_K} R_{Ke}^n |e|, \quad (18)$$

on all $K \in \mathcal{T}$, where the consistency error is $R_{Ke}^n = \frac{1}{|e|} (\bar{F}_{Ke}^n - \hat{F}_{Ke}^n)$; and

$$-(\bar{F}_{Ke}^n - F_{Ke}^n) + |e| R_{Ke}^n = (\bar{F}_{Le}^n - F_{Le}^n) - |e| R_{Le}^n = -(c_m T_e^n + S_e^n) |e| \\ - \left( \frac{c_m}{\Delta t} (\eta_e^n - \eta_e^{n-1}) + f(\bar{v}_e^{n-1}) - f(v_e^{n-1}) \right) |e|, \quad (19)$$
on all $e = K|L \in \mathcal{E}^{\Sigma}$, where the additional consistency errors are defined by

$$T_e^n = \frac{1}{|e|} \int_e \partial_t v(t^n, \cdot) - \frac{\bar{v}_e^n - \bar{v}_e^{n-1}}{\Delta t}, \quad S_e^n = \frac{1}{|e|} \int_e f(v(t^n, \cdot)) - f(\bar{v}_e^{n-1}). \quad (20)$$

For $n = 0$, we have $\eta_{\mathcal{E}\Sigma}^0 = \bar{v}_{\mathcal{E}\Sigma}^0 - v_{\mathcal{E}\Sigma}^0$. We multiply by $\epsilon_T^n$ and $\eta_{\mathcal{E}\Sigma}^n$ and obtain (like we obtained (13) and the coercivity inequality (15))

$$|(\epsilon_T^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\tau}^2 + \frac{c_m}{\Delta t} \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 \le \left(\frac{c_m}{\Delta t} \eta_{\mathcal{E}\Sigma}^{n-1} - (f(\bar{v}_{\mathcal{E}\Sigma}^{n-1}) - f(v_{\mathcal{E}\Sigma}^{n-1})), \eta_{\mathcal{E}\Sigma}^n\right)_{0,\Sigma} \\ + \sum_{e=K|L \in \mathcal{E}^*} |e| R_{Ke}^n (\epsilon_L^n - \epsilon_K^n) + \sum_{e=K|L \in \mathcal{E}^D} |e| R_{Ke}^n (-\epsilon_K^n) + \sum_{e=K|L \in \mathcal{E}^N} |e| R_{Ke}^n (-\epsilon_K^n) \\ + \sum_{e=K|L \in \mathcal{E}^\Sigma} |e| R_{Ke}^n (\epsilon_L^n - \epsilon_K^n - \eta_e^n) - (c_m T_{\mathcal{E}\Sigma}^n + S_{\mathcal{E}\Sigma}^n, \eta_{\mathcal{E}\Sigma}^n)_{0,\Sigma},$$

where $(\cdot, \cdot)_{0,\Sigma}$ denotes the natural scalar product on $\mathcal{E}^{\Sigma}$ associated to $\|\cdot\|_{0,\Sigma}$, and $T_{\mathcal{E}\Sigma}^n := (T_e^n)_{e \in \mathcal{E}\Sigma} \in \mathbb{R}^{N\Sigma}$ and $S_{\mathcal{E}\Sigma}^n := (S_e^n)_{e \in \mathcal{E}\Sigma} \in \mathbb{R}^{N\Sigma}$. We have $R_{Ke}^n = 0$ for $e \in \mathcal{E}^N$ (see def of $g_e^{N,n}$ above). We have $|f(\bar{v}_e^{n-1}) - f(v_e^{n-1})| \le \lambda \eta_e^{n-1}$ for all $e \in \mathcal{E}^{\Sigma}$, where $\lambda > 0$ is the Lipschitz constant for $f$, and then we obtain

$$\Delta t |(\epsilon_T^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\tau}^2 + c_m \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 \le c_m \left(1 + \frac{\lambda}{c_m} \Delta t\right) \|\eta_{\mathcal{E}\Sigma}^{n-1}\|_{0,\Sigma} \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} \\ + \Delta t R^n |(\epsilon_T^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\tau} + \Delta t (c_m \|T_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} + \|S_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}) \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}, \quad (21)$$

where $R^n = R^n(u_0, u_1) = \left(\sum_{e \in \mathcal{E}^* \cup \mathcal{E}^D \cup \mathcal{E}^\Sigma} \frac{|e|^2}{\tau_e} |R_{Ke}^n|^2\right)^{1/2}$. Using Young inequalities, we first prove that

$$\frac{\Delta t}{2} |(\epsilon_T^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\tau}^2 + \frac{c_m}{2} \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 \le \frac{c_m}{2} \left(1 + \frac{\lambda}{c_m} \Delta t\right)^2 \|\eta_{\mathcal{E}\Sigma}^{n-1}\|_{0,\Sigma}^2 \\ + \frac{\Delta t}{2} |R^n|^2 + c_m \Delta t \left( \|T_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} + \frac{1}{c_m} \|S_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} \right) \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}. \quad (22)$$

The last term is bounded for any $\alpha > 0$ as follows:

$$c_m \Delta t \left( \|T_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} + \frac{1}{c_m} \|S_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} \right) \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} \le \frac{c_m \alpha}{2} \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 \\ + \frac{c_m}{2\alpha} \Delta t^2 \left( \|T_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} + \frac{1}{c_m} \|S_{\mathcal{E}\Sigma}^n\|_{0,\Sigma} \right)^2. \quad (23)$$

We combine (22) and (23), choose $\alpha > 0$ such that $1 - \alpha = \frac{1}{1 + \frac{\lambda}{c_m} \Delta t}$, to obtain

$$\Delta t |(\epsilon_T^n, \eta_{\mathcal{E}\Sigma}^n)|_{1,\tau}^2 + \frac{c_m}{1 + \frac{\lambda}{c_m} \Delta t} \|\eta_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 \le c_m \left(1 + \frac{\lambda}{c_m} \Delta t\right)^2 \|\eta_{\mathcal{E}\Sigma}^{n-1}\|_{0,\Sigma}^2 \\ + \Delta t |R^n|^2 + \left(1 + \frac{\lambda}{c_m} \Delta t\right) \frac{2\Delta t}{\lambda} (c_m^2 \|T_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2 + \|S_{\mathcal{E}\Sigma}^n\|_{0,\Sigma}^2). \quad (24)$$
From (24) we extract an estimation on $\eta_{\varepsilon\Sigma}^n$, and a recurrence shows that

$$ \|\eta_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2 \le \left(1 + \frac{\lambda}{c_m}\Delta t\right)^{3n} \|\eta_{\varepsilon\Sigma}^0\|_{0,\Sigma}^2 + \frac{\Delta t}{c_m} \sum_{i=0}^{n-1} \left(1 + \frac{\lambda}{c_m}\Delta t\right)^{3i+1} |R^{n-i}|^2 \\ + \frac{2\Delta t}{\lambda c_m} \sum_{i=0}^{n-1} \left(1 + \frac{\lambda}{c_m}\Delta t\right)^{3i+2} (c_m^2 \|T_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2 + \|S_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2), $$

which proves that, for all $n = 1...N$,

$$ \|\eta_{\varepsilon\Sigma}^n\|_{0,\Sigma} \le \exp\left(\frac{3}{2}\frac{\lambda}{c_m}T\right) \left( \|\eta_{\varepsilon\Sigma}^0\|_{0,\Sigma}^2 + \frac{1}{c_m}R_{\mathcal{T}}^2 + \frac{2c_m}{\lambda}T_{\varepsilon\Sigma}^2 + \frac{2}{\lambda c_m}S_{\varepsilon\Sigma}^2 \right)^{1/2}, $$

with $R_{\mathcal{T}}^2 = \Delta t \sum_{n=1}^N |R^n|^2$, $S_{\varepsilon\Sigma}^2 = \Delta t \sum_{n=1}^N \|S_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2$, $T_{\varepsilon\Sigma}^2 = \Delta t \sum_{n=1}^N \|T_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2$.

In a second step, we start again the computation from inequality (21), and show, after some standard calculations, that

$$ \Delta t \sum_{n=1}^{N} |(\epsilon_{\mathcal{T}}^n, \eta_{\varepsilon\Sigma}^n)|_{1,\mathcal{T}}^2 \le c_m \|\eta_{\varepsilon\Sigma}^0\|_{0,\Sigma}^2 + 2c_m^2 T_{\varepsilon\Sigma}^2 + 2S_{\varepsilon\Sigma}^2 + 2(\lambda+1)T \max_{n=1...N} \|\eta_{\varepsilon\Sigma}^n\|_{0,\Sigma}^2. $$

It remains to estimate the consistency errors $R_{\mathcal{T}}$, $S_{\varepsilon\Sigma}$, and $T_{\varepsilon\Sigma}$. Assuming $C^2$ regularity, and the Lipschitz continuity of $f$, the estimates for $S_e^n$ and $T_e^n$ are obtained by the usual Taylor expansions (see Eq.(20)): $|S_e^n| \le \lambda\|v\|_{2,\infty} (h + \Delta t)$, and $|T_e^n| \le C\|v\|_{2,\infty} (h + \Delta t)$, where $\|v\|_{2,\infty}$ denotes it $C^2$ uniform norm. Standard results (see [2]) are used to estimate the terms $R_{Ke}^n$ in all cases except for $e = K|L \in \mathcal{E}^\Sigma$. In this case,

$$ \begin{aligned} \hat{F}_{Ke}^n &= \int_e \sigma_0 \nabla u_0(t^n, \cdot) \cdot n_{Ke} = |e| \sigma_0 \nabla u_0(t^n, x_e) \cdot n_{Ke} + |e| A_K \\ &= +\tau_{Ke}(u_0(t^n, x_e) - u_0(t^n, x_K)) + |e| (A_K + B_K) \quad (25) \\ &= -\tau_{Le}(u_1(t^n, x_e) - u_1(t^n, x_L)) + |e| (A_K - B_L), \end{aligned} $$

because $\sigma_0 \nabla u_0(t^n, x_e) \cdot n_{Ke} = -\sigma_1 \nabla u_1(t^n, x_e) \cdot n_{Le}$. Here, $|A_K| \le M\|u_0\|_{2,\infty}h$ and $|B_K| \le M\|u_0\|_{2,\infty}h$ (resp. $|B_L| \le M\|u_1\|_{2,\infty}h$) by Taylor expansion, with $M = \max(\sigma_0, \sigma_1)$. Instead of Eq. (11), we have that $\tau_{Ke}(u_0(t^n, x_e) - \bar{u}_K^n) + \tau_{Le}(u_1(t^n, x_e) - \bar{u}_L^n) = -|e|(B_K + B_L)$ and $\bar{v}_e^n = u_1(t^n, x_e) - u_0(t^n, x_e)$ (recall that $\bar{u}_K^n = u_0(t^n, x_K)$ and $\bar{u}_L^n = u_1(t^n, x_L)$). Hence we find that

$$ u_0(t^n, x_e) = \frac{\tau_{Ke}\bar{u}_K^n + \tau_{Le}\bar{u}_L^n - \tau_{Le}\bar{v}_e^n}{\tau_{Ke} + \tau_{Le}} - \frac{|e|}{\tau_{Ke} + \tau_{Le}}(B_K + B_L). \quad (26) $$

The discrete flux associate to $\bar{u}_K^n = u_0(t^n, x_K)$ and $\bar{u}_L^n = u_1(t^n, x_L)$ is

$$ \bar{F}_{Ke}^n = -\bar{F}_{Le}^n = \tau_{Ke}(\bar{u}_{K,e}^n - \bar{u}_K^n) = -\tau_{Le}(\bar{u}_{L,e}^n - \bar{u}_L^n), \quad (27) $$

where $\bar{u}_{K,e}^n$ and $\bar{u}_{L,e}^n$ are defined as in (11) (with $\bar{v}_e^n = u_1(t^n, x_e) - u_0(t^n, x_e)$). In view of Eq. (26), we have $u_0(t^n, x_e) = \bar{u}_{K,e}^n - \frac{|e|}{\tau_{Ke}+\tau_{Le}}(B_K + B_L)$.

At last, we substract (27) to (25), and use the previous remark, to obtain the last estimate: $R_{Ke}^n = \frac{\tau_{Ke}}{|e|}(u_0(t^n, x_e) - \bar{u}_{K,e}^n) + (A_K + B_K) = \frac{\tau_{Le}B_K - \tau_{Ke}B_L}{\tau_{Ke}+\tau_{Le}} + A_K \le \max(|B_K|, |B_L|) + |A_K| \le 2M \max(\|u_0\|_{2,\infty}, \|u_1\|_{2,\infty})h$.
# 5 Conclusion and perspectives

In conclusion, the TPFA generalizes naturally to the EMI model. An error estimate was obtained on the transmembrane voltage $v$ on the interface $\Sigma$, though under strong regularity assumptions. The convergence analysis, without error estimate, might be easily deduced from the existence proof from [1], and obtained with minimal regularity of the data. We believe that the extension to many cells ($N_C > 1$) is straightforward (as in [1]), and also that the consistency estimates generalize to functions $u_i \in H^2(\Omega_i)$ and $v \in H^1(0, T; L^2(\Sigma))$.

In practice, the scheme is intended to be used for simulating billions of cardiac cells, using HPC solutions developed within the EuroHPC MICROCARD consortium. In reality, the ODEs from Eq. (2) are coupled with large sets of nonlinear equations that need time-integration specific to cardiac models. The FVM simplifies the implementation of this coupling, by introducing the voltage $v$ as an explicit unknown. Moreover, a Dirichlet to Neumann operator on $\Sigma$ can be introduced to rewrite the problem as an evolution equation on the membrane only. This gives rise to approximation with boundary element methods. Again, the discrete Dirichlet to Neumann operator may be reconstructed with the FVM, as in Eq. (16).

On a short term, we plan to implement the scheme, and compare it to finite elements and a boundary element discretizations on simple 2D and 3D test cases. Afterwards, we would like to generalize the scheme to a FV-like method that is more robust with respect to the available meshes.

*Acknowledgement.* This project is funded by the European High-Performance Computing Joint Undertaking EuroHPC, under grant agreement No 955495.

# References

1. Becue, P.E.: Modélisation et simulation de l'électrophysiologie cardiaque à l'échelle microscopique. Theses, Université de Bordeaux (2018). URL https://theses.hal.science/tel-02019648

2. Eymard, R., Gallouët, T., Herbin, R.: Finite volume methods. Handbook of numerical analysis **7**, 713–1018 (2000)

3. Haïssaguerre, M., al: Localized structural alterations underlying a subset of unexplained sudden cardiac death. Circulation: Arrhythmia and Electrophysiology **11.7** (2018)

4. Mori, Y., Peskin, C.S.: A numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes. Communications in Applied Mathematics and Computational Science **4**(1) (2009). DOI 10.2140/camcos.2009.4.85

5. Tveito, A., Jæger, K.H., Kuchta, M., Mardal, K.A., Rognes, M.E.: A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Frontiers in Physics **5** (2017). DOI 10.3389/fphy.2017.00048