M Leguèbe

C Poignard

L Weynans
email: lisl.weynans@inria.fr

A second-order Cartesian method for the simulation of electropermeabilization cell models

Keywords: Finite differences on cartesian grids, Non-linear partial differential equations, Cell modeling Méthodes de différences finies sur grilles cartésiennes, Equations aux dérivées partielles non linéaires, Modélisation bio-cellulaire

In this paper we present a new finite-difference method to simulate electropermeabilization models, like the model of Neu and Krassowska [10] or the recent model of Kavian et al. [6]. These models are based on the evolution of the electric potential in conducting media. The main feature lies in the transmission of the voltage potential which is proportional to the electric flux (the normal derivative of the potential). An adapted scheme is necessary to accurately simulate the voltage potential in the whole cell, notably at the membrane separating the cell from the outer medium. Our finite difference scheme is inspired by the method of Cisternino and Weynans [3] for elliptic problems with immersed interfaces. This is a Cartesian grid method based on the accurate discretization of the fluxes at the interface, through the use of additional interface unknowns. The aim of the paper is to present the method and to study its stability and convergence properties. We prove the first-order convergence of the method for the dynamic non-linear model in one dimension, and the second order convergence in the stationary linear one-dimensional case. This results are corroborated by numerical experiments in two dimensions.

Introduction

In this paper, we aim at providing an efficient numerical method to compute the voltage potential in a biological cell submitted to high electric pulses. We propose a second order Cartesian method to simulate the electric cell model of Kavian et al. but our method is applicable for other problems involving specific transmission conditions at the interface.

Motivations

Electropermeabilization (also called electroporation) is a significant increase in the electrical conductivity and permeability of the cell membrane that occurs when pulses of large amplitude (a few hundred volts per centimeter) are applied to the cells: due to the electric field, the cell membrane is permeabilized, and then non-permeant molecules can easily enter the cell cytoplasm by diffusion through the electropermeabilized membrane areas. If the pulses are too long, too numerous, or if their amplitude is too high, the cell membrane is irreversibly destroyed and the cells are killed. However, should the pulse duration be sufficiently short (a few milliseconds or a few microseconds, depending on the pulse amplitude) then the cell membrane reseals within several tens of minutes: this is termed reversible electroporation, which preserves the cell viability and is used in electrochemotherapy to vectorize the drugs until the cell inside. Many clinical studies (phase II and phase III) have proven that electrochemotherapy of cutaneous or subcutaneous metastases or tumors, either with bleomycin or with cisplatin displays an objective response rate of more than 80%. Reduction of tumor size has been achieved with electrochemotherapy faster and more efficiently than in standard chemotherapy for both cutaneous and subcutaneous tumors. Therefore the efficacy of such a treatment has been proven for superficial metastases. However, despite cell electropermeabilization by micropulses is well-known from the experimental point of view, there is a lack of numerical results that are coherent with the experiments, preventing a systematic use of electropermeabilization in configurations far from the experiments, especially for deep tumors. This lack is mainly due to the complexity of the electropermeabilization models at the cell scale, which involve partial differential equations with non-standard transmission conditions through the cell membrane. Actually, the probably most achieved model has been derived in the 2000's by Neu, Krassowska et al. from physical considerations. However, the model involved many parameters, that prevents fitting to the experiments. Moreover the theoretical analysis of this model, which can help in the derivation of accurate numerical schemes, seems to be hardly achievable due to the unboundedness of several unknows. Roughly speaking, the current models provide a qualitative explanation of the electropermeabilization, but the problem of the quantitative description remains open. Very recently, Kavian et al. have provided a new phenomenological model involving few parameters, and they show the solvability of their model [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF].

The aim of this paper is to present an accurate numerical method to compute the voltage potential in a biological when electropermeabilizing electric is applied. The numerical method is a finite difference method on Cartesian grid inspired from the paper of Cisternino and Weynans [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] but adapted to the features of the electropermeabilization models.

Electropermeabilization model

We first recall the electropermeabilization model of Kavian et al. Due to its very thin thickness and its very low conductivity the cell membrane can be modeled as a surface electric material Γ with a capacity C and a surface conductivity S. The electroporation model of Kavian et al. describes the surface conductivity through a sliding-door type model [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF].

U (0, •) = 0 in O e ∪ O c , (1a)
and for t > 0

∆U = 0, in (0, T ) × (O e ∪ O c ) , (1b)
U (t, •) = g(t, •) on (0, +∞) × ∂Ω, (1c)
with the jump conditions

[σ∂ n U ] = 0, on (0, T ) × Γ, (1d)
C ∂ t [U ](t, •) + S(t, [U ])[U ] = σ c ∂ n U (t, •)| Γ -, on (0, T ) × Γ. (1e)
The surface membrane conductivity S is defined as an interpolation between the two values S ir and S L , which are respectively the surface conductivity of an irreversibly electropermeabilized region and the lipid surface conductivity:

S(t, λ) := S L + (S ir -S L )X(t, λ), (1f)
where the function X satisfies the differential equation for t > 0 (we set λ(t) := λ(t, s) := [U (t, s)] for a.e. s ∈ Γ):

   ∂X(t, λ) ∂t = max β(λ(t)) -X(t, •) τ ep ; β(λ(t)) -X(t, •) τ res . X(0, λ) = 0, (2)

RR n°8302

where β is an even regularized Heaviside function. For instance β can be chosen as

for all λ ∈ R, β(λ) := (1 + tanh(k ep (|λ| -V rev ))/2.
Despite jump condition (1d) is nothing but the flux continuity, the transmission condition on the jump of the potential (1d) reflects the influence of the thin resistive membrane. Actually, the jump of the potential is linked to the electric flux through a Kadam-Katchalsky type condition and therefore the solution U is not continuous across the interface Γ. Note that in condition (1e), it is necessary to implicit the flux ∂ n u| Γ -in the time discretization for stability reasons, otherwise even for the linear model (which consists in S(t, λ) ≡ S L ) a drastic Courant-Friedrichs-Lax condition appears, previously observed by Guyomarc'h et al. [START_REF] Guyomarc | A discontinuous Galerkin method for elliptic interface problems with application to electroporation[END_REF]. More precisely, in their paper they used a forward Euler time scheme (equation (26) pp 1005 of [START_REF] Guyomarc | A discontinuous Galerkin method for elliptic interface problems with application to electroporation[END_REF]), for which a very small time step has to be chosen. In [START_REF] Guyomarc | A discontinuous Galerkin method for elliptic interface problems with application to electroporation[END_REF], a discontinuous Galerkin method is proposed for the space discretization. We present here a finite difference method on a Cartesian grid, for which the flux of (1e) is implicit. Despite the use of Cartesian grid seems to less accurate than the meshing of the interface when dealing with one single cell, we emphasize that such Cartesian grids are useful for a large number of cells. Note that the main goal of this research is to fit the parameters of the model with the experiments, and in vitro experiments even for dilute cell solutions deal with few hundreds of cell per square centimeter, thus Cartesian grid and level-set method seem more appropriate than accurate meshing of each cell.

Outline of the paper

The paper is organized as follows. Section 2 describes the finite difference method used to solve numerically the electropermeabilization model. In section 3, we first prove that, in one dimension, this method is convergent with second order accuracy for the static linear model in space. We then study in section 4 the convergence of the method in the dynamical non-linear case. We prove that, in one dimension, the method converges with first order accuracy in both time and space. For the sake of clarity, we reduce the theoretical assessement to the one dimension. In section 5, we present two-dimensional numerical validations that corroborate the results of the previous sections. In particular we show that for the linear problem, the second-order accuracy is reached, while the method for the dynamical non-linear case converges with first order accuracy in time. We conclude the paper by perspectives on the applications of our numerical scheme.

Statement of the method

Description of the interface and classification of grid points

We choose to perform the discretization on a Cartesian grid covering O e ∪ O c . The points on this grid are defined by M ij = (x i , y j ) = (ih, jh), with h the grid spacing. We denote by u ij the approximation of U at the point (x i , y j ). If necessary, the additional subscript u n ij will be used to indicate the time t n = n dt, where dt is the time step. Similarly, if it is necessary to distinguish in which subdomain the grid point is, we will use the notations u c ij and u e ij . In order to describe accurately the geometric configuration in the vicinity of the interface we use the level set method introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]. We refer the interested reader to [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], [START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] for recent reviews of this method. The zero isoline of the level set function, defined by the signed function ϕ:

