O Kavian

M Leguèbe

C Poignard
email: clair.poignard@inria.fr

L Weynans

"Classical" Electropermeabilization Modeling at the Cell Scale

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

In this paper, we derive two models (a static and a dynamical model) based on the description of the electric potential in a biological cell in order to model the cell electropermeabilization. Existence and uniqueness results are provided for each model, and an accurate numerical method to compute the solution is described. We then present numerical results that corroborate the experimental results, which tends to justify the validity of our modeling. We emphasize that our new models involve very few parameters, compared with the most achieved model but they describe the same qualitative results. Moreover our numerical results are quantitavely close to the experimental data.

Introduction

The distribution of the electric potential in a biological cell is important for bio-electromagnetic investigations. A sufficiently large magnitude of the difference of transmembrane potential (denoted from now on by ∆TMP), which is the difference of the electric potentials between both sides of the cell membrane, leads to an increase of the membrane permeability [START_REF] Pucihar | Numerical determination of transmembrane voltage induced on irregularly shaped cells[END_REF][START_REF] Teissié | Mechanisms of cell membrane electropermeabilization: A minireview of our present (lack of?) knowledge[END_REF]. Molecules such as bleomycin can then diffuse across the plasma membrane. This phenomenon, called electropermeabilization, has already been used in oncology and holds promises in gene therapy [START_REF] Marty | Electrochemotherapy -an easy, highly effective and safe treatment of cutaneous and subcutaneous metastases: Results of esope (european standard operating procedures of electrochemotherapy) study[END_REF][START_REF] Serša | Application of electroporation in electrochemotherapy of tumors[END_REF], motivating precise assessments of the ∆TMP.

In this paper, we aim at studying theoretically and numerically two non-linear electrical models (a static and a dynamical model) of biological cells. These models are inspired from the static model of Ivorra et al. in [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF]. They describe the behavior of both electric potential and membrane conductivity when the cell is submitted to an electric pulse through a few parameters that will be fitted with the experiments. The inverse problem solving to calibrate the models, which is the main goal of this research, is not tackled in the present paper. This article is a first step, in which we present, and study theoretically and numerically new models with a few parameters, in order to simplify the forthcoming inverse problem solving. We emphasize these models are phenomenological in the sense that the membrane conductivity is described by an ad hoc law, which does not come from an homogenization of the nanoscale phenomena. Before stating the model we are going to study, we now detail the notation of the paper. Notation 1.1. Throughout this article, we shall use the following conventions and notation:

• We generically denote by n the normal to a closed smooth surface of R 3 (or a curve of R 2 ) outwardly directed from the inside to the outside of the domain enclosed by the surface.

• Let C be a surface embedded in R 3 , and let u be a sufficiently smooth function (in an appropriate sense) defined in a tubular neighborhood of C. We define u | C ± by ∀x ∈ C, u | C ± (x) = lim τ →0 + u(x ± τ n(x)).

The notation ∂ n u | C ± and ∂ t u | C ± stands for the normal and the tangential components of ∇u:

∀x ∈ C, ∂ n u | C ± (x) = lim τ →0 + ∇u(x ± τ n(x)) • n(x), ∂ t u | C ± (x) = ∇u(x) -∂ n u(x)n(x),
where the dot " • " denotes the Euclidean scalar product of R 3 . In the case of R 2 , the analogous notation is easily adapted.

• The jump [u] C of a function u defined in a neighborhood of the surface C is defined by

[u] C = u| C + -u| C -.

Electric potential in a biological cell

