# Phase-field model of bilipid membrane electroporation

Pedro Jaramillo-Aguayo, Annabelle Collin, Clair Poignard

► To cite this version:

Pedro Jaramillo-Aguayo, Annabelle Collin, Clair Poignard. Phase-field model of bilipid membrane electroporation. Journal of Mathematical Biology, 2023, 87 (18), 10.1007/s00285-023-01956-y. hal-04145549

HAL Id: hal-04145549

https://inria.hal.science/hal-04145549

Submitted on 29 Jun 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Phase-field model of bilipid membrane electroporation

Pedro Jaramillo-Aguayo, Annabelle Collin, Clair Poignard

## Abstract

This paper proposes a new model of membrane electropermeabilisation that combines the water content of the membrane and the transmembrane voltage. Interestingly, thanks to a well defined free-energy of the membrane, we somehow generalise the seminal approach of Chizmadzhev, Weaver and Krassowska, getting rid of the geometrical cylindrical assumption upon which most of the current electroporation models are based. Our approach is physically relevant and we recover a surface diffusion equation of the lipid phase proposed by Leguèbe *et al.* in a previous phenomenological model. We also perform a fine analysis of the involved nonlocal operators in two simple configurations (a spherical membrane and a flat periodic membrane) that enables us to compare the time constants of the phenomenon in spherical and flat membranes. An accurate splitting scheme combined with Fast Fourier Transforms is developed for efficient computations of the model. Our numerical results enable us to make a link between the molecular dynamics simulations of membrane permeabilisation and the experimental observations on vesicles and cells.

## 1 Introduction

Electroporation is a microscopic phenomenon that consists in imposing high and short electric pulses to biological cells in order to weaken the structure of the plasma membrane. Biologically active molecules that otherwise cannot diffuse through the membrane (e.g. hydrophilic compounds such as bleomycin or DNA) may then spread through the cell membrane. If the pulses are short enough the process is reversible: the membrane is not destroyed and reseals within minutes. The cell therefore can internalise external active molecules without losing its viability. The interest for the phenomenon has increased constantly until recently with the emergence of therapeutic strategies based on it, in oncology and also in cardiology.

Even though the phenomenon has been discovered in the late 60's, the ways the membranes become permeable by the effect of an intense electric field is still not well understood. In the late 90's, the teams of Chizmadzhev and Weaver [9, 39] proposed a description of the emergence of water pores in a sea of lipids, under very constraint geometries (the pores are cylindrical). Then Krassowska and Neu proposed a Smoluchowski equation for the population of cylindrical pores $n(r, t)$ of radius $r > 0$ at the time $t$, and they derived an asymptotic analysis to link the transmembrane potential to the total density of pores $N(t) = \int_{\mathbb{R}^+} n(r, t)dr$, which satisfies an ordinary differential equation [30, 10]. From these times, only slight modifications of the models have been proposed by Weaver [40, 34]. Kavian *et al.* proposed then a phenomenological version of the Krassowska and Debruin model, by limiting the number of parameter [18]. Leguèbe *et al.* proposed in [22] an extension by introducing the surface reaction-diffusion of oxidised lipids. This last model seems closer to the observations but it is phenomenological and lacks of any physical basis.
The aim of this paper is to propose a new model of membrane electroporation that combines membrane water content (initially 0 for non porated membrane) and transmembrane voltage. Our approach enables to make a link between the physics-based membrane free-energy approaches of Chizmadzhev and Weaver and the reaction diffusion model of the phenomenological model of Leguèbe et al. In particular, our model, which consists in Allen-Cahn equation for the water content of the membrane, and a nonlocal differential equation on the transmembrane voltage (TMV) is a generalisation of the previous approaches in the context of phase ordering kinetics.

More precisely, given the free-energy $\mathcal{E}$ of the membrane as a functional depending on the water content of the membrane $\phi$ and the transmembrane voltage $v$, our model of membrane electropermeabilisation is a system of time evolution non local equation on the membrane $\Gamma$:

$$ \left\{ \begin{array}{l} \partial_t \phi = -\alpha \frac{\partial \mathcal{E}}{\partial \phi} \quad \forall x \in \Gamma, \forall t > 0, \\ C_m(\phi) \partial_t v + (S_m(\phi) + \Lambda)v = G, \quad \forall x \in \Gamma, \forall t > 0, \\ v|_{t=0} = v_{\diamond}, \quad \phi|_{t=0} = \phi_{\diamond}, \text{ on } \Gamma, \end{array} \right. \quad (1) $$

where $\frac{\partial \mathcal{E}}{\partial \phi}$ is the Fréchet derivative of the functional $\mathcal{E}$ with respect to $\phi$, $\Lambda$ a pseudodifferential elliptic operator of order 1, which is a combination of Dirichlet to Neumann maps and $C_m$ and $S_m$ the capacitance and conductance of the membrane. This full model is presented in details in Section 3 but before that, we will start in Section 2 to present the choice of the free-energy $\mathcal{E}$ as a function of $\phi$ and $v$. In particular we somehow generalise the approach of Chizmadzhev and Weaver, getting rid of the geometrical cylindrical assumption. We then study the phase ordering Allen-Cahn model for the evolution of water content in the membrane. Section 3 focuses on the electric part of the model, that is the non local equation on the TMV.

We perform a fine analysis of the involved Dirichlet-to-Neumann nonlocal operators in two simple configurations (a spherical membrane and a flat periodic membrane) that enables us to compare the time constants of the phenomenon in spherical and flat membranes. These two simple geometrical configurations are of high importance because on the one hand, the shapes of cell in suspension and lipid vesicles in suspension are close to a sphere, while in some experiments flat film of lipid is used to investigate electroporation [12]. In addition most of the dynamic simulations on electroporation are performed in a flat geometrical setting [36, 37, 6, 35]. Section 4 is devoted to the linear stability analysis of the coupled problem. In Section 5, we present the choice of the parameters. Interestingly, we are able to make a link between our approach and the approach of Chizmadzhev and Weaver in terms of free-energy of the membrane. We also present the models of the membrane capacitance and conductance. Section 6 is devoted to the numerical strategy to solve the coupled problem in the setting of a flat membranes. The concluding section presents the perspectives of this research.

# 2 A mathematical model of membrane electropermeabilization

## 2.1 Free-energy of membrane subjected to a voltage

The membrane is considered as a combination of two mutually exclusive phases (lipid membrane and water filled pores) and its state is described through the continuous phase order parameter $x \in \Gamma \mapsto \phi(x) \in [0, 1]$, see Cahn [8, 5]. This order parameter is somehow related to the volume fraction of water that enters the membrane thanks to electroporation. The state $\phi = 0$ represents
the ideal pure lipid phase and $\phi = 1$ the ideal pure water phase. The free-energy $\mathcal{E}$ of the membrane subjected to a transmembrane voltage $v$ is the functional energy, the so-called Gibbs energy, given in [19]

$$ \mathcal{E}(\phi, v) = \frac{\kappa}{2} \int_{\Gamma} |\nabla \phi|^2 ds + \int_{\Gamma} \mathcal{W}_m(\phi) ds - \frac{1}{2} \int_{\Gamma} C_m(\phi) v^2 ds, \quad (2) $$

where $\kappa > 0$ is linked with the interfacial tension [15] and $\mathcal{W}_m$ is the double-well potential energy that describes the stable states of the membrane. The term $|\nabla \phi|^2$ regularises $\phi$ and helps control the thickness of the water-lipid interface, as explained by Bray [5].

From the physical view point, the gradient term in (2) describes the fact the interfaces water-lipid cost energy, and then two neighbouring pores will tend to merge to minimize the energy. This term enables thus to describe the interactions between neighbouring pores, which was not described by Weaver and Chizmadzhev approach. The term $\frac{1}{2}C_m(\phi)v^2$ is the electrostatic energy. It is worth noting that Weaver and Chizmadzhev considered the same electrostatic energy, which affects the potential energy to favor the phase $\phi = 1$ as described by Figure 1. Weaver and Chizmadzhev proposed to use a linear capacitance $C_m(\phi) = C_l(1 - \phi)C_w\phi$, $C_{l,w}$ being the capacitance of pure lipid and water phases respectively [1, 9, 39]. However the linear stability analysis provided in Section 4 shows that linear capacitance prevents the emergence of instabilities, as observed in a different context by Fraggedakis et al. in [15]. In this paper we choose to use the mixture model of Looyenga [26]

$$ C_m: \quad \phi \mapsto \frac{\epsilon_0}{h} \left( \left[ \epsilon_l^{1/3} + \phi(\epsilon_w^{1/3} - \epsilon_l^{1/3}) \right]^3 \vartheta_1(\phi) + \epsilon_w \vartheta_2(\phi) \right), \quad (3) $$

where

$$ \vartheta_i(\phi) = \frac{1 + \tanh(k_i(\phi - \phi_i^{\text{th}}))}{2}, \quad i = 1, 2, $$

are smooth cutoff functions. In the following, we set $k_1 = -15$, $k_2 = 13$ and $\phi_1^{\text{th}} = \phi_2^{\text{th}} = 1$.

As the water volume fraction is not conserved (defects can be created and disappear from the membrane), we consider the non-conserved dynamics associated to this energy functional (also called model-A by Bray [5]). The evolution of $\phi$ is determined by the $L^2$-gradient flow associated with the energy functional above. It corresponds to the Euler-Lagrange equation for the energy $\phi \mapsto \mathcal{E}(\phi, v)$. In other words, for a given kinetics coefficient $\alpha > 0$ – also called phase field mobility – the order parameter $\phi$ satisfies the following Allen-Cahn equation on the membrane surface $\Gamma$ [3, 5, 8]:

$$ \begin{aligned} \partial_t \phi - D_0 \Delta \phi &= -\alpha \frac{\partial \mathcal{W}_m}{\partial \phi}(\phi) + \frac{\alpha}{2} \frac{\partial C_m}{\partial \phi}(\phi) v^2, \quad \forall t > 0, \\ \phi(0, \cdot) &= \phi^0(\cdot), \end{aligned} $$

where $D_0 := \alpha\kappa$ is the lateral diffusion coefficient of lipids in $[\text{m}^2 \cdot \text{s}^{-1}]$, and $\alpha$ is in $[\text{m}^2 \cdot \text{J}^{-1} \cdot \text{s}^{-1}]$.

## 2.2 On the choice of the double well potential $\mathcal{W}_m$

The simplest double well potential<fnref n="1" /> considered in Bray’s book is $\phi \to \phi^2(1 - \phi)^2$, for which the pure phases $\phi = 0$ and $\phi = 1$ are global minima and the mixed phase $\phi = 1/2$ is the local maximum [5].

<footnote n="1">¹It is the polynomial with the lowest degree that describes the phase separation phenomenon.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Influence of a given constant transmembrane voltage (TMV) $v$ on the potential energy $W_m(\phi) - C_m(\phi)v^2/2$. One can see that the TMV tends to tilt the double-well potential to favor the phase $\phi = 1$, favoring thus the entry of water in the membrane. The parameters are given in Tables 1 and 2.</figcaption></figure>
The energy barrier to pass from the phase $\phi = 0$ to the phase $\phi = 1$ is then defined as the height $W_m(1/2) - W_m(0)$. Even though the qualitative behavior is included in this simple form of potential, it is necessary to introduce a parametrisable double well potential to make the link with the free-energy introduced by Weaver and Chizmadzhev [9, 39].

Throughout the paper, the double-well potential $W_m$ is set as :

$$W_m(\phi) := a_1\phi^2(1-\phi)^2 + a_2(\phi + \frac{1}{2})(\phi - 1)^2, \quad \forall \phi \in [0, 1], \qquad (5)$$

where $a_1 > \frac{3}{2}|a_2|$, so that $W_m$ satisfies

$$W_m(0) = \frac{a_2}{2}, \quad W_m(1) = 0, \quad W'_m(1) = W'_m(0) = 0, \quad W''_m(0) > 0, \quad \text{and} \quad W''_m(1) > 0.$$

The above energy potential describes the fact that the pure lipid ($\phi = 0$) and pure water ($\phi = 1$) phases are the only two stable phases. The energy barrier to pass from lipid to water equals $W_m(1/2) - a_2/2$. Note that unlike the symmetrical potential of Bray's book, the energy barrier to pass from the pure water phase to the pure lipid phase equals $W_m(1/2)$. Thus it requires more energy to pass from the pure water to the pure lipid phase than the opposite. Section 5 is dedicated to a fine analysis of the energy to choose the parameters $a_1$ and $a_2$ in adequation with Weaver and Chizmadzhev's work.

It is however important to keep in mind that $a_2$ mainly controls the value of the first local minimum ($W_m(0) = a_2/2$) and $a_1$ the size of the local maximum (at $\phi_0 := 1/2 - 3a_2/(4a_1)$) which is the height of the energy barrier to pass from 0 to 1.

The following lemma states the well-posedness of the phase-field model of the membrane submitted to a given transmembrane voltage.

**Lemma 1.** Let $s \ge 3$, $T > 0$, $C_m \in C^{s+2}(\mathbb{R})$, $v \in C([0, T[, H^s(\Gamma))$ and $\phi_{\diamond} \in H^s(\Gamma)$. Then there exists a unique mild solution $\phi \in C([0, T[, H^s(\Gamma))$ to

$$\partial_t \phi - D_0 \Delta \phi = -4\alpha a_1 \phi(\phi - 1)(\phi - 1/2 + \frac{3a_2}{4a_1}) + \frac{\alpha}{2} C'_m(\phi) v^2, \quad \text{in } ]0, +\infty[ \times \Gamma, \qquad (6a)$$

$$\phi|_{t=0} = \phi_{\diamond}, \quad \text{on } \Gamma. \qquad (6b)$$

Even though the result of the above lemma is quite standard, we give the proof to make the paper self-contained.

*Proof.* We use the semi-group theory combined with Picard fixed-point argument to prove this result. Let $S(t)$ be the contractive semigroup in $H^s(\Gamma)$ generated by the operator $f \mapsto D_0\Delta f$, then we can rewrite our problem in its mild formulation as

$$\phi(t) = S(t)\phi_{\diamond} + \int_{0}^{t} S(t - \tau)F(\tau, \phi) d\tau$$

where

$$F(\tau, \phi) := F_1(\phi) + F_2(\tau, \phi),$$

$$F_1(\phi) := -4\alpha a_1 \phi(\phi - 1)(\phi - \frac{1}{2} + \frac{3a_2}{4a_1}),$$

$$F_2(\tau, \phi) := \frac{\alpha}{2} C'_m(\phi) v^2.$$
We know that $\phi \mapsto F(\tau, \phi)$ is a Lipschitz function (uniformly on $[0, T]$) on bounded sets of $H^s(\Gamma)$.
In fact, $H^s(\Gamma)$ is an algebra given that $s \ge 3$ and $\dim(\Gamma) = 2$ and so $F_1$ is Lipschitz on bounded sets (it is a polynomial). A similar argument shows that $\phi \mapsto \frac{\alpha}{2}C'_m(\phi)v^2$ is Lipschitz if $\phi \mapsto C'_m(\phi)$ is Lipschitz on bounded sets of $H^s(\Gamma)$, which occurs thanks to the regularity of $C_m$.
Let $\mathcal{E} := C^0([0, T], H^s(\Gamma))$ provided with the following norm

$$ \|u\|_{\mathcal{E}} := \sup_{\tau \in [0, T]} e^{-\beta\tau} \|u(\tau)\|_{H^s(\Gamma)}, \quad \text{where } \beta > 0, $$

and consider the following function

$$ \Psi : \mathcal{E} \to \mathcal{E} $$

$$ \phi \mapsto S(t)\phi_{\diamond} + \int_{0}^{t} S(t - \tau)F(\tau, \phi)d\tau. $$

Let $C$ be a generic constant which takes other constants into account then for $\phi_1, \phi_2 \in B(0, R) \subset \mathcal{E}$,

$$ \begin{aligned} \|\Psi(\phi_1) - \Psi(\phi_2)\|_{H^s(\Gamma)} &\le \int_0^t \|F(\tau, \phi_1) - F(\tau, \phi_2)\|_{H^s(\Gamma)} d\tau, \\ &\le C \int_0^t e^{\beta\tau} e^{-\beta\tau} \|\phi_1 - \phi_2\|_{H^s(\Gamma)} d\tau, \\ &\le C \frac{1}{\beta} e^{\beta t} \|\phi_1 - \phi_2\|_{\mathcal{E}}, \end{aligned} $$

and so

$$ \|\Psi(\phi_1) - \Psi(\phi_2)\|_{\mathcal{E}} \le \frac{C}{\beta} \|\phi_1 - \phi_2\|_{\mathcal{E}}, $$

where the constant $C$ depends on $v$, $O_e$ and $O_c$. By taking $\beta$ large enough, $\Psi$ is a contractive map and we can conclude the proof. $\square$

## 2.3 Qualitative properties of the membrane order parameter at null TMV

In this section we summarize the existing mathematical results on the solution to Problem (6) (with $v \equiv 0$). There is an extensive literature studying the qualitative mathematical properties of (6). One can cite, for example [3, 13, 2]. All the following results are related to the following equation

$$ \partial_t \phi - D_0 \Delta \phi = -4\alpha a_1 \phi(\phi - 1)(\phi - 1/2 + \frac{3a_2}{4a_1}), \quad \text{in } ]0, +\infty[ \times \Gamma, \quad (7a) $$

$$ \phi|_{t=0} = \phi_{\diamond}, \quad \text{on } \Gamma. \quad (7b) $$

**Property 2.** If the initial state of the membrane $\phi_{\diamond}$ is smooth enough $(C^2(\Gamma)$ for example) and verifies $0 \le \phi_{\diamond}(s) \le 1$ for almost any $s \in \Gamma$ then, the solution $\phi$ to (7) satisfies the same bounds

$$ 0 \le \phi \le 1, \quad \text{a.e on } (0, +\infty) \times \Gamma. $$

The next property is due to Farina et al., see [14]. To reference it, we first need to introduce the following notion of stability.
**Definition 3.** We say that a stationary solution $\phi$ to (7) is stable if

$$ \int_{\Gamma} (\kappa |\nabla \xi|^2 + W''_{m}(\phi) \xi^2) dx \geq 0 $$

for every smooth function $\xi \in C^{\infty}(\Gamma)$. This quadratic form is called the second variation of the energy functional (2) (with $v \equiv 0$).

**Property 4.** Let $\Gamma$ be a smooth closed compact manifold (with nonnegative Ricci curvature), then the only stable solutions of (7) are constant.

**Remark 5.** This means that for the considered potential $W_m$, any porated membrane will either reseal to a non porated state ($\phi \equiv 0$) or disappear ($\phi \equiv 1$).

The following property is due to Alfaro et al. Before stating the result, we first introduce some notations.

**Notation 6.** Let $(a_n)_{n \in \mathbb{N}}$, $(b_n)_{n \in \mathbb{N}} \in \mathbb{R}$ be 2 sequences of $\mathbb{R}$.

• We write

$$ a_n \lesssim b_n $$

if there is a constant $C > 0$ such that for all $n \geq 0$, we have $a_n \leq Cb_n$.

• Similarly, we write

$$ a_n \ll b_n $$

if there is a sequence $\varepsilon_n \to_{n \to +\infty} 0$ such that $a_n \leq \varepsilon_n b_n$.

**Property 7.** Let $L > 0$ denote the characteristic length scale of $\Gamma$ and let $(D_{0,n})_{n \in \mathbb{N}}$, $(a_{1,n})_{n \in \mathbb{N}}$, $(a_{2,n})_{n \in \mathbb{N}}$ be three sequences strictly positive sequences. Let $\phi_n$ be the solution to the rescaled equation

$$ \partial_t \phi_n - \Delta \phi_n = \frac{-4a_{1,n}\alpha L^2}{D_{0,n}} \phi_n(\phi_n - 1) \left( \phi_n - \frac{1}{2} + \frac{3a_{2,n}}{4a_{1,n}} \right), \quad (8) $$

$$ \phi_n(t = 0, \cdot) = \phi_{\diamond}. \quad (9) $$

Assume that $(D_{0,n})_{n \in \mathbb{N}}$, $(a_{1,n})_{n \in \mathbb{N}}$ and $(a_{2,n})_{n \in \mathbb{N}}$ satisfy the following asymptotic behaviour

$$ \frac{3a_{2,n}}{4a_{1,n}} \lesssim \sqrt{\frac{D_{0,n}}{4a_{1,n}\alpha L^2}} \ll 1. $$

Then one has the following properties as $n \to \infty$:

1. Given any smooth initial condition $\phi_{\diamond} \in [0, 1]$, the typical profile of a solution to (8) involves different regions where $\phi_n \simeq 1$ and $\phi_n \simeq 0$ and transition interfaces between them of size

$$ \delta_n \lesssim \sqrt{\frac{D_{0,n}}{4a_{1,n}\alpha}} \quad [m]. \quad (10) $$

The time needed for the solution associated to $\phi_0$ to take above description is

$$ \mathcal{T}_n \lesssim \frac{1}{4a_{1,n}\alpha} \left| \log \left( \frac{D_{0,n}}{4a_{1,n}\alpha L^2} \right) \right| \quad [s]. \quad (11) $$

2. The evolution of the transition interfaces can be described by the evolution of the level set $\{x \in \Gamma | \phi(x) = 1/2\}$. In the asymptotic regime described above, the motion of this interface is primarily determined by mean-curvature flow [13, 2].

In the next section, we present the equation satisfied by the transmembrane potential.
# 3 Transmembrane voltage in a membrane

## 3.1 Electric field around membranes

Throughout the paper, $\Gamma$ is a closed 2D-surface of without boundary. We denote by $\Omega$ the bounded domain of $\mathbb{R}^3$ in which $\Gamma$ is embedded, and let $O_e$ and $O_c$ be the 2 connected components subsets of $\Omega$ and separated by $\Gamma$. Let $\sigma_e$ (resp. $\sigma_c$) be the constant conductivities of $O_e$ (resp. $O_c$), and let $C_m : \phi \to C_m(\phi)$, $S_m : \phi \to S_m(\phi)$ be the surface capacitance and conductance of $\Gamma$, strictly positive bounded smooth functions of the order parameter $\phi$.

Let $v_{\diamond}$ be a regular enough function of $\Gamma$. We assume that it belongs at least to $H^1(\Gamma)$. The electric potential $U$ around the membrane $\Gamma$ verifies the following partial differential equation (PDE)

$$ [U]_{\Gamma}(t=0, \cdot) = v_{\diamond}(\cdot), \quad \text{on } \Gamma, \qquad (12a) $$

and for any time $t > 0$:

$$ \left\{ \begin{array}{ll} \nabla \cdot (\sigma_e \nabla U) = 0, & \text{in } O_e, \\ \nabla \cdot (\sigma_c \nabla U) = 0, & \text{in } O_c, \\ U(t, \cdot) = g(t, \cdot), & \text{on } \partial\Omega, \\ \sigma_e \vec{n}_c \cdot \nabla U|_{\Gamma^-} + \sigma_e \vec{n}_e \cdot \nabla U|_{\Gamma^+} = 0, & \text{on } \Gamma, \\ \sigma_e \vec{n}_e \cdot \nabla U|_{\Gamma^+} = C_m(\phi) \partial_t [U]_{\Gamma} + S_m(\phi) [U]_{\Gamma}, & \text{in } \Gamma, \end{array} \right. \qquad (12b) $$

where $[U]_{\Gamma} = U|_{\Gamma^-} - U|_{\Gamma^+}$ and $U|_{\Gamma^{\pm}} : x \mapsto \lim_{\tau \to 0^+} U(x \mp \tau \vec{n}(x))$ and $\vec{n}(x) := \vec{n}_c(x)$ denotes the unitary normal vector to $\Gamma$ directed towards $O_c$ at any point $x \in \Gamma$, while $\vec{n}_e(x) = -\vec{n}_c(x)$ is the normal vector to $\Gamma$ directed towards $O_e$.

**Remark 8.** It is worth noting that only the initial value of the transmembrane voltage $[U]_{\Gamma}(t=0, \cdot)$ is required and not the value of the initial potential everywhere in the domain. This is due to the fact that in the domains $O_e$ and $O_c$, we have neglected the displacement currents and only the conductive currents are considered, leading to an elliptic equation in the inner and outer domains.

Interestingly, the above volume PDE can be rewritten in terms of Dirichlet-to-Neumann operators. Denote by $\Lambda_c$, $\Lambda_e$ and $\Lambda_o$ the three following Dirichlet-to-Neumann operators

$$ \Lambda_c : H^{1/2}(\Gamma) \to H^{-1/2}(\Gamma) \qquad \text{, where } v_c \text{ is the solution to} \quad \left\{ \begin{array}{ll} \nabla \cdot (\sigma_c \nabla v_c) = 0, & \text{in } O_c, \\ v_{c|_{\Gamma}} = f, & \end{array} \right. \qquad (13a) $$

$$ \Lambda_e : H^{1/2}(\Gamma) \to H^{-1/2}(\Gamma) \qquad \text{, where } v_e \text{ is the solution to} \quad \left\{ \begin{array}{ll} \nabla \cdot (\sigma_e \nabla v_e) = 0, & \text{in } O_e, \\ v_{e|_{\Gamma}} = f, & \end{array} \right. \qquad (13b) $$

$$ \Lambda_o : H^{1/2}(\partial\Omega) \to H^{-1/2}(\Gamma) \qquad \text{, where } v_b \text{ is the solution to} \quad \left\{ \begin{array}{ll} \nabla \cdot (\sigma_e \nabla v_b) = 0, & \text{in } O_e, \\ v|_{\Gamma} = 0, & \end{array} \right. \qquad (13c) $$
Following [18], it is equivalent to solve the volume equation (12) for $U$ or the following nonlocal equation on the surface $\Gamma$ for the transmembrane voltage (TMV), $v = [U]_\Gamma$,

$$C_m(\phi)\partial_t v + (S_m(\phi) + \Lambda)v = G, \quad (14a)$$

$$v(t = 0, \cdot) = v_{\diamond}(\cdot), \quad (14b)$$

where

$$\Lambda = \Lambda_c(\Lambda_e + \Lambda_c)^{-1}\Lambda_e \quad (14c)$$

$$G = \Lambda_c(\Lambda_e + \Lambda_c)^{-1}\Lambda_{og}. \quad (14d)$$

**Remark 9** (Invertiblity of $\Lambda_e + \Lambda_c$). The fact that the operator $\Lambda_e + \Lambda_c$ is invertible is proven in Kavianet al. (see Lemma 8 of [18]). Its comes from the fact that $\Lambda_e$ is invertible that $\Lambda_c$ is a non negative self-adjoint operator.

**Proposition 10.** Problem (14) is equivalent to Problem (12) in the following sense:

* If $U$ is solution to (12), then $[U]_\Gamma : \mathbb{R}^+ \times \Gamma \to \mathbb{R}$ is solution to (14).

* If $v$ solution to (14) then the piecewise function $\tilde{U} : \mathbb{R}^+ \times O_e \cup O_c \to \mathbb{R}$, defined on $O_e$ as the solution to

$$\left\{
\begin{array}{l}
\nabla \cdot (\sigma_e \nabla \tilde{U}) = 0, \forall x \in O_e, t \ge 0, \\
\tilde{U}_{|_{\partial\Omega}} = g, \forall t \ge 0, \\
\sigma_e \partial_{n_e} \tilde{U}_{|_{\Gamma^+}} = C_m(\phi) \partial_t v + S_m(\phi)v, \forall x \in \Gamma, t > 0, \\
\tilde{U}_{|_{\Gamma^+}}(0, \cdot) = -(\Lambda_e + \Lambda_c)^{-1}(\Lambda_0 g + \Lambda_c v_{\diamond}), \forall x \in \Gamma.
\end{array}
\right.$$

and defined on $O_c$ as the solution to

$$\left\{
\begin{array}{l}
\nabla \cdot (\sigma_c \nabla \tilde{U}) = 0, \forall x \in O_c, t \ge 0 \\
\tilde{U}_{|_{\Gamma^-}} = (\Lambda_c + \Lambda_e)^{-1}(\Lambda_e v - \Lambda_0 g), \forall x \in \Gamma, t \ge 0.
\end{array}
\right.$$

is solution to (12).

*Proof.* This proof is taken from [18]. Using the definition of our Dirichlet-to-Neumann operators we can see that $U$ is a solution to (12) if and only if it verifies the following equations

$$\left\{
\begin{array}{l}
-\Lambda_c U_{|_{\Gamma^-}} = C_m(\phi)\partial_t [U]_\Gamma + S_m(\phi)[U]_\Gamma, \\
\Lambda_c U_{|_{\Gamma^-}} + \Lambda_e U_{|_{\Gamma^+}} + \Lambda_0 g = 0.
\end{array}
\right.$$

Applying the operator $(\Lambda_e + \Lambda_c)^{-1}$ to the second equation and rearranging the terms results in the following relation

$$-U_{|_{\Gamma^-}} = -(\Lambda_c + \Lambda_e)^{-1}(\Lambda_e [U]_\Gamma - \Lambda_0 g).$$

Injecting this expression into the first equation shows the first point, with initial condition $v_{\diamond} := [U_0]_\Gamma$. The second point follows similarly. $\square$
### 3.1.1 Well-posedness of the transmembrane potential

The following result extends the well-posedness results of Kavian et al. [18] to time-varying and space dependent capacitance.

**Lemma 11.** Let $T > 0$, $s \ge 3$, $\beta_0, \beta_1 \in (0, 1]$, $G \in C^{\beta_0}([0, T[, H^s(\partial\Omega))$ and $\phi \in C^{1, \beta_1}([0, T[, H^s(\Gamma))$. Then for every initial condition $v_\circ \in L^2(\Gamma)$ there exists a unique classical solution

$$v \in C^1([0, T], L^2(\Gamma)) \cap C((0, T], D(\Lambda))$$

to

$$C_m(\phi)\partial_t v + (S_m(\phi) + \Lambda)v = G \\ v(t=0) = v_0.$$

*Proof.* This result is a direct application of Theorem 6.1 and Theorem 7.1 from Pazy book [31, Chapter 5]). In order to apply these results, let us first define

$$\underline{v} := \sqrt{C_m(\phi)}v,$$

which satisfies

$$\partial_t \underline{v} + \underbrace{\left( \frac{S_m(\phi) - \frac{1}{2} C'_m(\phi) \partial_t \phi}{C_m(\phi)} \right)}_{=:M(t)} \underline{v} + \underbrace{\frac{1}{\sqrt{C_m(\phi)}} \Lambda \frac{1}{\sqrt{C_m(\phi)}}}_{=:\tilde{\Lambda}(t)} \underline{v} = \frac{G}{\sqrt{C_m(\phi)}} \underbrace{=}_{=:f(t)}.$$

For all $t \ge 0$ the operator $\tilde{\Lambda}(t)$ is selfadjoint in $L^2(\Gamma)$ and monotone ($\Lambda$ is monotone, see [18]). The function $M(t)$ is bounded and so we can always change the unknown we are solving for by $\underline{v}_k(t) := e^{-kt}\underline{v}(t)$ for large enough $k > 0$ and just consider $M(t) > 2$ and bounded. To match the same notation used in [31], we define $A(t) = M(t) + \tilde{\Lambda}(t)$. This operator is also selfadjoint and monotone. Once this setting in place, we need to verify three conditions to apply both theorems:

* The domain of $D(A(t))$ is dense in $L^2(\Gamma)$ and is independent of $t \in [0, T]$.

* For $t \in [0, T]$, the resolvent $R(\lambda : A(t))$ exists for all $\lambda \in \mathbb{C}$, $Re(\lambda) \le 0$ and there exists $C > 0$ such that

$$\|R(\lambda : A(t))\|_{L^2(\Gamma)} \le \frac{C}{|\lambda| + 1}, \quad \forall t \in [0, T], \quad Re(\lambda) \le 0.$$

* There exists a constant $C'$ and $0 < \zeta \le 1$ such that

$$\|(A(t) - A(s))A(\tau)^{-1}\|_{L^2(\Gamma)} \le C'|t - s|^\zeta, \quad \forall s, t, \tau \in [0, T].$$

The first point only depends on the domain of $\tilde{\Lambda}(t)$ as $v \mapsto M(t)v$ is bounded on $L^2(\Gamma)$. As $\phi \mapsto C_m(\phi)$ is strictly positive and smooth, we obtain the following equality

$$D(A(t)) = D(\Lambda) = \{u \in H^{\frac{1}{2}}(\Gamma) | \Lambda u \in L^2(\Gamma)\},$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Two situations of membrane considered in this paper. On the left, the cell membrane separates the inner of the cell to the bath. The cell membrane $\Gamma$ is the sphere of radius $R_0$ and the boundary condition at $\partial\Omega$ describes the effect of the electric field applied. On the right, a flat bi-periodic patch of the membrane placed between two charged planes $\partial\Omega^{\pm}$.</figcaption></figure>

where the graph norm of $\Lambda$ is equivalent to the graph norm of $\tilde{\Lambda}(t)$ for all $t \in [0, T]$. In fact the domain of $\Lambda$ is $H^1(\Gamma)$ (see [18]) and so the norm equivalence results from the open mapping theorem. The second point results from the following calculation

$$ \| (\lambda - A(t))^{-1} \|_{L^2(\Gamma)} = \left\| \left( \lambda - 1 - (M(t) - 1) - \tilde{\Lambda}(t) \right)^{-1} \right\|_{L^2(\Gamma)} \le \frac{1}{\sqrt{|\lambda|^2 + 1 - \Re\lambda}} \le \frac{2}{|\lambda| + 1} $$

where we exploit the fact that $A(t) - 1$ is still a selfadjoint monotone operator. The third point is just a consequence of the Hölder continuity of $\partial_t\phi$, and the second point. Due the smoothness in time of $t \mapsto \phi(t)$, we can take $\zeta = 1$. $\square$

In the subsections, we characterize the operator $\Lambda$ in two important configurations: a spherical membrane and a flat torus, see Figure 2. These two cases are of importance, because cells in suspension are mostly round, while molecular dynamic simulations deal mostly with flat membrane. It is thus important to compare these two settings, in particular to understand the main differences in terms of order of magnitude of the time constants.

## 3.2 Spherical and flat membranes

As stated in the introduction, understanding the specificities of spherical and flat membranes is motivated by the fact that in the experiments cell and vesicle in suspension are mostly rounded, while molecular dynamic simulations are mostly performed in a flat periodic setting [12, 35]. In particular we propose in this section to perform a fine analysis of the nonlocal operator $\Lambda$ to provide quantitative criteria to compare the flat and the spherical settings, see Figure 2.
### 3.2.1 The operator $\Lambda$ in the case of a spherical membrane

We consider $\Gamma$ and $\partial\Omega$ two spheres of radius $R_0$ and $R_1$ respectively where $R_1 > R_0$ (see Figure 2). The interior of the cell $O_c = B(R_0)$ being the inner ball of radius $R_0$ and $O_e = B(R_1) \setminus (\Gamma \cup B(R_0))$ the rest of the domain.

The operator $\Lambda$ of (14)-denoted by $\Lambda^S$ in the spherical setting- can be explicitly diagonalised in terms of the eigenfunctions of the Laplace-Beltrami operator on the unit sphere denoted by $S^1$.

More precisely, the spectrum of the Laplace-Beltrami operator $-\Delta_S$ is the set $\{\ell(\ell + 1), \ell \in \mathbb{N}\}$, each eigenvalue $\ell(\ell + 1)$ being of multiplicity $2\ell + 1$, and the eigenfunctions are the well-known spherical harmonics (see Muller's book [29]).

As a result we get the following lemma.

**Lemma 12** (Eigenmodes of the operator $\Lambda^S$ for the sphere of radius $R_0$). The operator $\Lambda^S$ from Equation (14) is diagonalisable in the same basis as the sphere Laplace-Beltrami operator $-\Delta_{S^1}$. Let $\omega_{\lambda_{\Delta}^{S1}}$ be an eigenfunction of $-\Delta_{S^1}$ associated with the eigenvalue $\lambda_{\Delta}^{S1} \in \{\ell(\ell + 1), \ell \in \mathbb{N}\}$.

$$k_{\Delta}^{\pm} := \frac{-1 \pm \sqrt{1 + 4\lambda_{\Delta}^{S1}}}{2}, \quad \text{with } \lambda_{\Delta}^{S1} \in \{\ell(\ell + 1), \ell \in \mathbb{N}\},$$
$$\lambda^{c,S} := \frac{\sigma_c k_{\Delta}^{+}}{R_0},$$

$$\lambda^{e,S} := \frac{\sigma_e (k_{\Delta}^{+} - k_{\Delta}^{-} (\frac{R_1}{R_0})^{\sqrt{1+4\lambda_{\Delta}^{S1}}})}{R_0 (\frac{R_1}{R_0})^{\sqrt{1+4\lambda_{\Delta}^{S1}} - 1}},$$

$$\lambda^{o,S} := \frac{-\sigma_e}{R_0} \left(\frac{R_1}{R_0}\right)^{-k_{\Delta}^{-}} \frac{\sqrt{1 + 4\lambda_{\Delta}^{S1}}}{\left(\frac{R_1}{R_0}\right)^{\sqrt{1+4\lambda_{\Delta}^{S1}} - 1}},$$

then

$$\Lambda^S \omega_{\lambda_{\Delta}^{S1}} = \frac{\lambda^{e,S} \lambda^{c,S}}{\lambda^{c,S} + \lambda^{e,S}} \omega_{\lambda_{\Delta}^{S}}.$$

In other words, if $(\omega_{\lambda_{\Delta}^{S1}}, \lambda_{\Delta}^{S1})$ is the couple eigenvector/eigenvalue of $-\Delta_{S^1}$, then $(\omega_{\lambda_{\Delta}^{S1}}, \lambda^S)$ is a couple eigenvector/eigenvalue of the operator $\Lambda^S$ where

$$\lambda^S = \frac{\lambda^{e,S} \lambda^{c,S}}{\lambda^{c,S} + \lambda^{e,S}}.$$

**Remark 13.** As both operators $\Lambda^S$ and $-\Delta_S$ are diagonalisable in the same basis, which is a Hilbert basis of smooth functions in $L^2(S)$, the above list of eigenvalues for $\Lambda^S$ is exhaustive.

*Proof.* To prove the lemma, we diagonalise the three operators $\Lambda_o^S$, $\Lambda_e^S$ and $\Lambda_c^S$ and then conclude by using the above relations. First, we change to spherical coordinates. As $\sigma_c$ and $\sigma_e$ are constant coefficients, given an interval $I = (a, b)$ ($b > a \ge 0$), each Problem (13a)-(13b)-(13c) can be rewritten as

$$\left\{ \begin{array}{l} \left(\frac{2\partial_r}{r} + \partial_r^2\right)U + \frac{1}{r^2}\Delta_{S^1}U = 0, \quad \forall (r, \theta) \in I \times S^1, \\ U|_{\partial I \setminus \{r=0\}}(\theta) = f(\theta), \quad \forall \theta \in S^1 \end{array} \right\}$$
We take $f = \omega_{\lambda_{\Delta}^{\mathbb{S}^1}}$ and then we proceed by separation of variables. We assume the solution has the following ansatz $U(r, \theta) = h(r)\omega_{\lambda_{\Delta}^{\mathbb{S}^1}}(\theta)$. This implies that $h$ verifies the following second order ODE

$$ \frac{2}{r}h' + h'' - \frac{\lambda_{\Delta}^{\mathbb{S}^1}}{r^2}h = 0, $$

whose solutions are given by

$$ h(r) = Ar^{k_{\Delta}^{+}} + Br^{k_{\Delta}^{-}}, \text{ where } k_{\Delta}^{\pm} = \frac{-1 \pm \sqrt{1 + 4\lambda_{\Delta}^{\mathbb{S}^1}}}{2}, $$

for some coefficients $A, B \in \mathbb{R}$ to be determined. We determine the coefficients according to each Dirichlet boundary condition:

* By definition of the operator $\Lambda_c^{\mathbb{S}}$, $r \in (0, R_0)$ and $h(R_0) = 1$, hence $B = 0$ and $A = R_0^{-1}$ and therefore

$$ \Lambda_c^{\mathbb{S}}(\omega_{\lambda_{\Delta}^{\mathbb{S}^1}}) = \frac{\sigma_e k_{\Delta}^{+}}{R_0} \omega_{\lambda_{\Delta}^{\mathbb{S}^1}}. $$

* By definition of the operator $\Lambda_e^{\mathbb{S}}$, $r \in (R_0, R_1)$, $h(R_0) = 1$ and $h(R_1) = 0$, simple calculations lead to

$$ \Lambda_e^{\mathbb{S}}(\omega_{\lambda_{\Delta}^{\mathbb{S}^1}}) = \frac{\sigma_e (k_{\Delta}^{+} - k_{\Delta}^{-} (\frac{R_1}{R_0})^{\sqrt{1+4\lambda_{\Delta}^{\mathbb{S}^1}}})}{R_0 (\frac{R_1}{R_0})^{\sqrt{1+4\lambda_{\Delta}^{\mathbb{S}^1}}-1}} \omega_{\lambda_{\Delta}^{\mathbb{S}^1}}. $$

* By definition of $\Lambda_o^{\mathbb{S}}$, $r \in (R_0, R_1)$, $h(R_0) = 0$ and $h(R_1) = 1$, hence

$$ \Lambda_o^{\mathbb{S}}(\omega_{\lambda_{\Delta}^{\mathbb{S}}}) = \frac{-\sigma_e}{R_0} \left(\frac{R_1}{R_0}\right)^{-k_{\Delta}^{-}} \frac{\sqrt{1 + 4\lambda_{\Delta}^{\mathbb{S}^1}}}{\left(\frac{R_1}{R_0}\right)^{\sqrt{1 + 4\lambda_{\Delta}^{\mathbb{S}^1}} - 1}} \omega_{\lambda_{\Delta}^{\mathbb{S}}}. $$

This concludes the proof. $\square$

### 3.2.2 The operator $\Lambda$ in the case of a flat periodic membrane

We consider now the case of a flat periodic membrane $(\mathbb{R}/(LZ))^2$. More precisely, the domain $\Omega$ is defined by

$$ \Omega = (\mathbb{R}/(LZ))^2 \times (-H, H) \subset \mathbb{R}^3 $$

where $H$ and $L > 0$ are strictly positive constants.

The membrane is defined by $\Gamma = (\mathbb{R}/(LZ))^2 \times \{0\} \subset \Omega$ and $\partial\Omega = \partial\Omega^+ \cup \partial\Omega^-$ where

$$ \partial\Omega^+ = (\mathbb{R}/(LZ))^2 \times \{H\}, $$

$$ \partial\Omega^- = (\mathbb{R}/(LZ))^2 \times \{-H\}. $$

The exterior and interior of the cell are modeled as

$$ \mathcal{O}_e = \{(x, y, z) \in \Omega | z > 0\} $$

$$ \mathcal{O}_c = \{(x, y, z) \in \Omega | z < 0\} $$
respectively. In this setting, the electric potential around the membrane U satisfies (12) with the boundary condition

$$g(t, \cdot) = g^{\pm}(t, \cdot) \text{ on } \partial\Omega^{\pm},$$

where $g^+$ and $g^-$ are given smooth enough functions of $\partial\Omega^+$ and $\partial\Omega^-$ respectively. We also get the similar diagonalisation result as in the spherical sphere.

**Lemma 14** (Eigenmodes of the operator $\Lambda$ for the flat torus $(\mathbb{R}/(L\mathbb{Z}))^2$). The operator $\Lambda$ – denoted by $\Lambda^\mathbb{T}$ in the case of a flat torus – is diagonalisable in the same basis as the periodic Laplace operator $-\Delta_\mathbb{T}$. For $\mathbf{k} = (\mathbf{k}_1, \mathbf{k}_2) \in \mathbb{Z}^2$, let $4\pi^2\|\mathbf{k}\|^2/L^2$ and $\omega_\mathbf{k}(x, y) = \exp((2i\pi/L)(\mathbf{k}_1x + \mathbf{k}_2y))$ be the eigenpair of $-\Delta_\mathbb{T}$ in $\Gamma$ then

$$\Lambda^{\mathbb{T}} \omega_{\mathbf{k}} = \lambda_{\mathbf{k}}^{\mathbb{T}} \omega_{\mathbf{k}},$$

where

$$\lambda_{\mathbf{k}}^{\mathbb{T}} = \begin{cases} \frac{1}{H} \frac{\sigma_c \sigma_e}{\sigma_c + \sigma_e}, & \text{if } \mathbf{k} = 0, \\ \frac{2\pi |\mathbf{k}|}{L \tanh(\frac{2\pi}{L} H |\mathbf{k}|)} \frac{\sigma_c \sigma_e}{\sigma_c + \sigma_e}, & \text{otherwise.} \end{cases}$$

*Proof.* We proceed in the same way as in the case of the sphere. Denote by $\Lambda_c^\mathbb{T}$, $\Lambda_e^\mathbb{T}$, and $\Lambda_o^\mathbb{T}$ the operators (13) in this flat periodic setting. It is also convenient to split $\Lambda_o^\mathbb{T}$ into $\Lambda_o^{\mathbb{T},+}$ and $\Lambda_o^{\mathbb{T},-}$ defined as

$$\Lambda_o^{\mathbb{T},+}: H^{1/2}(\partial\Omega^+) \to H^{-1/2}(\Gamma), \quad g \mapsto -\sigma_e \partial_z v_{|_{\Gamma^+}}, \text{ where } v \text{ is the solution to } \begin{cases} \nabla \cdot \sigma_e \nabla v = 0, & \text{in } \mathcal{O}_e, \\ v|_{\Gamma} = 0, \\ v|_{\partial\Omega^+} = g^+ \end{cases}$$

$$\Lambda_o^{\mathbb{T},-}: H^{1/2}(\partial\Omega^-) \to H^{-1/2}(\Gamma), \quad g \mapsto \sigma_c \partial_z v_{|_{\Gamma^-}}, \text{ where } v \text{ is the solution to } \begin{cases} \nabla \cdot \sigma_c \nabla v = 0, & \text{in } \mathcal{O}_c, \\ v|_{\Gamma} = 0, \\ v|_{\partial\Omega^-} = g^- \end{cases}.$$

Then, from (14) we infer

$$\Lambda^{\mathbb{T}} = \Lambda_c^{\mathbb{T}}(\Lambda_e^{\mathbb{T}} + \Lambda_e^{\mathbb{T}})^{-1}\Lambda_e^{\mathbb{T}}, \text{ and } G = \frac{1}{\sigma_e + \sigma_c} (\sigma_c \Lambda_o^{\mathbb{T},+} g^+ - \sigma_e \Lambda_o^{\mathbb{T},-} g^-).$$

We diagonalise the above operators and conclude by means of the above relations. All the above problems associated to the operators $\Lambda_o^{\mathbb{T},\pm}$, $\Lambda_c^\mathbb{T}$ and $\Lambda_e^\mathbb{T}$ can be solved similarly due the symmetry of the domain. In all cases, we apply separation of variables, that is, we assume that solution has the following ansatz $v(x, y, z) = h(z)\omega_\mathbf{k}(x, y)$ (like for the sphere, $f = \omega_\mathbf{k}$). This results in the following ordinary differential equation (ODE)

$$\frac{d^2 h_{\mathbf{k}}}{dz^2} - \left(\frac{2\pi}{L}\right)^2 |\mathbf{k}|^2 h_{\mathbf{k}} = 0,$$
whose solution is of the form

$$h_k(z) = \begin{cases} A_k \sinh\left(\frac{2\pi}{L}z|\mathbf{k}|\right) + B_k \cosh\left(\frac{2\pi}{L}z|\mathbf{k}|\right), & \text{if } \mathbf{k} \neq 0, \\ A_0 z + B_0, & \text{if } \mathbf{k} = 0. \end{cases}$$

We determine the coefficients $A_k$ and $B_k$ in each case according to the definition of the operators. More precisely,

* In the case of the operator $\Lambda_e^T$, one has the conditions $h_k(0) = 1$ and $h_k(H) = 0$.

* In the case of the operator $\Lambda_c^T$, one has the conditions $h_k(0) = 1$ and $h_k(-H) = 0$.

* In the case of the operator $\Lambda_0^{T,+}$, one has the conditions $h_k(0) = 0$ and $h_k(H) = 1$.

* In the case of the operator $\Lambda_0^{T,-}$, one has the conditions $h_k(0) = 0$ and $h_k(-H) = 1$.

The simple calculations (left to the reader) lead to the expressions given by the lemma. $\square$

We use this setting to approximate what locally happens near the poles of the cell when we zoom in.

### 3.2.3 Spherical vs flat membranes

**Definition 15** (Time constant and (dis-)charging times the transmembrane voltage). We define the time constant of the membrane as

$$\tau^{\Lambda} = \frac{C_m(0)}{S_m(0) + \min_{\lambda \in \mathcal{S}(\Lambda) \setminus \{0\}}(\lambda)}, \quad (16)$$

where $\mathcal{S}(\Lambda)$ is the spectrum of $\Lambda$.

The charging (resp. discharging) time of the TMV subjected to a constant unidirectional electric field of magnitude $E = Ez$, that is $g^+ = EH = -g^-$, is defined as $5\tau^\Lambda$, corresponding in the linear regime to 99% of the complete charge (resp. discharge) of the membrane.

Simple calculations lead to the following property.

**Property 16** (Time constants of the TMV for spherical vs flat membranes). Using the expressions of the eigenvalues of the operators $\Lambda^S$ and $\Lambda^T$ given in Lemma 12 and Lemma 14 respectively, one gets

* In the case of a spherical membrane, one has

$$\min_{\lambda \in \mathcal{S}(\Lambda) \setminus \{0\}} (\lambda) = \lambda_1^S = \frac{1}{\frac{2R_0}{\sigma_c} + \frac{R_0 (R_1/R_0)^3 - 1}{1 + 2(R_1/R_0)^3}} \sim_{R_1 \to +\infty} \frac{\sigma_c}{R_0} \frac{2}{1 + 2\frac{\sigma_c}{\sigma_c + 2\sigma_c}}.$$

Therefore the time constant $\tau^S$ of a spherical membrane in $\mathbb{R}^3$ is given by

$$\tau^S \sim_{R_1 \to +\infty} \frac{C_m(0)}{S_m(0) + \frac{1}{R_0} \frac{2\sigma_c \sigma_c}{\sigma_c + 2\sigma_c}}.$$
• In the case of a flat torus,

$$ \tau^{\mathbb{T}} = \frac{C_m(0)}{S_m(0) + \frac{1}{H} \frac{\sigma_c \sigma_e}{\sigma_c + \sigma_e}}. $$

**Remark 17.** In order to consider the periodic flat membrane as a zoomed-in flat patch of the spherical cell model, it is thus necessary to match their time constants accordingly. This means that to compare both settings, we need to adjust the height H of the box as

$$ H = \frac{R_0}{2} \left( \frac{\sigma_c + 2\sigma_e}{\sigma_c + \sigma_e} \right). \qquad (17) $$

In molecular dynamics simulations of membrane electroporation, the dimension of the simulation boxes are a few nanometers thick (see for instance papers by Tieleman or Mounir [36, 37, 6, 35]), while the cell radius is about 5 to 10 µs. According to the above estimation, the corresponding time constant $\tau_{\text{MD}}^{\mathbb{T}}$ of the TMV in molecular dynamics simulation is thus about 3 order of magnitude smaller than the time constant $\tau_{\text{vesicle}}^{\mathbb{S}}$ of the TMV for spherical bilipid membrane (also call vesicle)

$$ \tau_{\text{MD}}^{\mathbb{T}} \sim 10^{-3} \tau_{\text{vesicle}}^{\mathbb{S}}. $$

Based on this observation, and due to the coupling between water content and TMV, we believe that this dramatic difference on the time constant of the TMV prevents molecular dynamics simulations to give quantitative information on the membrane electroporation, even though they provide interesting qualitative description of the phenomenon.

The following property provides a bound on the Fourier coefficients of the steady TMV in a membrane containing water.

**Property 18** (Steady TMV in porated membrane). Let $\bar{\phi}$ be a smooth phase order parameter constant in time on $\Gamma = (\mathbb{R}/(L\mathbb{Z}))^2$. Let $(p, q) \in (0, 1] \times [0, 1)$ be the proportion of pure water (resp. pure lipid) phases inside the flat membrane, and let $s_0$ be the mean of $S_m(\bar{\phi})$:

$$ p := \frac{1}{|\Gamma|} \int_{\{\bar{\phi}=1\}} dx, \quad q := \frac{1}{|\Gamma|} \int_{\{\bar{\phi}=0\}} dx, \quad \bar{S}_m = \frac{1}{|\Gamma|} \int_{\Gamma} S_m(\bar{\phi}) dx. $$

Denote by $\epsilon$ the term $\epsilon = 1 - p - q \in (0, 1)$. Let $v^{\infty}$ be the stationary solution of equation (14) with a constant source term $g^{+} = EH = -g^{-}$:

$$ (S_m(\bar{\phi}) + \Lambda)v^{\infty} = \frac{\sigma_e \sigma_e}{\sigma_c + \sigma_e} E. \qquad (18) $$

Assume $v^{\infty}$ that can be expanded in Fourier by

$$ v^{\infty}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^2} \xi_{\mathbf{k}} \exp \left( \frac{2i\pi}{L} \mathbf{k} \cdot \mathbf{x} \right), \quad \text{for a.e. } \mathbf{x} \in \Gamma. $$

Then one has the following estimate for $\mathbf{k} \in \mathbb{Z}^2$:

$$ |\xi_{\mathbf{k}}| \le \frac{\sqrt{p(1-p)(S_m(1) - S_m(0))^2 + \epsilon S_m^2(1)}}{\bar{S}_m + \frac{2\pi|\mathbf{k}|}{L \tanh(2\pi|\mathbf{k}|/H/L)} \frac{\sigma_e \sigma_e}{\sigma_e + \sigma_c}} \|v^{\infty}\|_{L^2, d\mu}, \qquad (19) $$

where $d\mu := dx/|\Gamma|$ is the probability measure so that

$$ \forall \psi \in L^2(\Gamma), \quad \|\psi\|_{L^2, d\mu} := \left( \frac{1}{|\Gamma|} \int_{\Gamma} |\psi(x)|^2 dx \right)^{1/2}. $$
**Remark 19** (Influence of the patch size). Letting L go to zero in (19), all the other parameter being fixed, shows that the Fourier coefficients $\xi_k$ the decreases linearly with L. Therefore small patches tend to flatten the TMV artificially even if we keep the same proportion of water in the membrane. This heuristic calculation suggests that large patches of flat membrane are needed so that flat torus can be seen as a zoom of the spherical setting.

*Proof.* The squared $L^2$ distance between $s_0$ and $S_m$ for the probability measure is bounded by

$$ \begin{aligned} \|S_m(\bar{\phi}) - \bar{S}_m\|_{L^2, d\mu}^2 &= \|S_m(\bar{\phi})\|_{L^2, d\mu}^2 - \bar{S}_m^2 \\ &\le S_m^2(1)p + S_m^2(0)q + S_m^2(1)\epsilon - \bar{S}_m^2 \\ &\le (S_m(1) - S_m(0))^2 p(1-p) + \epsilon S_m^2(1). \end{aligned} \quad (20) $$

Applying the operator $(\bar{S}_m + \Lambda)^{-1}$ to the static equation (18) leads to

$$ (\bar{S}_m + \Lambda)^{-1}(S_m(\bar{\phi}) - \bar{S}_m)v^\infty + v^\infty = \frac{2\sigma_c\sigma_e E}{\sigma_e + \sigma_c \bar{S}_m + \frac{1}{H} \frac{\sigma_c\sigma_e}{\sigma_c + \sigma_e}}. $$

Taking the $L^2$ dot-product for the probability measure $d\mu$ with $\exp(\frac{2i\pi}{L}\mathbf{k} \cdot \mathbf{x})$ and using (20) ends the proof. $\square$

# 4 The nonlinear coupled system and its stability analysis

Putting the TMV equation and the order parameter equation results in the following joint problem on $\Gamma$, we obtain

$$ \partial_t \phi - D_0 \Delta_{\Gamma} \phi = -\alpha W'(\phi) + \frac{\alpha}{2} C'_m(\phi) v^2, \quad \forall x \in \Gamma, \forall t > 0, \quad (21a) $$

$$ C_m(\phi)\partial_t v + (S_m(\phi) + \Lambda)v = G, \quad \forall x \in \Gamma, \forall t > 0, \quad (21b) $$

$$ v|_{t=0} = v_{\diamond}, \quad \phi|_{t=0} = \phi_{\diamond}, \text{ on } \Gamma. \quad (21c) $$

Let $(\bar{\phi}, \bar{v})$ be two constant solutions to (21) such that $\bar{\phi} < 0.5 - 3a_2/8a_1$. For small perturbations $\mu_{\diamond}, \varepsilon_{\diamond}$ of these initial conditions, we assume there exists a unique solution $(t, x) \mapsto (\phi(t, x), v(t, x))$. Let $\mu(t, x) := v(t, x) - \bar{v}$ and $\varepsilon(t, x) := \phi(t, x) - \bar{\phi}$, and set $\mu(0) = \mu_{\diamond}$ and $\varepsilon(0) = \varepsilon_{\diamond}$. Inserting these expressions into (21) and linearising results in the following system

$$ \partial_t \begin{bmatrix} \varepsilon \\ \mu \end{bmatrix} = \begin{bmatrix} D_0 \Delta_{\Gamma} - \alpha W''_m(\bar{\phi}) + \frac{\alpha}{2} C''_m \bar{v}^2 & \alpha C'_m(\bar{\phi}) \bar{v} \\ - \frac{S'_m(\bar{\phi}) \bar{v}}{C_m(\bar{\phi})} & - \frac{S_m(\bar{\phi}) + \Lambda}{C_m(\bar{\phi})} \end{bmatrix} \begin{bmatrix} \varepsilon \\ \mu \end{bmatrix}. \quad (22) $$

In both the case of $\Gamma = \mathbb{S}$ or $\Gamma = (\mathbb{R}/(L\mathbb{Z}))^2$ the operators $\Delta_{\Gamma}$ and $\Lambda$ are diagonalisable in the same $L^2$ basis. We can thus further simplify the problem by looking at the decomposition of $\varepsilon$ and $\delta$ in its Fourier modes (or spherical harmonics).

**Remark 20.** Due to $(\bar{\phi}, \bar{v})$ being constant, all of the operators in each entry of the square matrix in (22) commute.
We denote by $\varepsilon(t) = \sum \varepsilon_n(t)\omega_n$ and $\mu(t) = \sum \mu_n(t)\omega_n$ the Fourier (or harmonic) decomposition of the solution. Then each splitting the above linearized problem along each frequency $(\varepsilon_n, \mu_n)$ results in the following

$$ \partial_t \begin{bmatrix} \varepsilon_n \\ \mu_n \end{bmatrix} = \underbrace{\begin{bmatrix} -D_0 \lambda_n^\Delta - \alpha W_m''(\bar{\phi}) + \frac{\alpha}{2} C_m''(\bar{\phi}) \bar{v}^2 & \alpha C_m'(\bar{\phi}) \bar{v} \\ -\frac{S_m'(\bar{\phi}) \bar{v}}{C_m(\bar{\phi})} & -\frac{S_m(\bar{\phi}) + \lambda_n^\Lambda}{C_m(\bar{\phi})} \end{bmatrix}}_{:=A_n} \begin{bmatrix} \varepsilon_n \\ \mu_n \end{bmatrix} $$

where $\lambda_n^\Delta$ and $\lambda_n^\Lambda$ are the n-th eigenvalues of $-\Delta$ and $\Lambda$ respectively.

Linear instability occurs if there exists an integer $n$ (associated to non-constant eigenfunctions) such that $A_n$ has an eigenvalue with positive real part. This is the case if and only if $\det(A_n) < 0$ or $\mathrm{Tr}(A_n) > 0$ for some $n$. We assume that the conductivity of the membrane does not change much for $\bar{\phi} \in [0, 0.4]$ (see Section 5.4.1), we expect that $|S_m'(\bar{\phi})| \ll 1$. Therefore, the easiest condition to verify is $\det(A_n) < 0$ (for some $n$). In fact, we get exactly

$$ W_m''(\bar{\phi}) + \frac{D_0}{\alpha} \lambda_n^\Delta + \frac{S_m'(\bar{\phi}) C_m'(\bar{\phi})}{S_m(\bar{\phi}) + \lambda_n^\Lambda} \bar{v}^2 < \frac{C_m''(\bar{\phi})}{2} \bar{v}^2. $$

If $S_m'(\bar{\phi})$ is small enough for $\bar{\phi} \in [0, 0.4]$, we simplify this constraint by

$$ C_m''(\bar{\phi}) > \frac{2}{\bar{v}^2} W_m''(\bar{\phi}) + \frac{2}{\bar{v}^2} \frac{D_0}{\alpha} \lambda_n^\Delta + \rho, \quad \text{in units } [\mathrm{F.m}^{-2}], $$

for some (small) $\rho > 0$. Using the expression of $W_m$ given in (5), we infer

$$ C_m''(\bar{\phi}) > \frac{8a_1}{\bar{v}^2} \left( \left( 6\bar{\phi}^2 - \left(6 - \frac{3a_2}{a_1}\right)\bar{\phi} + \left(1 - \frac{3a_2}{2a_1}\right) \right) + \frac{D_0}{4a_1\alpha} \lambda_n^\Delta \right) + \rho. $$

In practice $a_2 \ll a_1$ and $\frac{D_0}{4a_1\alpha} \sim \delta^h$ which is the membrane thickness (see Section 5) this suggests the following simple sufficient condition allowing to study the effect of each parameter in our model

$$ C_m''(\bar{\phi}) > \frac{8a_1}{\bar{v}^2} \left( 6\bar{\phi}^2 - 6\bar{\phi} + 1 + (\delta^h)^2 \lambda_n^\Delta \right) + \rho. \quad (23) $$

**Remark 21.** Dividing the above inequality (23) by $C_m(\bar{\phi})$ shows that the emergence of instabilities is driven by $\frac{4a_1}{2C_m(\bar{\phi})\bar{v}^2}$, which is the ratio of the barrier energy of the double-well potential $W_m$ and the electrostatic energy of the membrane. In addition the parameter $\alpha$ impacts on the size of the initial instabilities created. This can be seen through the size of the highest eigenmodes (highest in the sense of the value of $|\vec{n}|$) that are linearly unstable.

* In the flat periodic case the eigenvalues of the surface Laplace-Beltrami read

$$ \lambda_{\vec{n}}^{\Delta} = \frac{4\pi^2|\vec{n}|^2}{L^2}, \quad \text{if } \Gamma = (\mathbb{R}/(L\mathbb{Z}))^2, $$

We take the largest eigenvalue $\lambda_{\vec{n}}$ such that $(\delta^h)^2 \lambda_{\vec{n}}^\Delta \sim 1$, which results in the following

$$ \frac{L}{|\vec{n}|} \sim 2\pi\delta^h. \quad (24) $$
As in the periodic case the $\vec{n}$-th eigenfunction is of wavelength $L/|\vec{n}|$, therefore linear instabilities reach eigenmodes of the lengthscale of the membrane thickness. This supports the idea that the size of pores created by the influence of the electric field are expected to initially be of the length scale of the membrane thickness.

• In the case of a spherical cell we get a similar result. Here the Laplace-Beltrami eigenvalues read

$$ \lambda_n^\Delta = \frac{n(n+1)}{R_0^2}, \quad \text{if } \Gamma = \mathbb{S}, \text{ the sphere of radius } R_0. $$

Here, the notion of wavelength of an eigenfunction is not directly applicable, however qualitatively, for higher values of n the eigenfunctions associated to $\lambda_n^\Delta$ oscillate more rapidly. Another way to see this is the following, for a given eigenvalue $\lambda_n^\Delta$, there is a spherical harmonic which divides the sphere in $2\lfloor n/2 \rfloor(1 + \lfloor n/2 \rfloor)$ regions with alternating signs (also known as nodal domains), see [23] for more information. Again, if we take the largest value of n such that $(\delta^h)^2\lambda_n^\Delta \sim 1$, then we get

$$ \frac{R_0^2}{n(n+1)} \sim (\delta^h)^2 \\ \frac{|\Gamma|}{4\pi n(n+1)} \sim (\delta^h)^2. $$

Therefore, we can see that there is at least one spherical harmonic with a nodal domain of size comparable to $(\delta^h)^2$ which is linearly unstable. We conclude similarly as in the flat membrane case, and we obtain a similar estimation of the initial size of pores.

# 5 Choice of membrane parameters

The aim of this section is to obtain reasonable estimates for the parameters in our model. Indeed, some of these parameters are not physiological, meaning that they cannot be directly measured. Therefore, we would like to infer a priori value ranges allowing realistic simulation. To this end, we are going to link our free-energy (2) to the energy model proposed by Chizmadzhev et al., from which the models of Krassowska, Weaver et al. are derived [20, 38, 40, 21].

## 5.1 Determination of the double-well potential coefficients

In the theory proposed by Abidor, Chizmadzhev and Weaver [1, 9, 39], the water pore is assumed to be a cylinder embedded in a sea of lipid. The introduced the notion of hydrophobic pores, defined as the pores with radius smaller that a pore parameter $r^*$. It was introduced to account for the fluctuation of lipids which does not generate water pathways in the membrane. Pores of radius above $r^*$ are defined as hydrophilic (see Figure 3) and enables the water to pass through the membrane.

Energy functional (2) (when $v = 0$) can be thought of as a generalisation of the standard energy model for hydrophilic pores (25). Interestingly, the use of reaction-diffusion on the water content enables us to describe the interaction between neighboring pores, which was not possible in the previous approach, since the pores were considered somehow isolated from each other. In order to
<figure><img src="image_4.png" /><figcaption>Figure 3: Energy needed to create a pore of radius $r$ according to Weaver and Chizmadzhev approach. A pore is considered to be hydrophilic if its radius $r$ is bigger than $r^*$, otherwise it is considered to be hydrophobic pore (see [39, 30]).</figcaption></figure>

compare the both approaches, it is instructive to consider the physical units of each term of the membrane free-energy (2).

As $\phi$ is dimensionless then the unit for $\kappa$ is the energy unit [J] and the coefficients $a_1$ and $a_2$ have surface tension units [J.m$^{-2}$]. We thus compare these terms with the linear and surface tension terms from the hydrophilic pore energy in [39]. Let $r^* > 0$ be the radius at which hydrophobic pores are more energetically favorable than hydrophilic pores (see [39, 30] for a description of $r^*$, in the literature its value is estimated to be $r^* := 0.5[\text{nm}]$), then for any radius $r > r^* > 0$ the energy needed to create a pore (see Figure 3) of that radius is

$$w(r) = 2\pi\gamma r - \pi\sigma r^2 + \frac{C}{r^4}, \quad \forall r > 0, \quad (\text{hydrophilic pore energy}), \quad (25)$$

where $\gamma$ is the linear tension of the edge of a pore and $\sigma$ the surface tension of the membrane. The last term represents the steric repulsion of the lipid heads in the edge of a pore. In practice it prevents (hydrophilic) pores from being too small which would be too energetically costly. In order to make the link with our energy functional more concrete, we are going to consider the case of an "ideal" smooth pore. This allows us to compare each term in (25) with the terms in our energy functional. Let $\phi_p$ be a smooth circular pore of radius $r_0 > 0$ with a small interface (where $\phi \in (0, 1)$) of size $\delta^h \ll r_0$ (see Figure 4a). The energy difference between an intact membrane $\phi \equiv 0$ and a membrane with such a pore $\phi_p$ is given by

$$\begin{aligned} \mathcal{E}(\phi_p, 0) - \mathcal{E}(0, 0) &= \frac{\kappa}{2} \int_{\Gamma} |\nabla \phi_p|^2 + \int_{\text{interface}} (W_m(\phi_p) - W_m(0)) dx - \int_{\text{interior}} W_m(0) dx, \\ &= \frac{\kappa}{2} \int_{\Gamma} |\nabla \phi_p|^2 + \pi((r + \delta^h)^2 - r^2)(\bar{W}_m - \frac{a_2}{2}) - \pi r^2 \frac{a_2}{2}, \end{aligned} \quad (26)$$

where $\bar{W}_m$ is the mean value of $x \in \Gamma \mapsto W_m(\phi_p(x))$ inside the pore interface. As $\phi_p$ takes values in $[0, 1]$ for all $x \in \Gamma$, we deduce that $\bar{W}_m \in [0, W_m(1/2 + 3a_2/2a_1)]$. Let $\gamma_0 := \bar{W}_m - a_2/2$, then
<figure><img src="image_5.png" /><figcaption>Figure 4: Ideal smooth circular pore $\phi_p$ of radius $r_0 > 0$ in a flat membrane $\Gamma$ centered at the origin. The rotational symmetry of the setting allows us to describe the pore as a function $r \mapsto \phi_p(r)$ where $r$ is the distance to the origin. The thickness of the lipid bilayer is denoted by $h > 0$. The interface is of size $\delta^h > 0$ such that $\delta^h \ll r_0$.</figcaption></figure>

(26) can be rewritten as

$$ \mathcal{E}(\phi_p, 0) - \mathcal{E}(0, 0) = 2\pi r_0 (\delta^h \gamma_0) - \pi r_0^2 \frac{a_2}{2} + \left( \frac{\kappa}{2} \|\nabla \phi_p\|_{L^2(\Gamma)}^2 + \pi (\delta^h)^2 \gamma_0 \right). \quad (27) $$

Looking at (27), the parallel between our energy functional and (25) becomes clear. We already know that $a_1$ and $a_2$ represent surface tensions and we therefore deduce that $\delta^h \gamma_0$ corresponds to the linear tension of the pore's edge $\gamma$ from (25), and similarly for the term $a_2/2$ and $\sigma$ from (25). This results in the following system of equations

$$ \bar{W}_m - \frac{a_2}{2} = \frac{\gamma}{\delta^h}, \quad \text{and} \quad a_2 = 2\sigma, $$

where $a_1$ and $a_2$ are the unknowns. The size $\delta^h$ of the interface for a typical pore can be roughly estimated. The size of $\delta^h$ intuitively should be comparable with half the membrane thickness $h/2$ (see Figure 4b). In order to determine the value of $a_1$ expressed in terms of the model proposed in [39], we need to invert the following nonlinear equation

$$ \bar{W}_m = \frac{\gamma}{\delta^h} + \sigma, $$

where $\bar{W}_m$ depends on $a_1$, $a_2$ and the shape of $\phi_p$ inside the interface.

As a first approximation, we estimate the value of $\bar{W}_m$ by assuming that $r \mapsto \phi_p(r)$ changes almost linearly inside the interface (meaning that $d\phi_p/dr(r) \sim 1/\delta^h$ inside the interface). This results in the following rough estimate

$$ \bar{W}_m \sim \int_0^1 a_1 \phi^2 (1-\phi)^2 + a_2 (\phi + 1/2)(\phi - 1)^2 d\phi = \frac{a_1}{30} + \frac{a_2}{4}. $$

Replacing this expression in the above equations results in the following values

$$ a_1 = 30 \frac{\gamma}{\delta^h} + 15\sigma, \quad \text{and} \quad a_2 = 2\sigma. $$

The last term in (27) will play a similar role to the steric repulsion term when a pore is too small. In the next paragraph, we continue to estimate all the values of the model before presenting the table summarizing the values used.
## 5.2 Parameters influencing membrane dynamics under null TMV

We can now estimate the values of $a_1$ and $a_2$, in addition to the kinetic coefficient $\alpha$. These values depend only on the cell membrane and so we consider $v \equiv 0$ throughout this section. To estimate the values of $a_1$ and $a_2$, we use the values of linear and surface tension from [20]. To estimate $\alpha$, we can roughly approximate it by using Property 7 in Section 2.3. We take the asymptotic relation bounding the size of the anti-phase interface ($\delta_n^h$ instead of $\delta^h$) from that property to be an equality and we obtain the following rough estimate

$$ \alpha = \frac{D_0}{4a_1(\delta^h)^2} \quad [\text{m}^2\text{J}^{-1}\text{s}^{-1}]. $$

In practice, given estimates for $a_1$, $a_2$, $D_0$ and $\delta^h$, we estimate $\alpha$ via the above formula and then verify the conditions from Property 7.

## 5.3 Summary of parameters

We summarize the choice of parameters used in our simulations in Table 1 and 2. The conductivity of an intact membrane (lipid conductivity $\sigma_l$) is estimated by looking at its characteristic time of charge

$$ \tau_m = \frac{C_m(0)}{S_m(0) + \frac{2\sigma_e\sigma_c}{R_0(\sigma_c + 2\sigma_e)}}, \quad \text{for an intact spherical membrane, with } C_m(0) = \varepsilon_0\varepsilon_l/h, $$

where $h$ is the membrane thickness, $\varepsilon_0\varepsilon_l$ the lipid permittivity, $\sigma_c$ and $\sigma_e$ the conductivities of the cell interior and exterior respectively, and $R_0$ is the cell radius. According to experimental observations (we refer the reader to [32, 33] and references therein), the value of $\tau_m$ ranges between 0.1 µs and 1 µs.

<table><caption>Table 1: Physical parameters of the cell electroporation model.</caption><thead><tr><th></th><th>Value</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td>$D_0$</td><td>2 × 10<sup>-12</sup> m<sup>2</sup>s<sup>-1</sup></td><td>lateral diffusion of lipids in the cell membrane</td><td>[25]</td></tr><tr><td>$\sigma$</td><td>10<sup>-6</sup> J.m<sup>-2</sup></td><td>surface tension of the lipid-bilayer</td><td>[20]</td></tr><tr><td>$\gamma$</td><td>1.8 × 10<sup>-11</sup> J.m<sup>-1</sup></td><td>pore edge energy</td><td>[20]</td></tr><tr><td>$\epsilon_0$</td><td>8.85 × 10<sup>-12</sup> F.m<sup>-1</sup></td><td>vacuum permittivity</td><td>[28]</td></tr><tr><td>$\epsilon_w$</td><td>80</td><td>relative water permittivity</td><td>[39]</td></tr><tr><td>$\epsilon_l$</td><td>2</td><td>relative lipid permittivity</td><td>[39]</td></tr><tr><td>$\sigma_c$</td><td>1 S.m<sup>-1</sup></td><td>interior medium conductivity</td><td>[33, 32]</td></tr><tr><td>$\sigma_e$</td><td>1.5 S.m<sup>-1</sup></td><td>exterior medium conductivity</td><td>[32, 33]</td></tr><tr><td>$\sigma_l$</td><td>3 × 10<sup>-5</sup> S.m<sup>-1</sup></td><td>lipid membrane conductivity</td><td>[11, 20]</td></tr><tr><td>$\sigma_w$</td><td>1 S.m<sup>-1</sup></td><td>pore conductivity</td><td>[11, 20]</td></tr><tr><td>$h$</td><td>5 × 10<sup>-9</sup> m</td><td>membrane thickness</td><td>[20]</td></tr></tbody></table>
<table><caption>Table 2: Estimated parameters specific to our model. We consider them up to order of magnitude, as they are taken to be rough approximations deduced from our analysis. That is why they may not verify the exact formulas used to defined them.</caption><thead><tr><th></th><th>Value</th><th>Description</th><th>Section</th></tr></thead><tbody><tr><td>δ<sup>h</sup></td><td>10<sup>-9</sup> m</td><td>pore edge size</td><td>(Section 5.1)</td></tr><tr><td>α</td><td>9.26×10<sup>6</sup> m<sup>2</sup>.J<sup>-1</sup>.s<sup>-1</sup></td><td>kinetic coefficient of our model</td><td>(Section 5.2)</td></tr><tr><td>a<sub>1</sub></td><td>5.44 × 10<sup>-1</sup> J.m<sup>-2</sup></td><td>measures energy barrier between the two stable states</td><td>(Section 5.1)</td></tr><tr><td>a<sub>2</sub></td><td>2.0 × 10<sup>-6</sup> J.m<sup>-2</sup></td><td>measures the surface tension of the membrane</td><td>(Section 5.1)</td></tr></tbody></table>

**Remark 22.** It is worth noting that with these choices of parameters, the coefficient $\sqrt{\frac{4a_1\alpha L^2}{D_0}}$ is large and $\frac{3a_2}{4a_1}$ small. Therefore, the asymptotic results, presented in Property 7, hold.

## 5.4 Choice of the conductance and capacitance model

In this section, we look at one behaviour of our model according to the model choices for $S_m$ and $C_m$. Previous models require the values of these functions for the lipid ($\phi = 0$) state or the pore state ($\phi = 1$). In this continuous setting, we have a modeling choice to make to extend these properties to "mixed" states between lipid and pore (for example $\phi = 1/2$). For both $S_m$ and $C_m$, the main approach will be to interpolate between the values of the function for the stable states. Of course, whatever the choice of interpolation, our model is only concerned by the values in a neighbourhood of $\phi \in [0, 1]$. We can for example consider linear interpolations and still consider $S_m$ and $C_m$ to be bounded if we just apply appropriate cutoff functions outside the region of interest.

### 5.4.1 Membrane conductance model

The choice of $S_m$ directly impacts on the amplitude of the transmembrane voltage $v$. In order to interpolate between the conductivity of the membrane (lipid conductivity) and the conductivity of a pore, we consider the following sigmoid function

$$S_m(\psi) = \frac{\sigma_m(\psi)}{h}, \quad \text{where} \quad \sigma_m(\psi) = \frac{1 + \tanh(k_0(\psi - \phi_{\text{th}}))}{2}(\sigma_w - \sigma_l) + \sigma_l, \quad \forall \psi \in \mathbb{R},$$

where $\sigma_l$ and $\sigma_w$ are the lipid and water conductivities respectively. The factor $k_0 > 0$ determines the steepness of the transition between these two conductivities and the value of $\phi_{\text{th}} \in [0, 1]$ is the water content threshold of the transition. The effect of electroporation on the cell membrane tends to suddenly and dramatically increases its conductivity. As $\phi$ evolves continuously in time and we do not expect a gradual change in membrane conductivity during its evolution, this translates into a steep transition for $S_m$. In other words $k_0$ cannot be too small if we wish to model this effect correctly. As for an appropriate value for $\phi_{\text{th}}$, this will depend on the value of $k_0$ to some extent, although it is intuitive to impose $\phi_{\text{th}} \sim 1/2$ as it is an instability mid point between the two stable states of the double well energy potential $W_m$. In the following we set $k_0 = 100$.

## 5.5 Influence of the capacitance model

In contrast to $S_m$, $C_m$ directly affects the order parameter $\phi$ as it is directly responsible (along with the TMV) for pore formation, which is why we must also consider the membrane dynamics in this analysis.
Looking at the Allen-Cahn equation (21a), to promote water entering the membrane, we need $C'_m(\phi) > 0$ (at least for $\phi \in [0, 1/2]$). Furthermore, we do not expect that the membrane order parameter $\phi$ uniformly increases as the TMV increases. The model of Looyenga given in (3) ensures that so that constant solutions to the joint problem (21) are (linearly) unstable and that $C'_m(\phi) > 0$.

# 6 Numerical simulations on a patch membrane

We present in this section the numerical simulations in the case of a flat periodic membrane, to take advantage of the Fast Fourier Transform (FFT), that enables to perform fast simulations [16].

## 6.1 Numerical scheme

In this section, we present the numerical schemes used to solve the coupled problem.

### 6.1.1 Scheme to solve the Allen-Cahn equation

We start with the scheme for the adimensionalised equation of the evolution of the order parameter. Let $L > 0$ be the characteristic length of $\Gamma$ and $\tau > 0$ the characteristic time. Then rescaling the equation on $\mathbb{T}^2 = \mathbb{R}^2/\mathbb{Z}^2$ leads to

$$ \begin{aligned} \partial_t \phi - \frac{D_0 \tau}{L^2} \Delta_{\mathbb{T}}^2 \phi &= -4\tau \alpha a_1 \phi(\phi - 1)(\phi - 1/2 + 3a_2/4a_1) + \frac{\tau \alpha}{2} C'_m(\phi) v^2, \quad \text{in } ]0, +\infty[ \times \tilde{\Gamma}, \\ \phi(t=0) &= \phi^0, \end{aligned} $$

where $x \in \mathbb{T}^2 \mapsto v(x)$ is constant in time. We solve this equation by using an operator splitting scheme (for an example see [24]). The surface $\mathbb{T}^2$ is discretised with a Cartesian grid.

We denote by $\Phi^n$ the spatial and temporal discretisation at time $t^n$ of $\phi(t^n, \cdot)$. At time $t^{n+1} := t^n + \delta t$ the numerical approximation $\Phi^{n+1}$ of $\phi(t^{n+1}, \cdot)$ is obtained by numerically solving the following ordinary differential equation (ODE) at each node of the square grid $(x_p, y_q)_{\{p,q\}} \in \tilde{\Gamma}$

$$ \left\{ \begin{array}{l} \frac{d\phi}{dt} = -4\tau\alpha a_1\phi(\phi-1)(\phi-1/2+3a_2/4a_1) + \frac{\tau\alpha}{2}C'_m(\phi)v^2, \quad t \in (t^n, t^{n+1}), \\ \phi(t^n, (x_p, y_q)) = (e^{\delta t \Delta} \Phi^n)(x_p, y_q), \end{array} \right. $$

where $e^{\delta t \Delta}$ (which is the heat semigroup) can be explicitly calculated by diagonalising the Laplacian in Fourier domain. Concretely, if

$$ \phi(t^n, \vec{x}) = \sum_{\vec{k} \in \mathbb{Z}^2} c_{\vec{k}} e^{2i\pi \vec{k} \cdot \vec{x}}, \quad \forall \vec{x} \in \tilde{\Gamma}, \qquad (28) $$

then

$$ e^{\delta t \Delta} \phi(t^n, \vec{x}) = \sum_{\vec{k} \in \mathbb{Z}^2} (e^{-\delta t 4\pi^2 |\vec{k}|^2 \frac{\tau D_0}{L^2}} c_{\vec{k}}) e^{2i\pi \vec{k} \cdot \vec{x}}. $$

The above calculation can be obtained explicitly for $\Phi^n$ by means of the Fast Fourier Transform (FFT). Once $e^{\delta t \Delta} \Phi^n$ is calculated, Equation (28) is numerically solved by means of a Taylor method
of order two. This approach is easy to implement and particularly fast to solve. It can be summarized as follows

$$
\left\{
\begin{aligned}
\Phi^* &= \text{FFT}^{-1} \cdot \left( e^{-\left(\frac{\delta t}{4\pi^2|\vec{k}|^2 \frac{\tau D_0}{L^2}}\right)} \right)_{\vec{k}} \cdot \text{FFT}(\Phi^n), \\
\Phi^{n+1} &= \Phi^* + \delta t F(\Phi^*)(1 + \frac{\delta t}{2} F'(\Phi^*)),
\end{aligned}
\right.
$$

with

$$
F : y \mapsto -4\tau\alpha a_1 y(y-1)(y-1/2 + 3a_2/4a_1) + \frac{\tau\alpha}{2} C'_m(y)v^2. \quad (29)
$$

Applying a Strang Splitting scheme to the above fractional step method gives order 2 accuracy in
this setting [27].

6.1.2 Scheme to solve the transmembrane voltage

The equation on the TMV is also numerically solved by a method of operator splitting. Following
the usual adimensionalisation (rescaling by L in space, and τ in time) as before, we obtain

$$
\begin{cases}
C_m(\phi)\partial_t v + \tau(S_m(\phi) + \Lambda^{\mathbb{T}})v = \tau G, & \forall t > 0, \quad \forall x \in \mathbb{T}^2, \\
v(0) = v_{\diamond},
\end{cases}
$$

where the operator $\Lambda^{\mathbb{T}}$ is described by the same eigenvalues as in Section 3, but with eigenfunctions

$$
\omega_{(k_1, k_2)}(x, y) := e^{2i\pi(k_1 x + k_2 y)}, \quad \forall (k_1, k_2) \in \mathbb{Z}^2, \quad \forall x \in \mathbb{T}^2.
$$

Assume that the function $x \in \mathbb{T}^2 \mapsto \phi(x)$ is constant in time, as well as $G$ (corresponding to the flat part of an ideal pulse). Let $V^n$ denote the spatial and temporal discretisation of $v(t^n, \cdot)$, then to obtain the numerical approximation $V^{n+1}$ of $v(t^{n+1}, \cdot)$ we numerically solve the following scheme

$$
\frac{V^* - V^n}{\delta t/2} + \frac{\tau}{2} \frac{1}{C_m(\phi)} (\Lambda^{\mathbb{T}} - \lambda_0^{\mathbb{T}}) (V^* + V^n) = 0, \quad (30)
$$

$$
V^{**} = e^{-\tau \delta t \left( \frac{S_m(\phi) + \lambda_0}{C_m(\phi)} \right)} V^* + \frac{G}{S_m(\phi) + \lambda_0^{\mathbb{T}}} \left( 1 - e^{-\tau \delta t \left( \frac{S_m(\phi) + \lambda_0^{\mathbb{T}}}{C_m(\phi)} \right)} \right), \quad (31)
$$

$$
\frac{V^{n+1} - V^{**}}{\delta t/2} + \frac{\tau}{2} \frac{1}{C_m(\phi)} (\Lambda^{\mathbb{T}} - \lambda_0^{\mathbb{T}}) (V^{n+1} + V^{**}) = 0, \quad (32)
$$

where $\Lambda^{\mathbb{T}} - \lambda_0^{\mathbb{T}}$ is a nonnegative symmetric operator. We remark that (31) corresponds to the ordinary differential equation (ODE) part of the splitting, while (30) and (32) correspond to the unbounded part of the splitting (where we have applied a Strang splitting scheme as in the section above). In order to invert equations (30) and (32) we use a Conjugate Gradient (CG) method [17]. For example, in the case of (30) this is possible by rewriting it as

$AY = B$

where

$$
Y := \sqrt{C_m(\phi)} V^*, \quad A := \left( \mathrm{Id} + \frac{\tau \delta t}{2} \frac{1}{\sqrt{C_m(\phi)}} \left( \Lambda^{\mathbb{T}} - \lambda_0^{\mathbb{T}} \right) \frac{1}{\sqrt{C_m(\phi)}} \right), \text{ and}
$$

$$
B := \left( \mathrm{Id} - \frac{\tau \delta t}{2} \frac{1}{\sqrt{C_m(\phi)}} \left( \Lambda^{\mathbb{T}} - \lambda_0^{\mathbb{T}} \right) \frac{1}{\sqrt{C_m(\phi)}} \right) \sqrt{C_m(\phi)} V^n.
$$
As A is symmetric positive definite operator (matrix) we can apply CG to solve this equation. Lastly, $V^*$ is then recovered by dividing the solution Y by $\sqrt{C_m(\phi)}$.

**Remark 23.** We gain the additional advantage of being able to use FFT in the CG algorithm. This greatly speeds up the method as the matrix A is generally a full matrix due to $\phi$ not necessarily being constant in space.

### 6.1.3 Joint Problem

A simple scheme to solve the joint problem is then obtained by combining the last two solvers. Meaning that for a discrete numerical approximations $(\Phi^n, V^n)$ of $(\phi(t^n, \cdot), v(t^n, \cdot))$ we calculate the approximate solution $(\Phi^{n+1}, V^{n+1})$ at time $t^{n+1}$ by applying the above solvers but replacing $v$ by $V^n$ in Equation (29) and $C_m(\phi)$ by $C_m(\Phi^n)$ in (30 - 32). Now we are ready to show numerical results.

## 6.2 Numerical results

We simulated a flat biperiodic membrane of size $L = 200$ [nm], see Figure 2-Right, under the influence of a 12 ns square pulse of a uniform electric field of intensity $|E| = 3.2 \cdot 10^6$ [V.m$^{-1}$]. After the membrane discharges its TMV we continue to simulate it for another 10 $\mu$s. We pick as the initial condition

$$ \phi_{\diamond} : (x, y) \mapsto \mu(x, y) \qquad (33) $$

where $|\mu| < 1.5 \cdot 10^{-2}$ is a Gaussian random noise (which has been smoothed for numerical stability purposes).

During the initial charge time the membrane evolves more or less uniformly in space until about $t = 12$ [ns]. At this time, the electric field is turned off but the membrane remains charged. It then discharges as a capacitor, slowly over time, while the order parameter continues to evolve under the influence of the remaining TMV. At around $t = 27$ [ns] the conductive effect of the membrane enters into play. At this point, the membrane is still charged but now discharges mainly through the high conductivity spots that result from the affected order parameter ($\phi > 1/2$). These places can be seen more easily in Figure 5e in the distribution of the TMV on the membrane. However, as the membrane charge is small, the membrane discharges within nanoseconds of this phenomenon taking place. This sudden discharge of the TMV leads us to consider that electroporation has occurred. Although, a clear membrane configuration as the one described in Property 7 does not take place but until about 1 $\mu$s later. This is in accordance with the time estimate described in that property. The initial perturbation of the membrane grows (between $t = 0$ [ns] and $t = 30$ [ns]) resulting in a stronger perturbation in accordance with linear instability analysis done in Section 4.

In this setting the membrane is almost completely discharged by $t = 33$ [ns]. From that point on, the transmembrane voltage does not play a significant role and only the Allen-Cahn equation drives the membrane dynamic.

We can see in Figure 5 the evolution of the TMV shortly after the initial electric pulse. The lifetime of pores depends on their size. After the pulse, the pores remains open several tens microsecond as shown in Figure 6. It is worth noting that the pore closure is driven by the mean curvature, as stated for the Allen-Cahn equation by Bellettini [4]. This has to be linked to Kroeger et al. paper in which a curvature-driven pore closure is proposed [21]. That is why we see the pores becoming rounder as they shrink in size in Figure 6.
<figure><img src="image_6.png" /><figcaption>Figure 5: Simulations of System (21). The upper row corresponds to the membrane order parameter $\phi$ while the bottom row are the associated TMV.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Simulations of System (21). The evolution of the membrane order parameter $\phi$ for 10 µs after the end of the pulse delivery.</figcaption></figure>
**Remark 24.** With respect to computational efficiency, we report that the entire simulation took 3.5 hours in a computer armed with a 2.4 GHz Dual-Core Inter Core i5 processor. The first 35 [ns] of simulation took 3.25 hours and the rest of the time was used to simulate the next 10 [µs]. For the first 35 [ns] we used a time step of δt = 0.02 [ns] and for the rest of the simulation δt = 15 [ns]. We highlight the fact that our implementation has not been fully optimized nor parallelized.

## 7 Conclusion

We have proposed a new model of membrane permeabilisation based on the free-energy of the membrane. The problem consists in a surface reaction-diffusion of the water content of the membrane coupled to the nonlocal equation on the transmembrane voltage as described by the system equation (21). Thanks to the analysis of the eigenvalue of the nonlocal operator on the transmembrane voltage, we have compared the time constants in the spherical and the flat periodic settings. This fine analysis enables us to state that in order to have comparable time constants in both spherical and flat settings, it is necessary to have a flat periodic setting with a height of the order of the sphere radius as shown in Equation (17).

This observation explains why the time constants obtained in the molecular dynamic simulations (in which the height box of simulation is a few tens of nanometers) are most likely too short compared with cell membrane electroporation.

We have also shown simulations in a flat periodic settings, thanks to a well-designed splitting scheme combined with Fast Fourier Transforms. It is worth noting that solving the problem in a geometry close to the cell geometry is still challenging.

Indeed, a pore edge is somehow comparable in size to the width of the cell membrane $h \sim 10^{-9}m$ and the characteristic length of a cell membrane is comparable with the cell radius $R_0 \sim 10^{-5}m$. Then to capture the evolution of pores we need a discretised membrane with a minimum amount of nodes of around

$$ \frac{R_0^2}{h^2} \sim 10^8. $$

Additionally, calculating Dirichlet-to-Neumann operators is particularly computing resource consuming even in the case where we know how to diagonalise it. Therefore such simulations require specific numerical strategies for general shape of the membrane. We are however confident that a well designed use of spherical harmonics should enable the fast simulations of our electropermeabilisation model in a near future.

From the modeling view point the free-energy of the membrane could be complexified. Indeed, Breton and Mir have shown that the intense electric field generates an oxidation of the lipids. This means that in addition to the pore, a new species of oxidised lipid is created. The interfaces oxidised-nonoxidised lipids could increase the permeability of the membrane, and thus could explain the long duration of the permeabilised state of the membrane [7]. A possible extension of our model to take into account lipid oxidation consists in associating a new energy potential $W_{\text{oxi}}$ to oxidised phospholipids and then describe the state of the membrane by means of an additional order parameter $\rho$, which describes the degree of oxidation of the membrane. For example we could conceive of having an energy functional of the form

$$ \mathcal{E}(\phi, \rho, v) = \frac{\kappa}{2} \int_{\Gamma} \|\nabla\phi\|^2 + \int_{\Gamma} W_{\text{oxi}}(\phi) f_{\text{sig}}(\rho) + W(\phi)(1 - f_{\text{sig}}(\rho)) - \frac{1}{2} \int_{\Gamma} C_m(\phi, \rho)v^2 \quad (34) $$
where $f_{\text{sig}} : \mathbb{R} \to [0, 1]$ is a sigmoid function to be chosen to interpolate between the two states of the membrane (oxidised or not). Additionally, the order parameter $\rho$ would follow a reaction-diffusion equation on the membrane. This would capture the oxidation process of phospholipids and the lateral diffusion of the oxidised phospholipids in the membrane. This bare-bones idea for an extension of our model is heavily inspired from [22] where a more complete, but phenomenological, model of electroporation is presented.

Moreover, as presented in the stability analysis section, there is an important interplay between the membrane potential $W_m$ and the modeling choice of $\phi \mapsto C_m(\phi)$. Seeing as we chose the $W_m$ for the sake of simplicity, a more physically sensible model should be possible with a more appropriate choice of membrane potential and maybe a different capacitance model. This choice will have little effect on the pore evolution once the membrane is discharged, so we can expect the same pore sealing behavior. These perspectives will be addressed in forthcoming work.

## Acknowledgements

This research has been partly granted by the Plan Cancer Inserm project MECI (2022-2025, no 21CM119-00).

## Conflicts of Interest

The authors have no conflict of interest to declare.

## References

[1] I. Abidor, V. Arakelyan, L. Chernomordik, Y. Chizmadzhev, V. Pastushenko, and M. Tarasevich. Electric breakdown of bilayer lipid membranes: I. the main experimental facts and their qualitative discussion. *Journal of Electroanalytical Chemistry and Interfacial Electrochemistry*, 104, 1979.

[2] M. Alfaro, D. Hilhorst, and H. Matano. The singular limit of the Allen-Cahn equation and the FitzHugh-Nagumo system. *Journal of Differential Equations*, 245(2), 2008.

[3] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. *Acta Metallurgica*, 27(6), 1979.

[4] G. Bellettini. Lecture optnotes on mean curvature flow: Barriers and singular perturbations. Lecture optnotes (Scuola Normale Superiore). Scuola Normale Superiore, Pisa, Italy, 2014 edition, 2014.

[5] A. J. Bray. Theory of phase-ordering kinetics. *Adv. Phys.*, 51(2), 2002.

[6] M. Breton, L. Delemotte, A. Silve, L. M. Mir, and M. Tarek. Transport of siRNA through lipid membranes driven by nanosecond electric pulses: an experimental and computational study. *J. Am. Chem. Soc.*, 134(34), 2012.

[7] M. Breton and L. M. Mir. Investigation of the chemical mechanisms involved in the electropulsation of membranes at the molecular level. *Bioelectrochemistry*, 119, 2018.
[8] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. i. interfacial free energy. J. Chem. Phys., 28(2), 1958.

[9] Y. Chizmadzhev, F. Cohen, A. Shcherbakov, and J. Zimmerberg. Membrane mechanics can account for fusion pore dilation in stages. Biophysical Journal, 69(6), 1995.

[10] K. A. DeBruin and W. Krassowska. Electroporation and shock-induced transmembrane potential in a cardiac fiber during defibrillation strength shocks. Ann. Biomed. Eng., 26(4), 1998.

[11] K. A. DeBruin and W. Krassowska. Modeling electroporation in a single cell. i. effects of field strength and rest potential. Biophys. J., 77(3), 1999.

[12] R. Dimova and C. Marques, editors. The giant vesicle book the giant vesicle book. Productivity Press, New York, NY, 2019.

[13] L. C. Evans, H. M. Soner, and P. E. Souganidis. Phase transitions and generalized motion by mean curvature. Communications on Pure and Applied Mathematics, 45(9), 1992.

[14] A. Farina, Y. Sire, and E. Valdinoci. Stable solutions of elliptic equations on riemannian manifolds. J. Geom. Anal., 23(3), 2013.

[15] D. Fraggedakis, M. Mirzadeh, T. Zhou, and M. Z. Bazant. Dielectric breakdown by electric-field induced phase separation. Journal of The Electrochemical Society, 167(11), 2020.

[16] M. Frigo and S. G. Johnson. The design and implementation of fftw3. Proceedings of the IEEE, 93(2), 2005.

[17] M. R. Hestenes, E. Stiefel, et al. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6), 1952.

[18] O. Kavian, M. Leguèbe, C. Poignard, and L. Weynans. “Classical” Electropermeabilization Modeling at the Cell Scale. Journal of Mathematical Biology, 68, 2014.

[19] D. Kondepudi and I. Prigogine. Modern thermodynamics. John Wiley & Sons, Nashville, TN, 2 edition, 2014.

[20] W. Krassowska and P. D. Filev. Modeling electroporation in a single cell. Biophys. J., 92(2), 2007.

[21] J. H. Kroeger, D. Vernon, and M. Grant. Curvature-driven pore growth in charged membranes during charge-pulse and voltage-clamp experiments. Biophys. J., 96(3), 2009.

[22] M. Leguèbe, A. Silve, L. Mir, and C. Poignard. Conducting and permeable states of cell membrane submitted to high voltage pulses: Mathematical and numerical studies validated by the experiments. Journal of Theoretical Biology, 360, 2014.

[23] J. Leydold. On the number of nodal domains of spherical harmonics. Topology, 35, 1996.

[24] Y. Li, H. G. Lee, D. Jeong, and J. Kim. An unconditionally stable hybrid numerical method for solving the allen-cahn equation. Computers & Mathematics with Applications, 60(6), 2010.
[25] G. Lindblom and G. Oradd. Lipid lateral diffusion and membrane heterogeneity. *Biochimica et Biophysica Acta (BBA) - Biomembranes*, 1788(1), 2009.

[26] H. Looyenga. Dielectric constants of heterogeneous mixtures. *Physica*, 31(3), 1965.

[27] R. I. McLachlan and G. R. W. Quispel. Splitting methods. *Acta Numerica*, 11, 2002.

[28] P. J. Mohr, B. N. Taylor, and D. B. Newell. CODATA recommended values of the fundamental physical constants: 2006. *Rev. Mod. Phys.*, 80(2), 2008.

[29] C. Muller. Spherical Harmonics. Lecture optnotes in mathematics. Springer, Berlin, Germany, 1966 edition, 1966.

[30] J. C. Neu and W. Krassowska. Asymptotic model of electroporation. *Phys. Rev. E*, 59, 1999.

[31] A. Pazy. Semigroups of linear operators and applications to partial differential equations. Applied mathematical sciences. Springer, New York, NY, 1983 edition, 2012.

[32] A. Silve, I. Leray, M. Leguèbe, C. Poignard, and L. M. Mir. Cell membrane permeabilization by 12-ns electric pulses: Not a purely dielectric, but a charge-dependent phenomenon. *Bioelectrochemistry*, 106, 2015.

[33] A. Silve, I. Leray, C. Poignard, and L. M. Mir. Impact of external medium conductivity on cell membrane electropermeabilization by microsecond and nanosecond electric pulses. *Sci. Rep.*, 6(1), 2016.

[34] K. C. Smith, R. S. Son, T. Gowrishankar, and J. C. Weaver. Emergence of a large pore subpopulation during electroporating pulses. *Bioelectrochemistry*, 100, 2014.

[35] M. Tarek. Atomistic simulations of electroporation of model cell membranes. *Adv. Anat. Embryol. Cell Biol.*, 227, 2017.

[36] D. P. Tieleman. The molecular basis of electroporation. *BMC Biochem.*, 5, 2004.

[37] P. T. Vernier, Z. A. Levine, Y.-H. Wu, V. Joubert, M. J. Ziegler, L. M. Mir, and D. P. Tieleman. Electroporating fields target oxidatively damaged areas in the cell membrane. *PLOS ONE*, 4(11), 2009.

[38] J. Weaver. Electroporation of biological membranes from multicellular to nano scales. *IEEE Transactions on Dielectrics and Electrical Insulation*, 10(5), 2003.

[39] J. C. Weaver and Y. Chizmadzhev. Theory of electroporation: A review. *Bioelectrochemistry and Bioenergetics*, 41(2), 1996.

[40] J. C. Weaver, K. C. Smith, A. T. Esser, R. S. Son, and T. Gowrishankar. A brief overview of electroporation pulse strength duration space: A region where additional intracellular effects are expected. *Bioelectrochemistry*, 87, 2012.