ϕ(x) = dist Σ (x) outside of the interface -dist Σ (x) inside of the interface (3)
Inria represents implicitely the interface Σ immersed in the computational domain. A useful property of the level set function is

n(x) = ∇ϕ(x) |∇ϕ(x)| , (4)
where n(x) is the outward normal vector of the isoline of ϕ passing on x. This allows us to compute the values of the normal to the interface. We say that a grid point is neighbouring the interface if the sign of ϕ changes between this point and at least one of its neighbours, see Fig. 2. On the contrary, grid points far from the interface are said regular grid points. If the intersection of the interface and the segment [M ij M i+1j ] exists, then we define the interface point I i,j,E = (x i,j,E , y j ) as this intersection. In the same way, if the intersection of the interface and [M i-1j M ij ] exists, then we define the interface point I i,j,W = (x i,j,W , y j ) as this intersection.

Similarly, the interface points I i,j,N = (x i , ỹi,j,N ) and I i,j,S = (x i , ỹi,j,S ) are respectively defined as the intersection of the interface and the segments

[M ij M ij+1 ] and [M ij-1 M ij ].
Let us remark that, with this notation, the same interface point is described in two different ways: for instance I i,j,S = I i,j-1,N or I i,j,E = I i+1,j,W .

At each interface point we create two additional unknowns, called interface unknowns, and denoted by ũi,j,γ with γ = E, W, N or S. The interface unknowns carry the values of the numerical solution on each side of the interface.

Global discretization of the system

We choose a first order semi-implicit temporal discretization of the equations on the electric potential U , and an explicit first order Euler discretization for X. The computational domain is a square of characteristic size 1, with N points in both x-and y-directions. We denote by (∆u) ij the finite-differences approximation of the Laplacian at point M i,j . The approximated normal derivative at the interface point I i,j,γ in the exterior (resp. interior domain), with γ = E, W, N or S, is denoted by (∂ n u e ) i,j,γ (resp. (∂ n u c ) i,j,γ ). We write the discretized system on the Cartesian grid:

• The Laplace operator on the grid points:

∀ (i, j) ∈ {1, .., N } 2 , ∀ n 0, (∆u) n+1 ij = 0. (5)
• The jump conditions (1d) :

∀ (i, j) ∈ {1, .., N } 2 , ∀ n 0, if I i,j,γ exists, then [σ(∂ n u) n+1 i,j,γ ] = 0. (6)
• The transmembrane potential evolution (1e):

∀ (i, j) ∈ {1, .., N } 2 , ∀ n 0, if I i,j,γ exists, then C [u] n+1 i,j,γ -[u] n i,j,γ dt + S n i,j,γ [u] n i,j,γ = σ c (∂ n u c ) n+1 i,j,γ . (7)
where S n i,j,γ , with γ = E, W, N or S, is given by

S n i,j,γ = S L + S ir X n i,j,γ ,
• The coefficient X:

∀ (i, j) ∈ {1, .., N } 2 , ∀ n 0, if I i,j,γ exists, then X n+1 i,j,γ = X n i,j,γ + dt max β([u] n i,j,γ ) -X n i,j,γ τ ep ; β([u] n i,j,γ ) -X n i,j,γ τ res . (8)
Initial conditions are u 0 = 0 and X 0 = 0. We impose Dirichlet boundary conditions on the outer boundary of the computational domain. What we call in the following "linear problem" corresponds to the case where S is kept constant in time, in which case no electropermeabilization happens.

Remark 2.1. It is also possible to discretize the transmembrane potential evolution as

C [u] n+1 i,j,γ -[u] n i,j,γ dt + S n i,j,γ [u] n+1 i,j,γ = σ c (∂ n u c ) n+1 i,j,γ ,
This option can be less convenient because it requires to solve at each time step a linear system with a different matrix, while the discretization that we have chosen uses the same matrix at each iteration. But on the other side, because it is more implicit, it may require weaker stability conditions on the time step. The analysis of this alternative discretization can be carried on in a similar way to the discretization that we have chosen.

In the following subsection we explain how we discretize the terms (∆u) n ij , (∂ n u e ) n i,j,γ and (∂ n u c ) n i,j,γ with γ = E, W, N or S. Once these discretizations are achieved we obtain a linear system that is solved with a standard linear solver such as restarted GMRES.

The discretization is inspired by the finite difference second order accurate method developed in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] for elliptic problems with interfaces. The second-order convergence has been validated in two dimensions on various geometries in the latter reference. This method is in the same spirit than the well-known Immersed Interface Method (IIM) of LeVeque and Li [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF]. The IIM relies on Taylor expansions of the solution on each side of the interface, with a local coordinate transformation near the interface to express the jump conditions in an appropriate frame. The elliptic operator is discretized on each grid point near the interface with formulas accounting for the jumps across the interface. In order to find these formulas a linear system with six unknowns needs to be solved for each of the concerned grid points. The method is second order accurate in maximum norm. Numerous developments of the IIM have been performed, among them: the fast IIM algorithm of Li [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] for elliptic problems with piecewise constant coefficients, the Explicit Jump Inria Immersed Interface Method (EJIIM) by Wiegmann and Bube [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF], the Decomposed Immersed Interface Method (DIIM) by Bethelsen [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF], and the MIIM (maximum principle preserving) by Li and Ito [START_REF] Li | Maximum principle preserving schemes for interfacxe problems with discontinuous coefficients[END_REF].

Another class of Cartesian method recently introduced by Zhou et al. is the Matched Interface and Boundary (MIB) method: [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF]. This method can provide finite-difference schemes of arbitrary high order. The solution on each side of the interface is extended on fictitious points on the other side. These fictitious values are computed by iteratively enforcing the lowest order interface jump conditions. Chern and Shu [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] also proposed a Coupling Interface Method, where the discretizations on each subdomain are coupled through a dimension by dimension approach using the jump conditions. All the methods cited above are second order accurate.

Contrary to these methods, the method introduced in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] is based on the use of additional interface unknowns rather than a non-standard discretization of the Laplacian directly accounting for the discontinuities near the interface. The interface unknowns are particularly useful in the case of the electropermeabilization model, since the value of the transmembrane potential is accurately needed. It is however possible, if necessary, to eliminate them from the linear system by using the Schur complement technique. Directly placing points on the interface is also useful to this model, as some data are only defined on the membrane, such as S. This property of the method allows for example to compute integrals on the interface with much efficacy, as no interpolation between grid values is needed.

Discrete elliptic operator

On the regular grid points, we use the standard centered finite differences scheme to approximate the Laplacian:

(∆u) ij = u i+1j -2u ij -u i-1j h 2 + u ij+1 -2u ij -u ij-1 h 2 . ( 9
)
This formula is second-order accurate:

(∆u) ij = ∆U (x i , y j ) + O(h 2 ).
For a grid point M ij neighbouring the interface, we compute (∆u) ij with the values on M ij and the closest points (grid or intersection points) to M ij in each direction. For instance for the situation illustrated on Fig. 3(a) we get:

(∆u) ij = ũi,j,E -u ij xi,j,E -x i - u ij -u i-1j h (x i,j,E -x i ) + h 2 + u ij+1 -2u ij -u ij-1 h 2 . ( 10
)
This discretization is first-order accurate: It is straightforward to compute a second order approximation of the x-derivative of the potential with three a priori non equidistant points. For example, we approximate the x-derivative on the left side of the interface with the points M i-1j , M ij and I i,j,E by: We thus have

(∆u) ij = ∆U (x i , y j ) + O(h).

Discretization of fluxes at the interface