A biological cell is a high contrast medium composed of a conducting cytoplasm O c surrounded by a thin and very insulating layer, embedded in a bath O e (see Fig. 1). The plasma membrane is a phospholipid bilayer, which is sprinkled over with proteins. Due to its thickness and its electrical properties, the membrane can be modeled as a surface electric material Γ with a capacity C m and a surface conductivity S m . We refer to the seminal articles of Hodgkin, Goldman, Katz, et Inria al. for the electric description of cell membranes [START_REF] Goldman | Potential, impedance and rectification in membranes[END_REF][START_REF] Hodgkin | The effect of sodium ions on the electrical activity of the giant squid axon[END_REF][START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Hodgkin | The influence of potassium and chloride ions on the membrane potential of single muscle fibres[END_REF]. Let σ be the conductivity of the medium, that is As described by Neu, Krassowska, and DeBruin [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] the electric potential in the whole cell is the discontinuous solution U to the following problem:

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

U (t, x) = g(t, x), on ∂Ω, (1c)
with the transmission conditions across the membrane Γ:

[σ∂ n U ] Γ = 0, (1d)
C m ∂ t [U ] Γ + S m [U ] Γ = σ c ∂ n U | Γ -. (1e)
Throughout, the paper we denote by Ω the Lipschitz domain

Ω = O e ∪ O c
and ∂Ω denotes its boundary. Moreover we define P H 1 (Ω) as

P H 1 (Ω) = u ∈ L 2 (Ω) : u| Oe ∈ H 1 (O e ), u| Oc ∈ H 1 (O c ) ,
and P H 1 0 (Ω) = u ∈ P H 1 (Ω) : u| ∂Ω = 0 .

Electropermeabilization phenomenon

When submitted to a high electric pulse -i.e. if the magnitude of the pulse g reaches a threshold value -the cell membrane permeability increases and large molecules that usually cannot diffuse through the plasma membrane (for instance, plasmids or bleomycin) enter inside the cytoplasm.

This phenomenon is called electroporation or electropermeabilization. For several years, different membrane models based on hydrodynamic, elasticity, hydroelasticity, viscoelasticity, or aqueous pore formation have been developed to describe the pore formation on the cell membrane (for more details, see the review of Pavlin et al. [START_REF] Pavlin | Electroporation of planar lipid bilayers and membranes[END_REF]). They all highlight a threshold value of the electric potential above which the electropermeabilization phenomenon occurs. However, the critical potential value changes with the models. Theoretical biophysicists consider the aqueous pore formation model as the most convincing current explanation. Nevertheless, the predictions of the model do not coincide with experiments either quantitatively or phenomenologically since no pore has ever been observed. In addition, the models based on the paper of Neu, Krassowska, et al. [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][START_REF] Debruin | Modelling electroporation in a single cell. II. Effects of ionic concentrations[END_REF] are too complex to be parameterized to fit the experiments. Roughly speaking, the current models provide a qualitative explanation of the electropermeabilization, but the problem of the quantitative description remains open.

Actually in vitro and in vivo experiments have never proved the electropore formation, which theoretically could reach detectable size (since macropores could be created according to Smith et al. [START_REF] Smith | Model of creation and evolution of stable electropores for DNA delivery[END_REF]); and it seems unclear whether electroporation results from holes punched in a lipid bilayer, as proposed in the current models (Teissié et al. [START_REF] Teissié | Mechanisms of cell membrane electropermeabilization: A minireview of our present (lack of?) knowledge[END_REF][START_REF] Teissié | In Vitro cell electropermeabilization[END_REF]). Moreover, the experimentally proved reversible process of the membrane electroporation is not clearly explained by the current models. In addition, and this is probably one of the main features of the experiments, the vectorization of large molecules requires both short time high-voltage pulses and long time low-voltage pulses (André et al. [START_REF] André | Efficiency of high -and low -voltage pulse combinations for gene electrotransfert in mucle, liver tumor and skin[END_REF]). Therefore, the presence of pores is still controversial despite structural changes of the membrane. For all these reasons, we prefer the term electropermeabilization to electroporation.

Modeling principle

The probably most achieved electropermeabilization model has been proposed by Neu and Krassowska [START_REF] Neu | Asymptotic model of electroporation[END_REF] and precisely described by DeBruin and Krassowska [START_REF] Debruin | Modelling electroporation in a single cell. I. Effects of field strength and rest potential[END_REF][START_REF] Debruin | Modelling electroporation in a single cell. II. Effects of ionic concentrations[END_REF]. It consists in adding an electroporation current I ep in the right-hand side of equation (1e):

C m ∂ t [U ] Γ + S m [U ] Γ = σ c ∂ n U | Γ -+ I ep .
The current I ep is given by a highly non-linear pore current i ep multiplied by the pore density N ep . The main drawback of such a model is its complexity since several parameters, such as the pore radius or the relative entrance length of the pore, which cannot be measured. Moreover, the mathematical well-posedness of the equations cannot be clearly established. Therefore, the inverse problem, that consists in fitting the parameters for each cell species, is hardly unsolvable numerically. In addition, the "philosophy" of the modeling is based on the pore creation, while as hinted above the very existence of these pores is controversial. For all these reasons, we choose to present here a new phenomenological model of electropermeabilization that could be fitted with the experiments. This model describes the membrane resealing and the memory of the applied pulses.

Our electropermeabilization modeling consists in describing the membrane permeabilization by choosing an appropriate function for the surface conductivity S m , instead of adding an electroporation current based on the pore creation as Neu, Krassowska, et al. did (see [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][START_REF] Smith | Model of creation and evolution of stable electropores for DNA delivery[END_REF][START_REF] Neu | Singular perturbation analysis of the pore creation transient[END_REF]). In addition, we propose two models: the first one is a "static" model that describes the electropermeabilization as being an instantaneous phenomenon for a single time-constant pulse. This can be seen as a preliminary model that describes the cell potential during the pulse. The second model is the time-dependent model of electropermeabilization. For each model, we present the theoretical results that ensure existence and uniqueness of the solution to the new problems Inria and then we present the numerical methods that allow the computation of the equations. We conclude by comparing our model with the model of Neu, Krassowska, et al.

The static equation

Based on the extensive review of Ivorra et al [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF], the surface conductivity S m is a function of the absolute membrane voltage, which tends to the value S L (the lipid surface conductivity) below a certain threshold V rev (the reversible electropermeabilization voltage) and tends to S ir (the surface conductivity of the irreversibly electropermeabilized region) above this threshold, with S ir being larger than S L . The "speed of the switch" between these two values is given by a parameter k ep . We may choose the following sigmoid function for S m :

∀λ ∈ R, S m (λ) = S L + (S ir -S L )[1 + tanh(k ep (|λ| -V rev ))]/2, (2)
however other functions with similar monotonicity properties can be considered. More precisely, in our model the function S m will satisfy the following condition:

S m ∈ C(R), λ → S m (λ) is even on R, 0 < S L ≤ S m (λ) ≤ S ir , S m is non decreasing on [0, +∞), lim λ→+∞ S m (λ) = S ir .        (3)
In particular, note that the mapping λ → λS m (λ) is increasing on R. Therefore, the static potential U satisfies the following problem:

∆U = 0, in O e ∪ O c , (4a)
[σ∂ n U ] Γ = 0, on Γ, (4b)
S m ([U ] Γ ) [U ] Γ = σ c ∂ n U | Γ -, on Γ, (4c)
U = g on ∂Ω. ( 4d
)
Remark 2.1. Model (4) can be seen as the limit of the model of Ivorra et al. [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF], when the membrane thickness tends to zero (we refer to [START_REF] Perrussel | Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer[END_REF] for asymptotic expansion of the voltage potential in high contrast medium with resistive thin layer).

In the following subsections, we study the non-linear problem (4). In particular, we emphasize that due to the non-linearity of the membrane conductivity, increasing numerically the thickness of the membrane (as performed in [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF]) leads to irrelevant results from the quantitative point of view. Therefore, we aim at providing efficient numerical methods in order to solve the above problem.

Existence and uniqueness of the static potential

This subsection is devoted to the proof of the following result. Theorem 2.2. Let g ∈ H 1/2 (∂Ω). There exists a unique U satisfying problem (4). This solution satisfies

U | Oe ∈ H 1 (O e ), U | Oc ∈ H 2 (O c ).
In order to prove this theorem, we proceed as follows. Denote by Λ c and Λ e the Dirichlet-to-Neumann operators on Γ (also called Steklov-Poincaré operators) for the Laplacian respectively in O c and in O e . More precisely, denote by n c (resp. n e ) the unitary outward normal to Γ directed from the inside to the outside of O c (resp.O e ). We define the operators Λ c and Λ e from H 1/2 (Γ) to H -1/2 (Γ) as:

∀f ∈ H 1/2 (Γ), Λ c (f ) := n c • σ c ∇v c | Γ -, where div(σ c ∇v c ) = 0 in O c , and v c | Γ = f, (5a)
Observe that using Wirtinger-Poincaré's inequality in the case of Λ c , or Poincaré's inequality in the case of Λ e , together with the continuity of the mapping

u → u |∂O from H 1 (O) into H 1/2 (∂O),
when O is sufficiently smooth, the following inequalities hold:

Λ c f, f = Oc σ c (x)∇v(x) • ∇v(x)dx ≥ C c f -M (f ) 2 H 1/2 (Γ) , (6)
Λ e f, f = Oe σ e (x)∇v e (x) • ∇v e (x)dx ≥ C e f 2 H 1/2 (Γ) , (7)
where

M (f ) = |Γ| -1 Γ f (τ ) dτ
is the mean value of f on Γ, and C e and C c are constants depending only on O e and O c respectively. Moreover, for a function g ∈ H 1/2 (∂Ω), we define Λ 0 (g) by:

Λ 0 (g) := n e • σ e ∇v | Γ + , where div(σ e ∇v) = 0 in O e , v | ∂Ω = g and v |Γ = 0. (8)
It is useful to recall that the operator Λ e is invertible, its inverse being given by another Consider the Hilbert space H defined by

H = H 1/2 (Γ) × H 1/2 (Γ),
with the norm

∀ u = (u e , u c ) ∈ H, u 2 H = u e 2 H 1/2 (Γ) + u c 2 H 1/2 (Γ) .
Problem (4) can be written on the manifold Γ with the help of the above Steklov-Poincaré operators. More precisely, problem (4) is equivalent to finding (u e , u c ) ∈ H such that

Λ e u e + S m (u e -u c )(u e -u c ) = -Λ 0 (g), Λ c u c -S m (u e -u c )(u e -u c ) = 0. ( 9
) Notation 2.3. Identifying the dual of L 2 (Γ) × L 2 (Γ) with L 2 (Γ) × L 2 (Γ)
, we denote by H the dual space of H and by • , • the duality between H and H .

The proof of Theorem 2.2 is an obvious application of the following theorem.

Theorem 2.4. Let G = (G e , G c ) ∈ H . There exists a unique u 0 = (u e , u c ) ∈ H such that Λ e u e + S m (u e -u c )(u e -u c ) = G e , Λ c u c -S m (u e -u c )(u e -u c ) = G c . (10)

Inria

Proof. We define the operator

Λ σ from H into H by ∀u ∈ H, Λ σ u = Λ e u e Λ c u c = Λ e 0 0 Λ c u e u c . (11)
Thanks to ( 6) and ( 7), we have

∀u ∈ H, Λ σ u, u ≥ C e u e 2
H 1/2 (Γ) + C c u c -M (u c ) 2 H 1/2 (Γ) . (12)
Since the function S m satisfies conditions (3), we introduce the function F defined by

∀ s ∈ R, F (s) = s 0 S m (z)z dz. Note that F is even, that is F (-s) = F (s). Let J 1 be the function defined on H by ∀u ∈ H, J 1 (u) = J 1 (u e , u c ) = Γ F u e (τ ) -u c (τ ) dτ.
One easily checks that J 1 is a C 1 function on H and

J 1 (u) ≥ 1 2 S L Γ |u e (τ ) -u c (τ )| 2 dτ.
Observe that, for any u ∈ H, the derivative J 1 (u) of J 1 at u is the linear map defined by

∀ h ∈ H, J 1 (u).h = Γ S m u e (τ ) -u c (τ ) (u e (τ ) -u c (τ ))(h e (τ ) -h c (τ ))dτ.
Let J be defined by

∀u ∈ H, J(u) = 1 2 Λ σ u, u + J 1 (u) -G, u .
J is of class C 1 on H and J is given by

∀ u ∈ H, J (u) = Λ e u e + S m (u e -u c )(u e -u c ) -G e Λ c u c -S m (u e -u c )(u e -u c ) -G c . (13)
In order to show that J is a monotone operator, we define the nonlinear operator B from H 2 into R by

∀(u, v) ∈ H 2 , B(u, v) = S m ([u])[u] 2 1 - [v] [u] 1 - S m ([v]) S m ([u]) [v] [u] , (14)
where for simplicity we denote by [u] = u e -u c for u = (u e , u c ) ∈ H. Taking into account the fact that S m satisfies (3), one checks easily that B(u, v) ≥ 0. According to ( 12) and ( 13), for (u, v) ∈ H 2 we have

J (u) -J (v), u -v = Λ σ (u -v), u -v + Γ B u(s), v(s) ds, (15)
from which we infer that J is a monotone operator. Therefore J is convex. In order to see the strict convexity of J, that is the strict monotonicity of J , we have to show that if for a given

u, v ∈ H we have J (u) -J (v), u -v = 0, then we have u = v. Observe first that, since λ → λS m (λ) is increasing, we have that B(u, v) = 0 =⇒ [u] = [v].
In particular, for these u, v, we have

B(u, v) = 0, which implies [u] = [v], hence u e -v e = u c -v c .
On the other hand, since

Λ e (u e -v e ), u e -v e = Λ c (u c -v c ), u c -v c = 0,
and since Λ e is coercive, we conclude that u e -v e = 0, which in turn implies that u c -v c = 0.

Finally, this shows that

J (u) -J (v), u -v = 0 =⇒ u = v.
Therefore, J is strictly monotone and J is strictly convex.

In order to show the coerciveness of J, that is J(u) → +∞ when u H → ∞, we proceed as follows. Observe first that

| G c , M (u c ) | ≤ G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) + G c H -1/2 (Γ) M (u e ) H 1/2 (Γ) , ≤ G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) + G c H -1/2 (Γ) u e H 1/2 (Γ) , hence | G, u | ≤ G e H -1/2 (Γ) + G c H -1/2 (Γ) u e H 1/2 (Γ) + G c H -1/2 (Γ) u c -M (u c ) H 1/2 (Γ) + G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) . Using Young's inequality (ab ≤ εa 2 + C(ε)b 2 , for ε > 0 and C(ε) := (4ε) -1 ) we infer | G, u | ≤ ε u e 2 H 1/2 (Γ) + ε u c -M (u c ) 2 H 1/2 (Γ) + 2 C(ε) G 2 H + G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) . (16)
On the other hand, using the fact that 2F (s) ≥ S L s 2 for s ∈ R, we deduce that

J 1 (u) ≥ 1 2 S L Γ |u e (τ ) -u c (τ )| 2 dτ ≥ 1 2 S L Γ |(u e -u c ) -M (u e -u c )| 2 dτ + Γ |M (u e -u c )| 2 dτ .
Using this inequality, together with ( 16), we obtain a lower bound for J(u) (here a(ε), C(ε) are positive constants depending on the arbitrary ε > 0, and b > 0 is a constant):

J(u) ≥ a(ε) u e 2 H 1/2 (Γ) + u c -M (u c ) 2 H 1/2 (Γ) + b Γ |M (u e -u c )| 2 dτ -G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) -C(ε) G 2 H . ( 17
) Inria Since M (u e -u c ) is a constant, we observe that one has M (u e -u c ) H 1/2 (Γ) = c * M (u e - u c ) L 2 (Γ)
for some positive constant c * independent of u. Consequentely, for any ε > 0 so that b -ε > 0, there exists a constant c(ε

) such that b Γ |M (u e -u c )| 2 dτ -G c H -1/2 (Γ) M (u e -u c ) H 1/2 (Γ) ≥ (b -ε) M (u e -u c ) 2 H 1/2 (Γ) -c(ε) G c 2 H -1/2 (Γ) .
Using this inequality in the lower bound [START_REF] Pavlin | Electroporation of planar lipid bilayers and membranes[END_REF], we conclude that

lim u H →+∞ J(u) = +∞
hence J achieves its minimum at a unique point u 0 ∈ H, which satisfies equation [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF].

The dynamical model

In this section, we focus on the dynamical description of the electropermeabilization. Our model is based on the description of two quantities: the time-dependent electric potential and the ratio of the electropermeabilized region over the total membrane area, which is also time-dependent. In subsection 3.1, we present the main considerations that lead to our model. We then study its solvability: existence and uniqueness results are presented in subsection 3.4.

Heuristics of the modeling

Experimental observations suggest that the permeabilization process at a certain location depends on whether the membrane conductivity is above a certain threshold or not. This leads us to define the surface membrane conductivity as an interpolation between the two values S ir and S L , the interpolation parameter ξ(t, s) ∈ [0, 1] being itself a function of time and of the point s on the membrane Γ. In our interpretation, the parameter ξ(t, s) measures in some way the likelihood that a given infinitesimal portion of the membrane is going to be electropermeabilized. More precisely, when ξ(t, s) equals 0 at a given point s ∈ Γ, the membrane conductivity equals the lipid conductivity at this point (thus there is no electropermeabilization), while for ξ(t, s) = 1 it corresponds to the maximal value of the membrane surface conductivity above which electropermeabilization is irreversible. Thus the time-dependent membrane conductivity, denoted by S m writes