(∂ x u) i,j,E = (u i-1j -ũi,j,E )(x i -xi,j,E ) h (x i-1 -xi,j,E ) - (u ij -ũi,j,E )(x i-1 -xi,j,E ) h (x i -xi,j,E ) . ( 11
) j j+1 j-1 i i-1 i+1 i+2 I i,j,E (
(∂ x u) i,j,E = ∂U ∂x (x i,j,E , y j ) + O(h 2 ).
The right x-derivative ∂U ∂x (x i+1,j,W , y j ) is approximated by (∂ x u) i+1,j,W which is computed similarly.

For the derivative along the y-direction we do not have unknowns located on the line parallel to the y-axis passing by I i,j,E . We therefore use a linear combination of (∂ y u) ij and (∂ y u) i-1j , defined respectively as second order approximations of the y-derivative on M ij and M i-1j . We obtain

(∂ y u) i,j,E = xi,j,E -x i-1 h (∂ y u) ij - xi,j,E -x i h (∂ y u) i-1j . ( 12
)
The resulting formula is second-order accurate:

(∂ y u) i,j,E = ∂U ∂y (x i,j,E , y j ) + O(h 2 ).
The formulas for (∂ y u) ij and (∂ y u) i-1j depend on the local configuration on the interface, but are based on the same principle as for [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]. The formulas [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] and ( 12) are consistent if the Inria point M i-1,j belongs to the same domain as M ij . Consequently they require that there are at least two adjacent points in each direction belonging to the same domain. Otherwise we use a first order discretization involving only three points of this side of the interface: see Fig. 3(c).

The normal derivatives at the interface are discretized by the scalar product of x-and yderivatives with the normal to the interface:

(∂ n u) i,j,γ = (∂ x u) i,j,γ n x + (∂ y u) i,j,γ n y , with γ = E, W, N, S, (13)
with (n x , n y ) an approximation of the vector normal to the interface at point I i,j,γ . Numerical instabilities may appear when the discretization of the flux at an interface point involves another intersection point closely located to a grid point. This situation is illustrated in Fig. 4(a). In these cases we use a decentered discretization: instead of using the second intersection point we use the closest grid point located in the opposite direction. For instance, on Fig. 4(b), the term (∂ y u c ) i,j,E is computed with the grid points (x i+1 , y j ), (x i+1 , y j+1 ) and (x i+1 , y j+2 ) instead of (x i+1 , y j ), (x i+1 , y j+1 ) and I i+1,j,S . 3 Second order convergence in linear static unidimensional case

In this section, we study the convergence of the method for the linear static model. In subsections 3.1 and 3.2 we respectively give the notations that will be used hereafter and detail the linear system arising from the discretization of the method. Then in subsection 3.3 we prove the monotonicity of the matrix of the linear system. At the end, in subsection 3.4 we use the monotonicity to estimate the value of the coefficients of the inverse matrix, and we prove the second-order convergence of the method.

Notations

Let us give some additional notations about the spatial configuration that we consider, cf Fig. 5: the interior domain and the exterior domain are

O c = (x L , x R ) and O e = (0, x L ) ∪ (x R , 1).
The outer boundary and the interface are composed of two points:

∂Ω = {0} ∪ {1} and Γ = {x L } ∪ {x R } .
The interface point x L is located between the grid points x kL-1 and x kL and let

x L = x kL -αh, with 0 ≤ α ≤ 1.
Similarly, the interface point x R is located between the grid points x kR and x kR+1 and we set x R = x kR + βh, There are 2 interface unknowns on each interface point, one interior and one exterior: u e L , u c L u c R and u e R . The computational domain is discretized on a uniform grid, where the grid points x i are defined by

with 0 ≤ β ≤ 1.
x i = i h, with h = 1/(N + 1).

Linear system

The fluxes at the interface are computed with the formulas:

(∂ n u c ) L = - 1 + 2α α(α + 1)h (u kL -u c L ) + α (1 + α)h (u kL+1 -u kL ), (14)
(∂ n u c ) R = 1 + 2β β(β + 1)h (u c R -u kR ) - β (β + 1)h (u kR -u kR-1 ), (15)
and

(∂ n u e ) L = - 3 -2α (1 -α)(2 -α)h (u e L -u kL-1 ) + 1 -α (2 -α)h (u kL-1 -u kL-2 ), (16)
(∂ n u e ) R = 3 -2β (1 -β)(2 -β)h (u kR+1 -u e R ) - 1 -β (2 -β)h (u kR+2 -u kR+1 ). (17)

Inria

All these flux formulas are on the same principle: if you have three points x 1 , x 2 and x 3 , then the coefficients allocated to these points in order to approximate with second-order accuracy the derivative in x 1 are respectively c 1 , c 2 and c 3 such that:

c 1 + c 2 + c 3 = 0, (18a)
x 1 c 1 + x 2 c 2 + x 3 c 3 = 1 or -1, (18b)
x 2 1 c 1 + x 2 2 c 2 + x 2 3 c 3 = 0. ( 18c
)
The sign in the above relationships depends on the orientation of the outer normal to the interface.

The jump conditions can be expressed as

σ c (∂ n u c ) L -σ e (∂ n u e ) L = 0, σ e (∂ n u e ) R -σ c (∂ n u c ) R = 0, S(u e L -u c L ) -σ c (∂ n u c ) L = 0, σ c (∂ n u c ) R -S(u e R -u c R ) = 0,
but we rather rewrite them as

S(u e L -u c L ) -σ e (∂ n u e ) L = 0, (19a)
S(u e R -u c R ) -σ e (∂ n u e ) R = 0, (19b)
σ c (∂ n u c ) L -S(u e L -u c L ) = 0, (19c)
σ c (∂ n u c ) R -S(u e R -u c R ) = 0, (19d)
so that the matrix of the linear system be monotone.

The discretization of the Laplacian operator reads

-u 0 + 2u 1 -u 2 h 2 = 0, -u 1 + 2u 2 -u 3 h 2 = 0, . . . . . . -u kL-3 + 2u kL-2 -u kL-1 h 2 = 0, -2 u e L -u kL-1 (1 -α)(2 -α)h 2 + 2 u kL-1 -u kL-2 (2 -α)h 2 = 0, -2 u kL+1 -u kL (1 + α)h 2 + 2 u kL -u c L α(1 + α)h 2 = 0, -u kL + 2u kL+1 -u kL+2 h 2 = 0, . . . . . . -u kR-2 + 2u kR-1 -u kR h 2 = 0, -2 u c R -u kR β(1 + β)h 2 + 2 u kR -u kR-1 (1 + β)h 2 = 0, -2 u kR+2 -u kR+1 (2 -β)h 2 + 2 u kR+1 -u e R (1 -β)(2 -β)h 2 = 0, -u kR+1 + 2u kR+2 -u kR+3 h 2 = 0, . . . . . . -u N -2 + 2u N -1 -u N h 2 = 0, -u N -1 + 2u N -u N +1 h 2 = 0.

Preliminary results

In the following, we denote by A the discretization matrix of the linear system described in the previous subsection. We first prove that A is monotone by showing that all the coefficients of its inverse are positive. Let u be an array of size N + 4 corresponding to N grid points and four interface unknowns such that Au 0, and let us show that u 0.

3.3.a Sign of the normal interface derivative

We begin by proving that if the minimum of u is located on an interior interface point, then the normal derivative at this interface point is negative. Similarly, if the minimum of u is located on an exterior interface point then the normal derivative at this interface point is positive. We do it for example for the normal derivative on x L in the exterior domain. The other cases would be treated exactly in the same way. The exterior normal derivative at x L is

(∂ n u e ) L = - 3 -2α (1 -α)(2 -α)h (u e L -u kL-1 ) + 1 -α (2 -α)h (u kL-1 -u kL-2 ).

Inria

By hypothesis, Au 0 hence

u kL-1 -u kL-2 h u e L -u kL-1 (1 -α)h , therefore (∂ n u e ) L - 3 -2α (1 -α)(2 -α)h (u e L -u kL-2 ) + 1 -α (2 -α) u e L -u kL-1 (1 -α)h , -3 + 2α + 1 -α (1 -α)(2 -α)h (u e L -u kL-1 ), - 1 (1 -α)h
(u e L -u kL-1 ), 0.

3.3.b Monotonicy of the matrix

The minimum of u can either be reached on a grid point in the interior or exterior subdomains, or on an interface point, interior or exterior.

• if the minimum is reached on a grid point in the interior domain:

In this case we denote by i min the index of the smallest component of u. We first assume that i min = 1, N . Using the Laplacian inequality on this point:

-u imin+1 + 2u imin -u imin-1 h 2 0,
we deduce that u imin-1 and u imin+1 are equal to u imin , and, using then all other Laplacian inequalities in the interior subdomain, we deduce that all values in the interior subdomain are equal to the minimum value u imin . Therefore, the approximate normal derivative is zero. Thus, using discrete relationships (19a) -(19d), the jump of the numerical solution and the normal exterior derivative are all negative:

0 = σ c (∂ n u c ) L S(u e L -u c L ) σ e (∂ n u e ) L , 0 = σ c (∂ n u c ) R S(u e R -u c R ) σ e (∂ n u e ) R .
We consider for instance the situation on the right interface point. From S(u e R -u c R ) 0 we deduce that u e R = u imin . We can also write

(∂ n u e ) R = 3 -2β (1 -β)(2 -β)h (u kR+1 -u e R ) - 1 -β (2 -β)h (u kR+2 -u kR+1 ) 0, thus 3 -2β (1 -β)(2 -β)h (u kR+1 -u e R ) 1 -β (2 -β)h (u kR+2 -u kR+1 ).
Moreover, due to the Laplacian inequality,

u kR+2 -u kR+1 u kR+1 -u e R 1 -β , thus 3 -2β (1 -β) (u kR+1 -u e R ) u kR+1 -u e R .
Therefore u kR+1 -u e R 0, from which we deduce that u kR+1 = u e R = u imin and thus u kR+2 = u e R = u imin also. Using the Laplacian inequality in the exterior subdomain we deduce that all point values in the right exterior subdomain are equal to u imin . The Laplacian inequality on the first grid point yields:

2u 1 -u 2 h 2 0,
therefore all the values of u are positive.

• if the minimum is reached on a grid point in the exterior domain:

Let us assume for instance that it is located on the left side. The right side would be treated the same way. Following the same reasoning as above, we deduce directly that all values in the left exterior subdomain are equal and positive. Therefore, as the minimum is positive, then all values of u are positive.

• if the minimum is reached on one exterior interface point.

Let us assume for instance that it is reached on u e R . We have u e R -u c R 0 and due to (19a) -(19d), we can write:

0 S(u e L -u c L ) σ e (∂ n u e ) R .
Now we use the same reasoning as in the case where the minimum is on a grid point in the interior subdomain and conclude that all values in the right exterior subdomain are positive and equal to the minimum value. Consequently all values of u are positive.

• if the minimum is reached on one interior interface point. Suppose for instance that it is located on x L . We then have that the normal derivative (∂ n u c ) L is negative. Due to (19a) -(19d), we can write:

0 σ c (∂ n u c ) L S(u e L -u c L ) σ e (∂ n u e ) L .
Therefore u e L = u c L , and thus (∂ n u e ) L 0. We conclude exactly with the same reasoning than in the previous that all values of u are positive.

We conclude that the matrix A is monotone and thus invertible.

3.3.c Estimation of

A -1 ∞
Since all the coefficients of A -1 are positive we can write A -1 ∞ = max

1 i N N j=1
a -1 i,j . We want to build an array v such that

Av =      1 1 . . . 1     
.

in order to have:

N j=1 a -1 i,j = v i .
We look for v as a quadratic function with respect to x.

Inria

For values on the left exterior subdomain, including the interface point values, we set

v i = -x 2 i /2 + Cx i + D,
for values on the interior subdomain, including the interface point values v i is such that

v i = -x 2 i /2 + C ′ x i + D ′ ,
and for values on the right exterior subdomain, including the interface point values, we write

v i = -x 2 i /2 + C ′′ x i + D ′′ .
We check that the discrete Laplacian equalities are satisfied for

(v i ) i=2,N -1 - v i+1 -2v i + v i-1 h 2 = 1 for all i ∈ {2, k L -2} ∪ {k L + 1, k R -1} ∪ {k R + 2, N -1} - v e L -v kL-1 (1 -α)h - v kL-1 -v kL-2 h (2 -α)h 2 = 1, - v kL+1 -v kL h - v kL -v c L αh 1 + α 2 = 1, - v c R -v kR βh - v kR -v kR-1 h 1 + β 2 = 1, - v kR+2 -v kR+1 h - v kR+1 -v e R (1 -β)h (2 -β)h 2 = 1.
In order to satisfy the same relationship for the lines v 1 and v N , corresponding to the two ends of the domain, we impose

D = 0, D ′′ = -C ′′ + 1/2.
For the lines corresponding to the jump conditions, inserting ( 14)-(17) in the relationships (19) we get

S(v e L -v c L ) + σ e 3 -2α (1 -α)(2 -α)h v e L -σ e 2 -α (1 -α)h v kL-1 + σ e 1 -α (2 -α)h v kL-2 = 1, S(v e R -v c R ) -σ e 2 -β (1 -β)h v kR+1 + σ e 3 -2β (1 -β)(2 -β)h v e R + σ e 1 -β (2 -β)h v kR+2 = 1, -σ c 1 + α αh v kL + σ c 1 + 2α α(α + 1)h v c L + σ c α (1 + α)h v kL+1 -S(v e L -v c L ) = 1, σ c 1 + 2β β(β + 1)h v c R -σ c 1 + β βh v kR + σ c β (β + 1)h v kR-1 -S(v e R -v c R ) = 1.
Due to the relationships (18), the above expressions simplify into:

S ((C -C ′ )x L + D -D ′ ) + σ e C = 1, S ((C ′′ -C ′ )x R + D ′′ -D ′ ) -σ e C ′′ = 1, σ c C ′ -S ((C -C ′ )x L + D -D ′ ) = 1, -σ c C ′ -S ((C ′′ -C ′ )x R + D ′′ -D ′ ) = 1.
With such coefficients C, D, C ′ , D ′ , C ′′ and D ′′ that do not depend on h we have

A -1 ∞ = max 1 i N v i max C 2 2 + D, C ′2 2 + D ′ , C ′′2 2 + D ′′ . RR n°8302

Proof of second-order convergence

Proposition 3.1. The numerical solution to the linear stationary problem converges with a second-order accuracy to the exact solution to the problem.

Proof. We denote by z the i -th row of A -1 , i.e. z j = A -1 i,j . We use the same indexing system for z as for the unknowns of the linear system. We define the Kronecker symbol for interface values as:

δ L,c i = 1 if i is the index of the left interior interface value, 0 otherwise, δ R,c i = 1 if i
is the index of the right interior interface value, 0 otherwise, δ L,e i = 1 if i is the index of the left exterior interface value, 0 otherwise, δ R,e i = 1 if i is the index of the right exterior interface value, 0 otherwise. We multiply the row vector z with the columns of A.

The columns 1 to k L -1 yield:

2z 1 -z 2 = h 2 δ 1 i , -z 1 + 2z 2 -z 3 = h 2 δ 2 i , . . . . . . -z kL-4 + 2z kL-3 -z kL-2 = h 2 δ kL-3 i , -z kL-3 + 2z kL-2 - 2 (2 -α) z kL-1 + σ e h 1 -α 2 -α z e L = h 2 δ kL-2 i , -z kL-2 + 2 (1 -α)(2 -α) z kL-1 + 2 (2 -α) z kL-1 -σ e h 2 -α 1 -α z e L = h 2 δ kL-1 i . The columns k L + 1 to k R -1 yield: -z kL+1 + 2 (1 + α) z kL + 2 α(1 + α) z kL -σ c h 1 + α α z c L = h 2 δ kL i , 2z kL+1 -z kL+2 - 2 (1 + α) z kL + σ c h α 1 + α z c L = h 2 δ kL+1 i , -z kL+1 + 2z kL+2 -z kL+3 = h2δ kL+2 i , . . . . . . -z kR-3 + 2z kR-2 -z kR-1 = h 2 δ kR-2 i , -z kR-2 + 2z kR-1 - 2 (1 + β) z kR + σ c h β β + 1 z c R = h 2 δ kR-1 i , -z kR-1 + 2 β(1 + β) z kR + 2 (1 + β) z kR -σ c h 1 + β β z c R = h 2 δ kR i .
Inria and the columns k R + 1 to N :

2 2 -β z kR+1 + 2 (1 -β)(2 -β) z kR+1 -σ e h 2 -β 1 -β z e R -z kR+2 = h 2 δ kR+1 i , 2z kR+2 -z kR+3 - 2 2 -β z kR+1 + σ e h 1 -β 2 -β z e R = h 2 δ kR+2 i , -z kR+2 + 2z kR+3 -z kR+4 = h 2 δ kR+3 i , . . . . . . -z N -2 + 2z N -1 -z N = h 2 δ N -1 i , -z N -1 + 2z N = h 2 δ N i .
The columns corresponding to interface values yield:

S + σ e 3 -2α (1 -α)(2 -α)h z e L - 2 (1 -α)(2 -α)h 2 z kL-1 -Sz c L = δ L,e i , -Sz e L + S + σ c 1 + 2α α(α + 1)h z c L - 2 α(1 + α)h 2 z kL = δ L,c i , (20)
- 2 β(1 + β)h 2 z kR + S + σ c 1 + 2β β(β + 1)h z c R -Sz e R = δ R,c i , (21)
-Sz c R + S + σ e 3 -2β (1 -β)(2 -β)h z e R - 2 (1 -β)(2 -β)h 2 z kR+1 = δ R,e i .
We now want to prove that all the coefficients of z corresponding to the interior subdomain, that is, z j with k L j k R , are O(h). The lines (20) and (21) can be rewritten:

S + σ c 1 + 2α α(α + 1)h z c L = 2 α(1 + α)h 2 z kL + Sz e L + δ L,c i , S + σ c 1 + 2β β(β + 1)h z c R = 2 β(1 + β)h 2 z kR + Sz e R + δ R,c i .
Because we have proved that A -1 ∞ is bounded independently of h we can write that z e L , z c L , z c R and z e R are all O(1). We also know, since the matrix A is monotone, that all the coefficients z j are positive. Consequently, z kL and z kR are O(h). Now, if we sum the lines corresponding to columns k L + 1 to k R -1 we get:

z kL+1 - 2 1 + α z kL + σ c h α 1 + α z c L + z kR-1 - 2 1 + β z kR + σ c h β β + 1 z c R = kR-1 i=kL+1 h 2 δ j i .
Then

z kL+1 + σ c h α 1 + α z c L + z kR-1 + σ c h β β + 1 z c R = 2 1 + α z kL + 2 1 + β z kR + kR-1 i=kL+1 h 2 δ j i ,
from which we infer that z kL+1 and z kR-1 are O(h). Similarly, doing successive sums, we prove that all the coefficients z j with k L j k R are O(h). We would prove similarly the same result for the others coefficients z j with 1 j k L -1 or k R + 1 j N but for the sake of brevity we do not write the proof here.

We eventually conclude that all the coefficients of the row i of A -1 at the exception of z e L , z c L , z c R and z e R are of order O(h), without any coefficient of the type 1 α , 1 1-α , 1 β or 1 1-β that may be singular if α or β tend to 0 and lead to instabilities in the convergence. Consequently the array containing the difference between the numerical solution and the exact solution on each discretization point, which is equal to the product of the array of the truncation errors on each point by A -1 , is of order O(h 2 ), which ensures the second-order accuracy of the method.

Convergence result for the one-dimensional dynamical model

In subsection 4.1 we rewrite the linear problem to solve as an interface problem. Then in subsection 4.2 we prove that in one dimension the matrix appearing in the formulation of the interface problem is symmetric and its spectral radius is bounded by one. Theses properties are actually crucial to prove the stability of the numerical method. Then we prove that the scheme is converging with first order accuracy, using techniques similar to those used to prove the existence and uniqueness of the solution to the model of Kavian et al. [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF].

Formulation of the discrete problem as an interface problem

In order to analyze the convergence of the numerical method to the exact solution to the electropermeabilization model, we need to rewrite the linear system as an interface problem. To this purpose we rewrite the whole linear system by separating the unknowns in the interior domain from the unknowns in the exterior domain. We denote by u e g and u c g the respective unknowns on grid points in the exterior and interior domains, and u e p and u c p the unknowns on interface points in the exterior and interior domains. This leads to The two first block lines represent the discretization of the Laplacian operator on the exterior and interior grid points. ∆ e g and ∆ e p are the block matrices corresponding to the discretization of the Laplacian in the exterior subdomain, and ∆ c p and ∆ c g the block matrices corresponding to the discretization of the Laplacian in the interior subdomain. Remark that ∆ e g and ∆ c g are invertible since they correspond to the discretization of the Laplacian on grid points in the exterior and interior domain. The source term g contains the terms accounting for boundary conditions.

The third block line represents the discretization of the flux equality across the interface. More precisely,

Φ e g u e g + Φ e p u e p = σ e ∂ n u e , Φ c g u c g + Φ c p u c p = -σ c ∂ n u c ,
where ∂ n u e and ∂ n u c are vectors containing respectively the discretization of the normal derivative to the interface in the exterior and interior domains. The last block line represents the Inria evolution of the jump of the solution across the interface. S n is the diagonal matrix containing the values of (S) n i,j,γ at each interface point. Using the invertibility of ∆ e g and ∆ c p we write

u e g = -(∆ e g ) -1 (∆ e p u e p -g), u c g = -(∆ c g ) -1 ∆ c p u c p .
Injecting these expressions in the discrete approximations of fluxes yields

σ e ∂ n u e = -Φ e g (∆ e g ) -1 (∆ e p u e p -g) + Φ e p u e p , σ c ∂ n u c = Φ c g (∆ c g ) -1 ∆ c p u c p -Φ c p u c p ,
and we define the matrices Λ c , Λ e and Λ 0 as

σ c Λ c = Φ c g (∆ c g ) -1 ∆ c p -Φ c p , σ e Λ e = -Φ e g (∆ e g ) -1 ∆ e p + Φ e p , σ e Λ 0 = Φ e g (∆ e g ) -1 .
Let us notice that Λ c and Λ e are discrete Dirichlet-to-Neumann operators: indeed

Λ c : v → Λ c v = -Φ c g v g -Φ c p v with ∆ c g v g = -∆ c p v.
The vector v g corresponds to grid values in the interior domain, and is the solution to the discrete Laplacian in the interior domain with v as Dirichlet boundary conditions on the interface.

Λ e : v → Λ e v = Φ e g v g + Φ e p v with ∆ e g v g = -∆ e p v.
The vector v g corresponds to grid values in the exterior domain, and is the solution to the discrete Laplacian in the exterior domain with v as Dirichlet boundary conditions on the interface and homogeneous Dirichlet conditions on the outer boundary.

We have

σ e Λ 0 g + σ e Λ e u e p + σ c Λ c u c p = 0. (23)
We do not write the proof here for the sake of brevity, but one can prove that Λ e is invertible. Let us notice also that the matrix σ e Λ e + σ c Λ c is invertible. Indeed, it is the linear application that associates to a vector w p representing the values at the interface points (both exterior and interior, assuming that there is no discontinuity across the interface for these values), the jump of the exterior and interior normal derivatives to the interface of the solution to the linear problems

∆ c g w c g = -∆ c p w p , ∆ e
g w e g = -∆ e p w p .

Consequently, the matrix

Id + σc σe Λ -1 e Λ c = 1 σe Λ -1
e σ e Λ e + σ c Λ c is also invertible. Therefore, using (23) and denoting [u] = u e p -u c p we write

u e p = 1 σ e Λ -1 e (-σ e Λ 0 g + σ c Λ c u c p ), u c p = u e p -[u] = - σ c σ e Λ -1 e Λ c u c p + Λ -1 e Λ 0 g -[u],
Therefore

u c p = Id + σ c σ e Λ -1 e Λ c -1 -[u] + 1 σ e Λ -1 e Λ 0 g . (24)
The matrix -

Id + σc σe Λ -1 e Λ c -1
is the linear application that associates to an array [v], representing a jump across the interface, the array v c p representing the interior values at the interface, satisfying the linear problem:

  ∆ e g 0 ∆ e p 0 0 ∆ c g 0 ∆ c p Φ e g Φ c g Φ e p Φ c p       u e g u c g u e p u c p     =   0 0 0   . with v e p -v c p = [v]. The temporal evolution of [u] can be expressed as C dt [u] n+1 -σ c Λ c u c p = C dt Id -S n [u] n . (25)
The evolution of [u] can be rewritten using formula (24) as

Id + dt C σ c Λ c Id + σ c σ e Λ -1 e Λ c -1 [u] n+1 = Id - dt C S n [u] n + dtF, (26)
with

F = -σc C Λ c Id + σc σe Λ -1 e Λ c -1
Λ -1 e Λ 0 g. In the following, for the sake of brevity, we will note

M = Id + dt C σ c Λ c Id + σ c σ e Λ -1 e Λ c -1 -1
.

The relationship (26) can thus be re-written Note that ρ(Id

[u] n+1 = M Id - dt C S n [u] n + dt M F. (27)
-dt C S n ) < 1 if |1 -dt C S| < 1 for all values of (S) i,j,γ with γ = E, W, S, N, which is true notably if dt < 2 C SL+Sir . Remark 4.2. The expression [u] n+1 = M Id - dt C S n [U ] n + dtF,
is linearly equivalent to the initial linear system (25), which is consistent with the equations to solve. Therefore, by injecting the exact solution ũ in the latter formula, one can write

[ũ](t n+1,x ) = M Id - dt C S(X(t n , x)) [ũ](t n , x) + dtF + O(dt) , ∀x ∈ Γ.
It is proven in [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF] that the exact solution is bounded in C([0, T ], L 2 (Γ)), therefore we deduce that F 2 = O(1).

Inria

Convergence result in one-dimension

We address the convergence of the numerical solution to the exact solution to the non-linear problem. We denote by [u] n the vector containing the values of the jump of the exact solution at the interface points at time t n , and [u] n h the vector containing the values of the jump of the numerical solution at the same points and same time. X n and X n h are defined similarly as the exact and approximated vector values of X at these interface points. We denote the error in discrete L 2 norm at time step n by

E n = [u] n -[u] n h 2 .
Because our method is formulated in the finite-differences formalism, we need to assume that the solutions to the electropermeabilization model considered are smooth enough to have the truncation errors tending to zero when the grid spacing h tends to zero.

The geometrical configuration and notations are the same as in section 3, see Fig. 5.

Lemma 4.3. In one-dimension, the matrix

Λ c Id + σc σe Λ -1 e Λ c -1
is symmetric and all its eigenvalues are positive.

Proof. We begin by proving that Λ c is symmetric and has only positive eigenvalues, i.e.

∀v p , w p , (Λ c v p , w p ) = (Λ c w p , v p ) and (Λ c v p , v p ) ≥ 0.
We denote by v L , and w L two values on x L , v R and w R two values on x R . Let v g and w g be the solutions to the linear problems

∆ c g v g = -∆ c p v L v R and ∆ c g w g = -∆ c p w L w R .
The arrays -Λ c v p and -Λ c w p represent the discrete normal derivatives at the interface. We thus can write: Φ c g v g + Φ c p v p = -Λ c v p and Φ c g w g + Φ c p w p = -Λ c w p . In one dimension we can write:

v g =      v kL v kL+1 . . . v kR     
, and

w g =      w kL w kL+1 . . . w kR      .
We sum over every grid point i inside the interior domain the product (∆v) i w i where (∆v) i is the discrcete Laplacian of v defined by (∆ c g v g + ∆ c p v p ) i . For the indexes k L and k R , the Laplacian term is also multiplied by (α + 1)/2 or (β + 1)/2 in order to obtain the same denominator as the interior terms, which will simplify the computations. Then we make a summation by parts.

(v kL+1 -v kL ) h 2 w kL - (v kL -v L ) αh 2 w kL + kR-1 i=kL+1 (v i+1 -v i ) -(v i -v i-1 ) h 2 w i + (v R -v kR ) βh 2 w kR - (v kR -v kR-1 ) h 2 w kR = 0, - (v kL -v L ) αh 2 w kL - kR-1 i=kL (v i+1 -v i ) h 2 (w i+1 -w i ) + (v R -v kR ) βh 2 w kR = 0.
We use the fact that the discrete Laplacian is zero at points x kL and

x kR : v R -v kR βh - v kR -v kR-1 h = 0, v kL+1 -v kL h - v kL -v L αh = 0.
Furthermore, the discrete second order normal derivative to the interface at points x L and x R are

(∂ n v) L = -(∂ x v) L = 1 + 2α α(α + 1)h (v kL -v L ) - α (1 + α)h (v kL+1 -v kL ), (∂ n v) R = (∂ x v) R = 1 + 2β β(β + 1)h (v R -v kR ) - β (β + 1)h (v kR -v kR-1 ).
With the help of these latter equations we prove that

(∂ n v) L = - (v kL -v L ) αh , (28)
(∂ n v) R = (v R -v kR ) βh . (29)
Thus,

(∂ n v) L h w kL + (∂ n v) R h w kR = kR-1 i=kL (v i+1 -v i ) h 2 (w i+1 -w i ). (30)
We can consequently write

(∂ n v) L h (w L -αh(∂ n w) L ) + (∂ n v) R h (w R -βh(∂ n w) R ) = kR-1 i=kL (v i+1 -v i ) h 2 (w i+1 -w i ),
which leads to

(∂ n v) L h w L + (∂ n v) R h w R = kR-1 i=kL (v i+1 -v i ) h 2 (w i+1 -w i ) + α(∂ n v) L (∂ n w) L + β(∂ n v) R (∂ n w) R . (31)
We conclude that

(∂ n v) L h w L + (∂ n v) R h w R = (∂ n w) L h v L + (∂ n w) R h v R , which is equivalent to (Λ c v p , w p ) = (Λ c w p , v p ),
and therefore Λ c is symmetric. Furthermore, taking w i = v i for all i in (31) we obtain (Λ c v p , v p ) ≥ 0 which proves that Λ c has only positive eigenvalues. We would prove similarly the same result for Λ e . v.

Now we consider again the matrix

Λ c Id + σc σe Λ -1 e Λ c
Λ c Id + σ c σ e Λ -1 e Λ c -1 v • v = Λ c Id + σ c σ e Λ -1 e Λ c -1 v • Id + σ c σ e Λ -1 e Λ c Id + σ c σ e Λ -1 e Λ c -1 v, = Λ c w • Id + σ c σ e Λ -1 e Λ c w, = Λ c w • w + σ c σ e Λ e (Λ -1 e Λ c w) • (Λ -1 e Λ c ) ≥ 0.
Inria Lemma 4.4. If dt < τ ep , then for all n ≥ 0 and on each interface point, we have 0 ≤ X n h ≤ 1. Proof. The value of X n+1 h is computed from X n h and [u] n h with a first order explicit Euler scheme. The formula reads

X n+1 h = X n h + dt max β([u] n h ) -X n h τ ep , β([u] n h -X n h τ res .
We use the notation

X n+1 h = X n h + dt β([u] n h ) -X n h τ n ,
with τ n = τ ep or τ res depending on the sign of β([u] n h ) -X n h . We define (X n h ) -= max(0, -X n h ) and multiply the latter equation by (X n+1 h ) -+ (X n h ) -. We obtain:

(X n+1 h ) -+ (X n h ) -X n+1 h -X n h dt = (X n+1 h ) -+ (X n h ) -β([u] n h ) -X n h τ n , X n+1 h (X n+1 h ) --X n h (X n h ) - dt = (X n+1 h ) -+ (X n h ) -β([u] n h ) -X n h τ n + X n h (X n+1 h ) - dt - X n+1 h (X n h ) - dt , - ((X n+1 h ) -) 2 -((X n h ) -) 2 dt = (X n+1 h ) -+ (X n h ) -β([u] n h ) τ n - X n h (X n h ) - τ n + X n h (X n+1 h ) -1 dt - 1 τ n - X n+1 h X n h dt .
The terms on the first line of the right hand-side are positive, notably because β 0. We assume that dt τ n . If X n h > 0, or if X n h < 0 and X n+1 h < 0 then the term on the second line of the right hand-side is also positive and (X n+1 h ) -< (X n h ) -. On the contrary, if X n h > 0 and X n+1 h > 0 then one has also (X n+1 h ) -< (X n h ) -. Consequently the sequence ((X n h ) -) n 0 is decreasing. Because (X n h ) -0 for all n 0 and (X 0 h ) -= 0 one concludes that X n h 0 for all n 0.

With the same reasoning, by considering (X n h -1) + = max(0, X n h -1) and using the fact that β -1 0, we can prove that X n h 1 for all n 0. Lemma 4.5. For all T > 0 there exists a constant K(T ) depending on T , β and τ ep but not on other parameters, such that on every interface point the following inequality is satisfied for all

0 n T dt sup i n X i -X i h K(T ) sup i n [u] i -[u] i h + O(dt). (32)
Proof. The function β ′ is bounded on [-∞, +∞], and we denote by

β ′ max = sup |β ′ (x)|. X n+1 -X n+1 h |X n -X n h | + |ε n 2 | + dt max β([u] n ) -X n τ ep , β([u] n ) -X n τ res -max β([u] n h ) -X n τ ep , β([u] n h ) -X n τ res + dt max β([u] n h ) -X n τ ep , β([u] n h ) -X n τ res -max β([u] n h ) -X n h τ ep , β([u] n h ) -X n h τ res .
Let us observe that

max β([u] n ) -X n τ ep , β([u] n ) -X n τ res -max β([u] n h ) -X n τ ep , β([u] n h ) -X n τ res |β([u] n ) -β([u] n h )| τ ep , β ′ ∞ |[u] n -[u] n h | τ ep ,
and similarly

max β([u] n h ) -X n τ ep , β([u] n h ) -X n τ res -max β([u] n h ) -X n h τ ep , β([u] n h ) -X n h τ res |X n -X n h | τ ep . Consequently X n+1 -X n+1 h (1 + dt τ ep ) |X n -X n h | + dt β ′ ∞ τ ep |[u] n -[u] n h | + |ε n 2 | .
Now we use a discrete Gronwall lemma to obtain

|X n -X n h | exp n dt τ ep X 0 -X 0 h + n-1 i=0 exp n dt -i dt τ ep dt β ′ ∞ τ ep [u] i -[u] i h + ε i 2 , n-1 i=0 exp T τ ep dt β ′ ∞ τ ep [u] i -[u] i h + ε i 2 , sup i n X i -X i h β ′ ∞ T τ ep exp T τ ep n-1 i=0 dt [u] i -[u] i h + n-1 i=0 exp T τ ep ε i 2 ,
from which we deduce the result of the lemma. is symmetric and has only positive eigenvalues. As we have proved in lemma 4.4 that 0

X 1 then S > 0. Therefore M = Id + dt C σ c Λ c Id + σc σe Λ -1 e Λ c -1 -1
is also symmetric and:

M 2 = ρ(M ) 1.
Using (27) we write

[u] n+1 h 2 ρ(M ) [u] n h 2 + dt ρ(M ) F 2 , [u] 0 h 2 + dt n i=0 F 2 [u] 0 h 2 + T F 2 .
Therefore [u] n h 2 is bounded independently of n and h. In one-dimension the interface values are only two, therefore the boundedness in L 2 norm implies the boundedness in L ∞ norm.

Inria

Proposition 4.7. If the ratio dt/h is kept constant, if dt < τ ep and 1 -dt C S < 1 for all values of (S) i,j,γ with γ = E, W, S, N, then in one-dimension the numerical method converges with first-order accuracy to the exact solution.

Proof. Let us recall the formula used to compute

[u] n+1 h and X n+1 h from [u] n h and X n h : [u] n+1 h = M Id - dt C S(X n h ) [u] n h + dtF , (33)
X n+1 h = X n h + dt max β([u] n h ) -X n h τ ep , β([u] n h ) -X n h τ res (34)
We know that the exact solution satisfies

[u] n+1 = M Id - dt C S(X n ) [u] n + dtF + ε n 1 ( 35
)
with ε n 1 related to the truncation error of the equation. If we assume that the ratio dt/h is kept constant, we know that ε n 1 = O(dt 2 ) because the discretization in first-order accurate in time.

[

u] n+1 -[u] n+1 h = M Id - dt C S(X n ) [u] n + dtF -Id - dt C S(X n h ) [u] n h + dtF + ε n 1 , = M ([u] n -[u] n h ) - dt C M (S(X n )[u] n -S(X n h )[u] n h ) + ε n 1 , [u] n+1 -[u] n+1 h 2 [u] n -[u] n h 2 + dt C (S L + S ir ) [u] n -[u] n h 2 + dt C [u] n h ∞ S(X n ) -S(X n h ) 2 + ε n 1 2 , [u] n+1 -[u] n+1 h 2 1 + dt C (S L + S ir ) [u] n -[u] n h 2 + dt C K(T )S ir [u] n h ∞ sup i n [u] n -[u] n h 2 + ε n 1 2 + O(dt 2 ). [u] n+1 -[u] n+1 h 2 1 + dt C K(T )S ir sup i n [u] i h ∞ + dt C (S L + S ir ) sup i n [u] i -[u] i h 2 + ε n 1 2 + O(dt 2 ), sup i n+1 [u] i -[u] i h 2 1 + dt C K(T )S ir sup i n [u] i h ∞ + dt C (S L + S ir ) 1+dt K(T ) sup i n [u] i -[u] i h 2 + ε n 1 2 + O(dt 2 ).
With a discrete Gronwall lemma, we deduce that

sup i n [u] i -[u] i h 2 e T K(T ) [u] 0 -[u] 0 h 2 =0 + n i=0 e (ndt-idt) K(T ) ε i 1 2 + O(dt 2 ) =O(dt)
.

We conclude that the method converges with a first-order accuracy.

RR n°8302 5 Numerical validations

In this section, we study numerically the convergence of the method in two dimensions, for both static and dynamic cases, with and whithout electropermeabilization (corresponding to nonlinear and linear cases respectively). For an analysis of the "physical" behavior of the model we refer to the numerical studies in [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF].

Parameters used in the simulations

Table 1 presents the parameters used for the following simulations. Some of them are biological parameters, chosen in accordance with reference studies in the literature, the other ones are specific to the model described in subsection 1.2. Figure 6 provides an example of results that are usually obtained in our simulations, for an irregularly shaped cell in this case. The electric field is visualized with lines while the membrane conductivity is depicted by a white to black scale and the thickness of the membrane. Note that this thickness is only a handful visualization artifact to locate the permeabilized regions of the membrane.

Convergence study for the linear static problem

In order to determine the convergence order of our method, we derive the analytic solution of the static problem as reference in a comparison with our results. When considering a circular cell or radius R 1 , centered on the origin, inside a concentric domain Ω of radius R 2 , R 2 > R 1 , and by choosing a constant conductivity S = S L , one can express such a solution of the linear problem. Let the source g = ER 2 cos θ with E tuning the amplitude of the electric field. with α e , β e , and α c given by

α c = σ c S L R 1 + 1 + σ c σ e R 2 + σ c S L R 1 + 1 - σ c σ e R 2 1 R 2 -1 g, α e = 1 2
σ c S L R 1 + 1 + σ c σ e α c , β e = 1 2
σ c S L R 1 + 1 - σ c σ e α c R 2 1 .
This solution is used as a Dirichlet condition on the boundary of our simulation box, which is a square contained in the disk of radius R 2 , and containing the disk of radius R 1 (cf. figure 7).

The error between the result u of our simulation and the exact solution U is measured on both grid and interface points. In order to compute a discretized L 2 norm on the interface, the latter is considered linear between intersection points. If M m-1 , M m and M m+1 designate three consecutive interface points, the length element ∆l m associated to the point M m is chosen as:

∆l m = 1 2 ( M m+1 -M m + M m -M m-1 )
.

RR n°8302

• u e i-2j

• u e i-1j

• u e ij • u e ij+1 • u e ij+2 • u c ij • u e i-1jE • u e ij+1S A e ij A c ij O e O c

RR n°8302 6 Conclusion

In this paper we have presented a new finite differences method to simulate electropermeabilization models based on elliptic PDE's describing the evolution of the electric potential in conducting media. It is a Cartesian grid method based on the accurate discretization of the fluxes at the interface between the exterior domain and the cell, representing the cell membrane. This accurate discretization is performed with additional interface unknowns, located at the intersection between the interface and the grid axes. The use of interface unknowns is particularly relevant in the context of electropermeabilization, because the jump in electric potential across the cell membrane is explicitly given in the models and is measured by experiments.

We have studied the stability and convergence of the method in one dimension, and proved that in one-dimension, it converges with a first-order accuracy for the dynamic nonlinear problem, and with a second-order accuracy for the linear stationary problem. The first order-accuracy in the non-linear case is due to the fact that a first-order scheme is employed for the time evolution. Finally we have presented numerical results in two dimensions corroborating the theoretical study.

In the future, we intend to use this method to perform simulations of the electropermeabilization of high numbers of cells in the same computational domain. To this purpose we will develop a parallelized version of the method, with the help of the PETSc library. Furthermore, we plan to modify the method to take into account a possible variation of the volume of the electropermeabilized cells, which is highlighted by experiments. This will be handled by computing the temporal evolution of the level-set function, following the equation of evolution for the interface that will be added to the model.

4

  Convergence result for the one-dimensional dynamical model 4.1 Formulation of the discrete problem as an interface problem . . . . . . . . . . . . 4.2 Convergence result in one-dimension . . . . . . . . . . . . . . . . . . . . . . . . . 5 Numerical validations 5.1 Parameters used in the simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Convergence study for the linear static problem . . . . . . . . . . . . . . . . . . . 5.3 Convergence study for the linear dynamic problem . . . . . . . . . . . . . . . . . 5.4 Convergence study with electropermeabilization . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Geometry of the problem. The cell O c is imbedded in the bath O e . The whole domain Ω is defined by Ω = O e ∪ O c .

Figure 2 :

 2 Figure 2: Example of geometrical configuration, with regular grid points in black, points neighbouring the interface in white, and interface points with both possible notations.

Fig. 3 (

 3 Fig. 3(b) and 3(c) provide examples of interface configurations.It is straightforward to compute a second order approximation of the x-derivative of the potential with three a priori non equidistant points. For example, we approximate the x-derivative on the left side of the interface with the points M i-1j , M ij and I i,j,E by:

Figure 3 :

 3 Figure 3: Examples of discretization of the elliptic operator (Fig. 3(a)), of second (resp. first) order flux discretization (Fig. 3(b)) (resp. Fig. 3(c)) at the interface point I i,j,E .

Figure 4 :

 4 Figure 4: Centered and decentered stencils. In Fig. 4(a) the discretization of the flux on the point I i,j,E involves the intersection point I i+1,j,S and is ill-conditionned. Fig. 4(b) shows the decentered stencil: the discretization of the flux on the point I i,j,E involves grid point (i+1, j +2) instead of the intersection point I i+1,j,S .

Figure 5 :

 5 Figure 5: Geometrical configuration in one dimension.

Proposition 4 . 1 .

 41 Problem (22) and problem (27) are equivalent.

- 1 .

 1 Let v be an array of the same size as [u]. We denote by w = Id + σc σe Λ -1 e Λ c -1

Figure 6 :

 6 Figure 6: Numerical electric field and permeabilization of a cell with an irregular shape. The iso-values of the electric field are visualized with colored lines while the membrane conductivity is depicted by a white to black scale and the thickness of the membrane.

Figure 8 :

 8 Figure 8: L 2 norm calculation for grid points. u e ij is an extrapolation of the solution from the values of u e in O e to the point M ij in O c , involving interface points as well.

Figure 9 :Figure 10 :

 910 Figure 9: Error vs relative grid spacing h = 1/N x = 1/N y for the linear static case, with N x and N y being the number of points in the x and y directions.

Figure 11 :

 11 Figure 11: Error E h versus grid spacing h, in both static (11(a)) and dynamic (11(b)) case.

  Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Electropermeabilization model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Outline of the paper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Description of the interface and classification of grid points . . . . . . . . . . . . 2.2 Global discretization of the system . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Discrete elliptic operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Discretization of fluxes at the interface . . . . . . . . . . . . . . . . . . . . . . . Proof of second-order convergence . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents
	1 Introduction
	1.1 2 Statement of the method
	2.1

3 Second order convergence in linear static unidimensional case 3.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Linear system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.a Sign of the normal interface derivative . . . . . . . . . . . . . . . . . . . . 3.3.b Monotonicy of the matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.c Estimation of A -1 ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4

  This condition is satisfied for instance if we take dt < 2 C SL+Sir . Proof. We first have to prove that [u] n h 2 is bounded independently of n and h.

	σe Λ -1 e Λ c Because of lemma 4.3, the matrix Λ c Id + σc	-1

Lemma 4.6. The sequences [u] n h 2 and [u] n h ∞ are bounded independently of n and h in onedimension if 1 -dt C S < 1 for all values of (S) i,j,γ with γ = E, W, S, N.

Table 1 :

 1 Parameters set to fit to the results given by[START_REF] Neu | Asymptotic model of electroporation[END_REF][START_REF] Debruin | Modelling electroporation in a single cell. I. Effects of field strength and rest potential[END_REF]. (EP stands for electropermeabilization, and EPd stands for electropermeabilized).

	Variable	Symbol Value	Unit
	Biological parameters:
	Extracellular conductivity	σ e	5	S/m
	Intracellular conductivity Capacitance	σ c C	0.455 9.5 × 10 -3 F/m 2 S/m
	Membrane surface conductivity	S L	1.9	S/m 2
	Cell radius	r	50	µm
	Membrane thickness	δ	5	nm
	Specific parameters of the model:
	EP threshold	V rev	1,5	V
	EP switch speed	k ep	40	V -1
	EP characteristic time	τ ep	1 × 10 -6	s
	Resealing characteristic time EPd membrane surface conductivity	τ res S ir	1 × 10 -3 2.5 × 10 8	s S/m 2

RR n°8302

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

This research has been partly granted by the French national agency throughout the research projects INTCELL (2010-BLAN-916-04) and MEMOVE (2011-BS01-006-01).

Domain of calculation of the analytic solution

Computational domain

Figure 7: Numerical estimation of the order of accuracy of the method. The analytic solution to the linear problem is calculated in concentric circular domains. The restriction of this solution to the boundary of the computational domain (dashed line) provides the Dirichlet data for the numerical solution.

The error on the interface is defined by:

where N l is the number of unknowns at interface.

For grid points, a special treatment is needed for points neighbouring the interface for a second order L 2 norm calculation. Indeed, as shown in Fig. 8, considering the solution constant on a virtual cell of size h x × h y around each grid point is not valid, as these cells may intersect with both inner and outer domains. Let A c ij (resp. A e ij ) be the area of the part of the cell which is in O c (resp. O e ), which can be null if the cell is entirely in one domain. If the point M ij is located in O c (resp. O e ), the solution in the exterior domain O e (resp. O c ), as well as the solution U , are extrapolated on M ij , with respective values u e ij and U e ij (resp. u c ij and U c ij ). The error on grid points is then defined by:

The total error reads:

For the static problem, a second-order of convergence is achieved, as presented in Fig. 9.

Inria

Convergence study for the linear dynamic problem

The solution of the dynamical problem has the same form as the solution (36) of the static equation, with time-dependent coefficients:

Then the jump on the interface [U ] satisfies the following ordinary differential equation (ODE):

Using the continuity of the flux (1d) we infer: that

Therefore the jump is given by:

from which we infer all the coefficients. The error is computed at the final time of the simulation T f = 100µs, similarly as for the static equation. Fig. 10 shows that a second-order of convergence is obtained for this case as well.

Convergence study with electropermeabilization

When electropermeabilization is taken into account, i.e. when the conductivity S depends on the TMP difference [U ], we are not able to compute an analytic solution. In order to compute the order of convergence p of our method, a solution u h computed with a grid spacing h is compared to u h/2 , computed with a refined grid spacing h/2. The error E h is then of the order of h p , with:

E h is only computed on grid points, since interface points may not match on both grids. On Fig. 11 are presented the convergence results in the static and dynamic case. In the static case, the numerical solution is obtained with a pseudo-time iterative scheme described in [START_REF] Kavian | Classical" electropermeabilization modeling at the cell scale[END_REF]. In the dynamic case, we take k ep = 10 rather than 40, because the electropermeabilization happens early with the choice k ep = 10, helping us to perform the convergence study in time. E h is measured after that the electropermeabilization has begun, at the final time T f = τep 2 . The ratio dt dx is kept constant, equal to τ ep . Let us notice that it is smaller than the one used in the linear case, due to the non-linearity. Convergence is achieved with an order one accuracy, which is coherent with the fact that the temporal integration scheme is only order one.

Inria