∀(t, s) ∈ (0, +∞) × Γ, S m (t, s) = S L + ξ(t, s)(S ir -S L ). (18)
On the other hand, the changes in the conductivity at a certain location s ∈ Γ depend on the transmembrane voltage. Denoting by [u] := u e -u c the jump in the potential between the outside and the inside of the cell for u := (u e , u c ) ∈ H (as we did in the previous sections), we therefore assume that

ξ(t, s) = X(t, [u(t, s)]), (19)
where the function (t, λ) → X(t, λ) will be defined below.

The main idea of the modeling consists in writing a differential equation that describes the dynamics of (t, λ) → X(t, λ) similarly to a sliding door model. Let β be a function satisfying

β ∈ W 1,∞ (R), λ → β(λ) is even on R, λ → λβ (λ) belongs to L ∞ (R), 0 ≤ β(λ) ≤ 1, β is non decreasing on (0, +∞), lim λ→+∞ β(λ) = 1.            (20)

RR n°8005

An example of such a function would be

∀λ ∈ R, β(λ) := (1 + tanh(k ep (|λ| -V rev ))/2.
Let us describe now the evolution of X. For λ 0 ∈ R, which will stand for a given ∆TMP, and for an initial value X 0 ∈ [0, 1] of X, set β 0 := β(λ 0 ). Then we consider the following two possibilities:

• Either β 0 -X 0 is positive, in which case the electric pulse is sufficiently high to enlarge the electropermeabilized region, with a characteristic time of electropermeabilization of order τ ep .

• Or β 0 -X 0 is negative, in which case we consider that the pulse is not high enough to increase the electropermeabilization. Therefore, the membrane tries to reseal with a characteristic resealing time of order τ res . Since experimental observations suggest that this phenomenon takes much more time than the electropermeabilization process, we assume that τ res > τ ep .

Remark that when a cell is at rest, X 0 equals zero, but if high voltage pulses have been applied earlier than the initial time, X 0 might not be equal to zero.

Based on these considerations, for any λ 0 ∈ R, we assume that X(•, λ 0 ) satisfies the following differential equation:

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

Statement of the mathematical problem

We first write the equation satisfied by the potential U defined on the domain O e ∪ O c . We assume that before the imposition of the electrical pulses g on the external boundary ∂Ω, the cell potential is at rest and given by U 0 ∈ H 1 (Ω). This resting potential translates the ionic exchanges through the membrane. According to equalities ( 18)-( 19), the membrane conductivity S m is defined by: S m (t, λ)

:= S L + (S ir -S L )X(t, λ). (22)
We seek the solution (U, X) to the following system of equations: for (U 0 , X 0 ) given, the pair of functions (U, X) satisfy:

U | t=0 = U 0 , and for any t > 0,

∆U = 0, in (0, T ) × (O e ∪ O c ) , U (t, •) = g(t, •) on (0, +∞) × ∂Ω, (23a) [σ∂ n U ] = 0, on (0, T ) × Γ, (23b)
C m ∂ t [U ](t, •) + S m (t, [U ])[U ] = σ c ∂ n U (t, •) | Γ -, on (0, T ) × Γ, (23c)
where, writing λ = [U ](t, s), the function X appearing in [START_REF] Teissié | In Vitro cell electropermeabilization[END_REF] satisfies the differential equation

∂X(t, λ) ∂t = max β(λ) -X(t, λ) τ ep ; β(λ) -X(t, λ) τ res , t > 0, (24a)
X(0, λ) = X 0 . ( 24b
) Inria Remark 3.1.
Observe that if the source g does not depend on the time t, the stationary point (U * , X * ) of the system (23)-(24a) is the unique solution to

∆U * = 0, in (O e ∪ O c ) , U * | ∂Ω = g on ∂Ω, [σ∂ n U * ] = 0, on Γ, (S L + (S ir -S L )X * ) [U * ] = σ c ∂ n U * | Γ -, on Γ,
where

X * = β([U * (s)]),
which coincides with the static model (2).

Properties of the function X

Let us state the following lemma regarding the solution to equation (24):

Lemma 3.2. Let T > 0 be fixed and let β satisfy condition [START_REF] Serša | Application of electroporation in electrochemotherapy of tumors[END_REF]. For T > 0 and any λ ∈ R, the following differential equation

   ∂X(t, λ) ∂t = max β(λ) -X(t, λ) τ ep ; β(λ) -X(t, λ) τ res , ∀t ∈ (0, T ), X(0, λ) = X 0 ∈ [0, 1], (25)
has a unique solution X(•, λ) ∈ C 1 ([0, T ]). Moreover one has 0 ≤ X(t, λ) = X(t, -λ) = X(t, |λ|) ≤ 1, ∀t ∈ [0, T ], ∀λ ∈ R.
In addition, there exists a constant K(T ) > 0 such that for any λ 1 , λ 2 ∈ R we have

∀t ∈ [0, T ], |λ 1 X(t, λ 1 ) -λ 2 X(t, λ 2 )| ≤ K(T ) |λ 1 -λ 2 | . (26)
Proof. The mapping

X → max([β(λ) -X]/τ ep , [β(λ) -X]/τ res ), defined from R into itself, is
clearly Lipschitz for all fixed λ. Therefore, the differential equation (25) has a unique solution X ∈ C 1 ([0, T ]), for any given X 0 ∈ R. Assuming now that 0 ≤ X 0 ≤ 1, multiplying the equation (25) by X -:= max(-X, 0), and using the fact that β ≥ 0, one gets

1 2 ∂ t |X -| 2 ≤ 0. Thus X -(•, λ) ≡ 0 on [0, T ].
Similarly, using the fact that β -1 ≤ 0, and multiplying the equation by (X -1) + = max(X -1, 0), one sees that

∂ t |(X -1) + | 2 ≤ 0,
and finally 0 ≤ X(t, λ) ≤ 1, for any t ∈ [0, T ].

In order to show the estimate (26), we proceed as follows. For λ ∈ R denote by Y (t, λ) := λX(t, λ). Since for any (a, b) ∈ R, we have max(a, b) = a + (a -b) -, assuming for instance that τ ep < τ res we observe that

max(a/τ ep , a/τ res ) = a τ ep + τ -1 ep -τ -1 res a -,
and using the equality a -= (|a| -a) /2 we get

max(a/τ ep , a/τ res ) = a τ -1 ep + τ -1 res /2 + |a| τ -1 ep -τ -1 res /2.
From this we induce that Y satisfies the following O.D.E

   ∂Y (t, •) ∂t = H(λ, Y ), ∀t ∈ (0, T ), Y (0, λ) = Y 0 = λX 0 , (27)
where the function H is defined by

H : (λ, Y ) → H(λ, Y ) = (λβ(λ) -Y ) τ -1 ep + τ -1 res /2 + λ |β(λ) -Y /λ| τ -1 ep -τ -1 res /2. ( 28
)
In order to see that the mapping λ → H(λ, Y ) is globally Lipschitz on R, since by our assumptions on the function β we know that λβ (λ) is uniformly bounded on R, the function λ → λβ(λ) is Lipschitz on R, and thus we have only to verify that the function h

: λ → h(λ, Y ) = λ |β(λ) -Y /λ| is Lipschitz on R.
Indeed, this is the case as one can see by a simple computation that

∂h ∂λ = |β(λ) -Y /λ| + (λβ (λ) + Y /λ) β(λ) -Y /λ |β(λ) -Y /λ| , = (β(λ) + λβ (λ)) β(λ) -Y /λ |β(λ) -Y /λ| , and therefore ∂h ∂λ (λ, Y ) ≤ |β(λ) + λβ (λ)| .
From this, one infers clearly that the mapping (λ,

Y ) → H(λ, Y ) is globally Lipschitz on R × R. Writing Y (t, λ 1 ) -Y (t, λ 2 ) = t 0 (H(λ 1 , Y 1 ) -H(λ 2 , Y 2 )) (t)dt + (λ 1 -λ 2 )X 0 ;
we conclude that ∀t ∈ (0, T ),

|Y (t, λ 1 ) -Y (t, λ 2 )| ≤ K(T ) |λ 1 -λ 2 | ,
thanks to an invocation of Gronwall lemma.

Existence and uniqueness of the dynamical potential

Since the non-linearity of problem [START_REF] Teissié | Mechanisms of cell membrane electropermeabilization: A minireview of our present (lack of?) knowledge[END_REF] appears in the transmission condition (23c), we are going to rewrite it on the surface Γ using the Steklov-Poincaré operators, in the same manner as we did in the previous section for the static model. We first prove the following property:

Lemma 3.3. The operator Λ e + Λ c is positive, selfadjoint and invertible from H 1/2 (Γ) into H -1/2 (Γ). The operator B := Id +Λ -1 e Λ c is therefore invertible, from H 1/2 (Γ) into itself.
Moreover, define the domain D(Λ c B -1 ) as

D(Λ c B -1 ) = ϕ ∈ H 1/2 (Γ) : Λ c B -1 ϕ ∈ L 2 (Γ) .
The operator

Λ c B -1 , D(Λ c B -1 ) is m-accretive (more precisely D(Λ c B -1 ) = H 1 (Γ)).

Inria

Proof. That Λ e + Λ c is invertible is an easy consequence of the fact that the operator Λ e is a positive, selfadjoint and invertible operator while Λ c is non-negative and selfadjoint. Thus, Λ e + Λ c is also a positive selfadjoint invertible operator. Let ϕ ∈ D(Λ c B -1 ). Then, by definition (and invertibility) of B (i.e. Id = Id +Λ -1 e Λ c B -1 ) we have

Λ c B -1 ϕ, ϕ = Λ c B -1 ϕ, ϕ , = Λ c B -1 ϕ, B -1 ϕ + Λ c B -1 ϕ, Λ -1 e Λ c B -1 ϕ , = Λ c B -1 ϕ, B -1 ϕ + Λ e Λ -1 e Λ c B -1 ϕ, Λ -1 e Λ c B -1 ϕ ≥ 0. Λ c B -1 is therefore accretive. Let f ∈ L 2 (Γ)
, let λ > 0 and let U be the unique solution in P H 1 0 (Ω) to the following problem:

-∆U = 0, in O e ∪ O c , U | ∂Ω = 0, σ e ∂ n U | Γ + = σ c ∂ n U | Γ -, λσ c ∂ n U | Γ -+ U | Γ --U | Γ + = f. Therefore, setting v := U | Γ --U | Γ + , v satisfies v + λΛ c B -1 v = f.
In addition, since

Λ c B -1 is nonnegative , we have v L 2 (Γ) ≤ f L 2 (Γ) . Therefore, Λ c B -1 is m-accretive.
Lemma 3.4. Let U 0 and X 0 be two enough regular functions defined respectively in Ω and on Γ, and recall that S m (t, λ) := S L + (S ir -S L )X(t, λ) is defined in [START_REF] Teissié | In Vitro cell electropermeabilization[END_REF]. Finding the solution (U, X) to problem (23)-( 24), if it exists, is equivalent to finding (u e , u c , X), with u e = U |Γ + and u c = U |Γ -satisfying:

u e = u c -v, (29)
u c = B -1 v -Λ -1 e Λ 0 g , ( 30
)
where v is the solution to

C m ∂ t v + Λ c B -1 v + S m (t, v)v = G, v(0, •) = ϕ, (31)
with ϕ and G being defined as

ϕ = U 0| Γ + -U 0| Γ -, G := Λ c B -1 Λ -1 e Λ 0 g,
and where, writing λ = v(t, s), a.e.(t, s)

∈ (0, T ) × Γ, X satisfies    ∂ t X(t, λ) = max β(λ) -X(t, λ) τ ep ; β(λ) -X(t, λ) τ res , ∀t > 0, X(0, λ) = X 0 . (32)

RR n°8005

Proof. The lemma is a straightforward consequence of the definition of the Steklov-Poincaré operators Λ c , Λ e , Λ 0 , and of the invertibility of Λ e . Indeed, condition (23b), that is the continuity of the flux across Γ, boils down to Λ e u e + Λ 0 g + Λ c u c = 0, from which, thanks to the invertibility of Λ e , we infer that

u e -u c = -Bu c + Λ -1 e Λ 0 g = -v.
In addition, by definition of v, we have

Λ c u c = Λ c B -1 v -Λ c B -1 Λ -1 e Λ 0 g = Λ c B -1 v -G.
Then, the transmission condition (23c) (multiplied by (-1)) reads (31), provided we recall that β is an even function and that S m is defined by [START_REF] Teissié | In Vitro cell electropermeabilization[END_REF].

We now show that the evolution equations appearing in lemma 3.4 have a unique solution.

Theorem 3.5. Assume that β satisfies (20), G ∈ L p ((0, T ); L 2 (Γ)) for some p > 1, and that ϕ ∈ L 2 (Γ) is given. Let X 0 ∈ L ∞ (Γ) such that 0 ≤ X 0 ≤ 1 on Γ, and let S m be defined as in [START_REF] Teissié | In Vitro cell electropermeabilization[END_REF].

Then, there exists a unique function v ∈ C([0, T ]; L 2 (Γ)), mild solution to the system

C m ∂ t v + Λ c B -1 v + S m (t, v)v = G, ∀t ∈ (0, T ), v(0) = ϕ, (33)
where S m is given by [START_REF] Teissié | In Vitro cell electropermeabilization[END_REF], and writing λ = v(t, s), a.e (t, s) ∈ (0, T ) × Γ, one has

   ∂ t X(t, λ) = max β(λ) -X(t, λ) τ ep ; β(λ) -X(t, λ) τ res , ∀t ∈ (0, T ), X(0, λ) = X 0 . (34)
Moreover, if ϕ ∈ H 1 (Γ) and G ∈ W 1,1 ((0, T ); L 2 (Γ)), the above mild solution is a classical solution to (33), in the sense that

v ∈ C([0, T ]; H 1 (Γ)) ∩ C 1 ([0, T ]; L 2 (Γ)).
Proof. To simplify the notations we denote by A the operator

A = 1 C m Λ c B -1 .
In a first step, we are going to show existence and uniqueness of a mild solution to (33) in C([0, T ]; L 2 (Γ)). Thanks to lemma 3.2, we know that for any λ ∈ R, the solution to

   ∂ t X(t, λ) = max β(λ) -X(t, λ) τ ep ; β(λ) -X(t, λ) τ res , ∀t > 0, X| t=0 = X 0 ,
Inria exists and belongs to C 1 ([0, T ]). Moreover, 0 ≤ X ≤ 1, and the mappings λ → X and λ → λX are Lipschitz. Now, upon setting

F(t, v) := -S m (t, v)v + G = -(S L + (S ir -S L )X(t, v)) v + G,
solving equations ( 33) is equivalent to finding v ∈ C([0, T ]; L 2 (Γ)) solution to the following equation, which is the mild version of equation ( 33):

v = e -tA ϕ + 1 C m t 0 exp (-(t -τ )A) F(τ, v)dτ. (35)
Thanks to inequality (26) of lemma 3.2, it is clear that the mapping

v → F(•, v),
is Lipschitz from the space E := C([0, T ]; L 2 (Γ)) into itself. This means that there exists K 1 > 0 such that for any v, w ∈ C([0, T ]; L 2 (Γ)) we have

F(•, v) -F(•, w) L ∞ ((0,T );L 2 (Γ)) ≤ K 1 v -w L ∞ ((0,T );L 2 (Γ)) .
We shall endow the space E with the norm

ψ E := sup t∈[0,T ] e -αt ψ(t, •) L 2 (Γ) ,
for some α > 0, which will be chosen below. If we set

Φ(U c )(t) := e -tA ϕ + 1 C m t 0 exp (-(t -τ )A) F(τ, U c (τ ))dτ,
then Φ : E -→ E is a continuous mapping. We shall check that, upon choosing α appropriately, it is a strict contraction. Thus, it has a unique fixed point, providing the unique solution of (35). Indeed

Φ(U c )(t) -Φ(V )(t) = t 0 exp - t -τ C m Λ c B -1 [F(τ, U c ) -F(τ, V )] dτ.
Since the operator Λ c B -1 is m-accretive, the operator A generates a contraction semi-group and the following estimate holds

e -(t-τ )A [F(τ, U c ) -F(τ, V )] L 2 (Γ) ≤ [F(τ, U c ) -F(τ, V )] L 2 (Γ) .
Therefore, we infer

Φ(U c )(t) -Φ(V )(t) L 2 (Γ) ≤ t 0 F(τ, U c ) -F(τ, V ) L 2 (Γ) dτ, ≤ K 1 U c -V E t 0 e ατ dτ, ≤ C(p)K 1 α -1 e αt U c -V E ,
from which we conclude that

Φ(U c ) -Φ(V ) E ≤ α -1 C(p)K 1 U c -V E .

RR n°8005

This implies that for α large enough, the mapping Φ is a strict contraction on E. Equation (33) has thus a unique mild solution in C([0, T ]; L 2 (Γ)).

Suppose now that G belongs to W 1,1 ((0, T ); L 2 (Γ)) and that ϕ ∈ H 1 (Γ). The mild solution given by formula (35) belongs to C 1 ([0, T ]; L 2 (Γ)) hence we infer

∂ t v + S m (•, v)v ∈ C([0, T ]; L 2 (Γ)).
Therefore, setting U c = U |Oc , we have

∆U c (t, •) = 0, in O c , σ c ∂ n U c (t, •) ∈ L 2 (Γ),
from which we infer that U c (t, •) ∈ H 3/2 (O c ). Hence we deduce that since the domains are smooth U c (t, •) |Γ -belongs to H 1 (Γ). Similar reasoning for U e = U |Oe implies that U e (t, •) |Γ + belongs to H 1 (Γ), and therefore the jump v belongs to C([0, T ]; H 1 (Γ)), which ends the proof.

Numerical simulations

In this section, we provide some numerical results that show the consistency of our models. In order to solve both static and dynamical problems, we first present the finite-difference method on a cartesian grid adapted from the second-order scheme of Cisternino and Weynans [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF]. Note that it is not the scope of this paper to prove rigorously that our scheme is of order two, however numerical simulations of Fig 3(b) seems to confirm this order. The rigorous numerical analysis of the method will be performed in a forthcoming work. Since the long-term goal of this work is to fit the models with the experimental data, and since experiments are performed with several thousands of cells per cm 2 , parallel computing will be necessary to solve the inverse problem. Therefore, we choose to use accurate schemes on cartesian grids, such as the scheme of [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF], for which the accuracy has been shown and the parallelization has be already performed, despite Galerkin formulations might be more adapted to the single cell problem, since variational formulation holds.

Spatial discretization

We perform the discretization on a cartesian grid covering the domain Ω = O e ∪ O c , which is a square domain of length L (L equals 200 µm for the computations, see Table 1). The interface is described by a level-set function [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF], which separates the extra-and intra-cellular domains by the use of a signed distance function ϕ. The normal to the interface n(x) outwardly directed from the inner to the outer of the cell is directly obtained by computing numerically ∇ϕ(x).

The grid spacing is denoted by h, and N is the number of points such as

N = L/h.
For any (i, j) ∈ N 2 we denote by M ij the grid points defined by

M i,j = (x i , y j ), where x i = ih x , y j = jh y , ∀(i, j) ∈ N 2 .
The numerical approximation of the solution to the static or to the dynamical model at the point (x i , y j ) is generically denoted by u ij . Standard approximation of the Laplacian is used in the discretized domains O e and O c far from the interface, which is the cell membrane Γ. However, due to the jump conditions, a special treatment of the approximation of the Laplacian and of the computation of the fluxes is needed at the points nearing the cell membrane.

Inria

If the intersection of the interface and [M ij M i+1j ] exists, then we define the interface point I i+1/2, j = (x i+1/2 ,j , y j ) as this intersection. We create two additional unknowns at this interface point, called interface unknowns, and denoted by ũe i+1/2 ,j and ũc i+1/2 ,j . The interface point I i, j+1/2 = (x i , ỹi, j+1/2 ) is similarly defined as the intersection of Γ and the segment

[M ij M ij+1 ].
An example of the discretization method is given by Figure 2. On regular grid points, that are not neighboring the interface, the Laplacian is discretized with a standard centered second-order finite-difference scheme. A specific five points stencil including the interface points is used for neighboring points, as shown in Figure 2(a).

Figure 2(b) provides an example of the discretization of ∇U on both sides of the interface. The x-derivative of U can be computed with second-order accuracy using a one-sided formula involving three grid points. For example we approximate the flux on the left (for instance exterior) side of the interface with the points M i-1j , M ij and I i+1/2 ,j by:

∂U ∂x (x, y j ) ≈ (u i-1j -ũe i+1/2 ,j )(x i -x) h x (x i-1 -x) - (u ij -ũe i+1/2 ,j )(x i-1 -x) h x (x i -x) , (36)
where for the sake of brevity, we have replaced xi+1/2 ,j by x. The y-derivative cannot be obtained in the same way, since there are no grid points aligned with the interface point in the y-direction. 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

∂U e ∂y (x, y j ) ≈ x -x i-1 h x (∂ y u) ij - x -x i h x (∂ y u) i-1j . (37)
The formulas for (∂ y u) ij and (∂ y u) i-1j depend on the local configuration on the interface, but they are based on the same principle as for (36). The scheme is stabilized by using a shifted y-stencil if two interface points are involved in the same flux discretization, as illustrated by Figure 2(c).

j j+1 j-1 i i-1 i+1 i+2 I i+1/2,j (a)

Accuracy of the finite difference method

In order to show the accuracy of the numerical method, we compare both explicit and numerical solutions to the linear static problem (4), i.e without electroporation. This means that in (4c), S m is constant equal to S L .

Consider a domain composed of two concentric disks. The cell is the disk of radius R 1 , centered at 0. The domain Ω is the disk with the same center as the cell, and whose radius R 2 is stricly greater than R 1 . The boundary data g equals ER 2 cos θ, such that the cell is embbeded in a uniform electric field of magnitude E in the x-direction. The exact solution U to (4) is then explicitely given by

∀(r, θ) ∈ (R 1 , R 2 ) × [0, 2π],
U e (r, θ) = (α e r + β e r -1 ) cos θ,

∀(r, θ) ∈ (0, R 1 ) × [0, 2π], U c (r, θ) = α c r cos θ, (38)
where α e , β e , and α c are 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 analytic solution is projected on the edges of the square Q L centered at 0, and whose characteristic length

L satisfies R 1 < L/2 < L √ 2/2 < R 2 (see Fig 3(a)
). In order to verify the accuracy of the space discretization, at least in this configuration, we solve numerically problem (4) in the square Q L , with the trace U | ∂Q L of the analytic solution U on the edges of Q L as Dirichlet boundary condition for different grid spacings h. We set the parameters equal to R 2 = 150 µm, R 1 = 50 µm, L = 200 µm, E = 400 V/cm, the electric parameters being given by Table 1. These parameters come from the values given by DeBruin and Krassowska, see Table 1 of [START_REF] Debruin | Modelling electroporation in a single cell. I. Effects of field strength and rest potential[END_REF].

We then compare the numerical solution U h to U . The relative error is computed using both grid and interface points:

E( U , U h ) := ||U h -U || L 2 (Ω) + |U h -U | L 2 (Γ) || U || L 2 (Ω) + | U | L 2 (Γ) . (39)
We denote by h r the relative grid spacing defined by

h r = h/L. Figure 3(b)
shows the behavior of E( U , U h ) with respect to h r : the space discretization is of order two for this specific case. Note Cisternino and Weynans [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] have shown the second order accuracy of their method, therefore despite we have adapted the scheme to our problem, we are confident in the good accuracy of the numerical method.

Computation of the non-linear static model

When solving the static equation ( 4), if one uses the following naive iterative 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. Relative error between the two solutions with respect to the grid spacing is plotted in Figure 3(b). one finds that the iterative scheme oscillates between two values. This might be a consequence of the fact that the membrane conductivity takes its extreme values S L and S ir instead of reaching an intermediate state. Another issue might be that the mapping L g defined by

S m ([U n ])[U n+1 ] = σ c ∂ n U n+1
L g : v → u, where u ∈ P H 1 (Ω) is solution to ∆u = 0, in O e ∪ O c , u | ∂Ω = g, [σ∂ n u] Γ = 0, S m [v] Γ [u] Γ = σ c ∂ n u| Γ ,
is not a contraction, since the Lipschitz constant of S m is of order k ep S m 1. We use a modified mapping L ρ,g :

L ρ,g : v → u, where u ∈ P H 1 (Ω) is solution to ∆u = 0, in O e ∪ O c , u | ∂Ω = g, [σ∂ n u] Γ = 0, [u] Γ + ρS m [v] Γ [u] Γ -ρσ c ∂ n u| Γ = [v] Γ ,
where ρ is a small positive parameter chosen so as L ρ,g be a contractive operator. The following iteration process is used: U 0 given, and for any n ≥ 0 U n+1 = L ρ,g (U n ).

For the simulations we set the stopping criterion of the scheme at γ = 10 -12 , meaning that the numerical solution is obtained when the relative error E(U n+1 , U n ) is smaller that γ.

Influence of the parameter S ir

The parameter S ir , which is the conductivity of the fully "electroporated" membrane is hardly measurable by the experiments. It is therefore important to investigate its influence on the model. Figure 4 shows that S ir has a little influence on the membrane conductivity S m , as the value of X counter-balance the variation of S ir . Therefore, the numerical criterion to define the electroporation should involve the parameter S m .

Comparison with the model of Ivorra, Mir and Villemejeane

We compare our results with the simulations of Ivorra et al. by studying the influence of the extracellular medium conductivity on the membrane conductivity. In Figure 5, we show results similar to those presented in Figure 7 of [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF]. Note that to perform their simulations, Ivorra et al. have multiplied by ten the membrane thickness. This is the reason why their permeabilizing field is of order of magnitude of 2MV/m, which is much higher than the magnitude used in the experiments, typically [START_REF] Mir | Therapeutic perspectives of in vivo cell electropermeabilization[END_REF][START_REF] Mir | Electroporation of cells in tissues. Methods for detecting cell electropermeabilisation in vivo[END_REF][START_REF] Gowrishankar | Microdosimetry for conventional and supra-electroporation in cells with organelles[END_REF] to the range 20 to 30 kV/m. In contrast, our model provides more realistic conditions of electropermeabilization.

Computation of the dynamical problem

Time-discretization of the model

The time-derivative ∂ t [U ] of ( 23) is discretized using the following scheme:

C m [U ] n+1 -[U ] n dt -σ c ∂ n U n+1 c + S m (t n , [U ] n )[U ] n = 0. (40)
The Runge-Kutta method of order 4 is used to compute the variable X, with time steps dt. show the numerical results at t = 100 µs using the parameters of Table 1. In order to visualize the membrane electropermeabilization, we depict it with boxes, which are colored and sized according to the values of S m at each point of Γ and at t = 100 µs. We emphasize this is a visualization artefact: in our model, the cell membrane is a surface without any thickness. 1. Fig. 6(b) shows that the electropermeabilized regions depend on the shape and orientation of Γ.

Main parameters influence

The key parameters of the model define the electropermeabilization coefficient S m , that is k ep , V rev , S ir and the characteristic times τ ep and τ res . A numerical sensitivity analysis was led to determine how the behavior of the solution with respect to a variation of each specific parameter, as shown on Figure 7. All the parameters defining the function β have a very small influence on the average X of X over the cell membrane. Even for small values of k ep , the values of X are only modified by a factor 2 (Figure 7(a)). On the other hand, the "fully electroporated" membrane conductivity S ir , which was first taken as (σ c + σ e )/(2δ), affects greatly the order of magnitude of X, changing from 10 -6 to 10 -2 (Figure 7(d)). Therefore, as for the static case the relevant quantity to observe the phenomenon is S m .

Comparison with the model of Neu, Krassowska, and Debruin

The main difference between the model of Neu, Krassowska, et al. and ours resides in the addition of an electroporation current I ep = N ep i ep , instead of a direct description of the variations of the surface membrane conductivity S m .

The equation satisfied by the transmembrane voltage in the model of Neu, Krassowska, et al. [START_REF] Debruin | Modelling electroporation in a single cell. I. Effects of field strength and rest potential[END_REF] reads

-σ c ∂ n U c = C m ∂ t [U ] + S L [U ] + N ep i ep , (41)
where the ionic reversal currents have been neglected. N ep is the pore density, obeying the 1

E = 10kV/m 1 E = 25kV/m 1 E = 40kV/m E = 10kV/m E = 25kV/m E = 40kV/m (a) Influence of EP switch coefficient kep. V rev (V) S m (S/cm 2 ) 0 0.5 1 1.5 2 2.5 3
1.25 2.5

3.75

5 1 E = 10kV/m 1 E = 25kV/m 1 E = 40kV/m E = 10kV/m E = 25kV/m E = 40kV/m (

Long-time behavior of the numerical solution to the dynamical model

In this paragraph, we compare the long-time behavior of the solution U dyn to the dynamic model for a constant pulse with the solution U stat to the static model.

Simulations are done in order to reach the time scale of the resealing characteristic time τ res (Figure 10). A constant pulse g is applied until the steady state of the dynamical system is reached. Observe that the stationary value of S m , which is reached after about 500 µs is lower than the value of S m a few tens of microseconds after the begining of the pulse delivery. For pulses of short duration, it is therefore important to simulate our new dynamical model instead of the steady model of Ivorra et al. [START_REF] Ivorra | Electrical modeling of the influence of medium conductivity on electroporation[END_REF].

Conclusion

In this paper, we introduce two models describing the electropermeabilization of a single cell. We first study the static model, which is inspired by Ivorra Inria uniqueness results. We then derive a new dynamic model of cell electropermeabilization, which takes the permeabilizing time into account. We studied mathematical properties of this new model. We then provided an accurate finite-difference method on cartesian grid to compute these models, and we eventually presented numerical simulations for both static and dynamic models, that corroborate the results of the most achieved model of Neu, Krassowska, et al.

The main feature of our models lies in the fact that without loss of accuracy it is composed by a small number of parameters (mainly 4 parameters: S ir , V rev , and τ ep , and τ res for the dynamical system) compared with the sophisticated models with tens of hardly measurable parameters of Neu, Krassowska, et al. Therefore, a forthcoming fitting of our models with the experimental data seems feasible, which is hardly the case for models with a large number of parameters.

From the biological point of view, we highlight the fact that the static model can be used for very long pulses (around 1 ms) but, for short pulses around 10 µs and below, the dynamics of the phenomenon have to be considered.

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 .

Λ

  e (f ) := n e • σ e ∇v e | Γ + , where div(σ e ∇v e ) = 0 in O e , v e | ∂Ω = 0 and v e | Γ = f.

  Steklov-Poincaré operator (or what is sometimes called a Neumann-to-Dirichlet operator). Namely, for any ψ ∈ H -1/2 (Γ) given, one has Λ -1 e (ψ) = v |Γ where v ∈ H 1 (O e ) satisfies the equation div(σ e ∇v) = 0 in O e , v |∂Ω = 0 and n e • σ e ∇v |Γ + = ψ.

  Discretization of the Laplacian on the points at the interface. Discretization of ∇U at the interface: non-stabilized stencil. Discretization of ∇U at the interface: stabilized stencil.

Figure 2 :

 2 Figure 2: Discretizations of the Laplacian and of the gradient of U at the interface. The first y-derivative stencil on the right side is shifted to avoid an ill-conditioned discretization.

  Domains of calculation of the analytic and of the numerical solutions.h (log scale)E(u -U ) (log scale)

  Log-log diagram of the error E( U , U h ) given by (39) with respect to hr.

Figure 3 :

 3 Figure3: 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. Relative error between the two solutions with respect to the grid spacing is plotted in Figure3(b).

Figure 4 :

 4 Figure 4: Membrane surface conductivity S m at the pole of the cell (θ = 0) obtained by solving the static equation with different values of S ir for pulse magnitudes from 10 to 50 kV/m (the numerical values of the other parameters are given in Table1).

Figure 5 :

 5 Figure 5: Non-linear membrane conductivity of the static model for 4 different extracellular conductivity: σ e = 1 S/m (•), σ e = 0.1 S/m (•), σ e = 0.01 S/m ( ), σ e = 0.001 S/m ( ).

Figure 6 (

 6 Figure6(a) and 6(b) show the numerical results at t = 100 µs using the parameters of Table1. In order to visualize the membrane electropermeabilization, we depict it with boxes, which are colored and sized according to the values of S m at each point of Γ and at t = 100 µs. We emphasize this is a visualization artefact: in our model, the cell membrane is a surface without any thickness.

  (a) Solution to the dynamical model for a circular cell. (b) Solution to the dynamical model for a smooth non-convex cell shape.

Figure 6 :

 6 Figure 6: Solution to the dynamical problem with two different cell shapes at t = 100 µs. The numerical parameters are given by Table 1. Fig. 6(b) shows that the electropermeabilized regions depend on the shape and orientation of Γ.

1 E 1 E

 11 b) Influence of threshold voltage Vrev. τ ep (log) = 25kV/m 1 E = 40kV/m E = 10kV/m E = 25kV/m E = 40kV/m (c) Influence of the characteristic time τep. S ir (S/cm 2 ) (log) = 10kV/m 1 E = 25kV/m 1 E = 40kV/m E = 10kV/m E = 25kV/m E = 40kV/m (d) Influence of the EPd surface conductivity S ir .

Figure 7 :

 7 Figure 7: Influence of each parameter on the mean value S m of S m (Fig. 7(a)-7(b)-7(c)), and on the mean value X of X (Fig. 7(d)) at t = 100 µs. Three magnitudes of electric pulses are considered: 10, 25 and 40 kV/m.

Figure 8 :

 8 Figure 8: Comparison of the ∆TMP obtained respectively with our model (solid lines) and with the model of Neu, Krassowska, et al. (dashed), with parameters from table 1.

  et al., and we show existence and Polar angle θ (rad)

Figure 9 :Figure 10 :

 910 Figure 9: Current density through the membrane of the model of Neu, Krassowska, et al., that is J ep = N ep i ep (in dashed line), compared with the membrane current density of our model (in solid line), J ep = S m [U ], along the the cell membrane at 100 µs.

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	σ c	0.455	S/m
	Capacitance	C m	9.5 × 10 -3 F/m
	Membrane surface conductivity	S L	1.9	S/m
	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	τ res	1 × 10 -3	s
	EPd membrane surface conductivity	S ir	2.5 × 10 8	S/m
	Numerical parameters:
	Simulation box size	L	200	µm
	Grid points (each side)	N	50
	Time step	∆t	20	ns
	Pulse duration	T p	100	µs
	Duration of simulation	T f	150	µs
	Number of time steps	N T	7500
				Inria

RR n°8005

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

The authors thank very warmly the biologists of the CNRS research team "UMR 8203 Vectorology and anti-cancerous therapeutics" at the Institut Gustave Roussy. Particulary, the discussions with Dr. A. Silve and Dr. L.M. Mir were very helpful in the understanding of the cell electropermeabilization. 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).