Cell scale modeling of electropermeabilization by periodic pulses.

: In this paper, we focus on the behaviour of periodic solutions to a cell-scale elec-tropermeabilization model previously proposed by Kavian et al. (JMB,2012). Since clinical per-meabilization protocols mostly submit cancer cells to trains of periodic pulses, we investigate on parameters that modify signiﬁcantly the resulting permeabilization. Theoretical results of existence and uniqueness of periodic solutions are presented, for two diﬀerent models of membrane electric conductivity. Numerical simulations were performed to corroborate these results and illustrate the asymptotic convergence to periodic solutions, as well as the dependency on biological parameters such as the cell size and the extracellular conductivity.


Introduction
Electropermeabilization or electroporation is a phenomenon that occurs when a biological cell or a lipid vesicle is submitted to a high electric field.If a sufficiently large potential difference is applied to the membrane, its strucure is altered and molecules that are usually not able to enter the cytoplasm can diffuse inside the cell.Reversible electropermeabilization is already used to improve delivery of drugs such as bleomycin in oncology [1,4], and also to make the cell permeable to very large molecules for gene transfer [7].
This paper is an extent of a previous modeling work [5] in which we proposed a phenomenological approach of electropermeabilization, compared to the previous most achieved models of Neu, Debruin and Krassowska [9,3].
Experiments show that a single pulse is often not sufficient to achieve cell electropermeabilization, with a wide variety of pulses, from nanopulses to millipulses [8,12,11].The applied pulses are usually trains of several pulses, and can be considered as a periodic source.Here we focus theoretically and numerically on the behaviour of solutions to our equations in this particular case of sources.The objective of this paper is to show that the membrane conductivity reaches a stable periodic state when submitted to a periodic voltage.
When very intense electric pulses are clinically applied, healthy tissues in the vicinity of the electrodes are completely damaged due to the high electric field.Moreover, thermal effects induce cell destruction if pulses last for too long.It is therefore important to determine pulse characteristics (voltage, duration) that are sufficient to permeabilize the targeted tissue, minimizing at the same time the caused damage.Considering that, we will focus on the amount of time needed for the cell to reach a periodic conductive state, depending on which model is chosen to describe the membrane conductivity.
We will show, in particular, that using a static model of conductivity results in a high convergence speed of the solution to a periodic state.A first conclusion should be that a few pulses are needed to reach it.However, results given by a fully dynamical model infer that more pulses are required to actually achieve electropermeabilization.
After recalling in the first section the model on the electric potential in a biological cell, and especially across its membrane, we will justify the existence and uniqueness of periodic solutions in the second section.In the final section, we will illustrate how a cell at rest submitted to a periodic source reaches a periodic state.We will validate these results with numerical simulations.

The electric potential in a biological cell
A biological cell is considered as an homogeneous medium O c , which is separated from an exterior bath O e by a phospholipidic membrane Γ (see fig. 1).Due to the high resistivity and the thickness of the membrane, it is considered as a single surface interface between the two domains.The capacity of this electric material is denoted by C m , and its surface conductivity by S m .Let σ be the conductivity of the medium, considered piecewise constant: The electric potential U in both intra-and extracellular domains follows Poisson's law: with the transmission conditions on the flux and the transmembrane potential across the interface Γ: where g represents the potential that is applied by external electrodes, and I ep is the total current due to electropermeabilization.In the mostly used models of electroporation proposed by Neu, Krassowska and De Bruin [9,3], this current was defined as the product of a pore density N ep by the current flowing through a single pore i ep .The pore current was defined by a highly non-linear function of the transmembrane potential.
Derived from a linearization of this total current I ep , we proposed in our previous work [5] to replace this definition by the product of the surface conductivity of the membrane S m by the transmembrane potential difference [U ] Γ .The transmission condition now holds : In order to take into account electropermeabilization, we write S m as the sum of the lipid conductivity at rest S 0 and a non-linear part S ep depending on [U ] Γ : The next section introduces the different conductivity models we chose.

Models of membrane conductivity
Throughout this paper, we will focus on two models of membrane conductivity: • a static conductivity model, which has a sigmoid profile: Inria where S 1 is another conductivity constant, much larger than S 0 .k ep defines the speed of the switch between the rest state (S m = S 0 ) and the fully permeabilized state (S m = S 0 + S 1 ).V rev designates the voltage threshold that needs to be crossed to permeabilize the membrane.
We will refer to this model as the β-model.
• a dynamical model which will be called X-model: for t 0, where β is the same function as the β-model and τ ep is the characteristic time constant of formation and resealing of pores.The X variable follows a sliding-door model around the electropermeabilization threshold V rev .Using this model makes the dynamic of electropermeabilization intrisic to the cell and not dependent of the applied pulse.We showed in [5] that X verifies: In addition, we showed that there exists K > 0 such that: All these models fall in the hypothesis of theorem 10 in [5], so existence and uniqueness of solutions to problem (2) holds.
The main topic is now to know if given a periodic source term g, there exists a unique solution that is also periodic.Also, starting from any initial condition, does the solution converge to the periodic state, and at which speed ?
3 Existence and uniqueness of periodic solutions Notation For any functionnal space H and T > 0, we will designate by C T (H) the subspace of functions of C([0, +∞), H) that are T -periodic :

And let
This section will be dedicated to the proof of the following theorem : Let T > 0 and g be a T -periodic function of C T (H 1 (∂Ω)) ∩ W 1,1 (0, +∞), L 2 (∂Ω) .There exists a unique T -periodic solution U to the following problem: Also, let V 0 ∈ P H 1 (Ω) and V be the solution to where C is a real constant.
Proof: Existence and uniqueness of the solution V is given by [5].In order to prove existence and uniqueness of periodic solutions to problem (9), we introduce the Dirichlet-to-Neumann (or Steklov-Poincaré) operators Λ c and Λ e on the interface for the Laplacian.Denote by ν c (resp.ν e ) the unitary normal to Γ directed from the inside to the outside of O c (resp.O e ).Λ c and Λ e are defined by: Using these operators, the transmission condition (2d) of the transmembrane potential can be rewritten on the manifold Γ.For the sake of readability, we now denote by u the transmembrane potential difference [U ] Γ across the membrane.The condition now reads: where e Λ 0 g is the external source term which has been brought to Inria the interface via another Steklov-Poincaré operator: In [5], we have already proven that the operator (H 1 (Γ), A) is m-accretive.
In order to prove the existence of periodic solutions, it is necessary to verify firstly the boundedness of the solution u to the following problem, equivalent of (2): where A is a m-accretive operator and S ep the conductivities previously defined in eqs.( 5) to (6).We have, for all t > 0, It is clear that for all conductivity models S ep (t, u)u, u is positive.Hence thanks to Young's inequality: We can now invoque Gronwall's lemma to prove boundedness of solutions in finite times: ∀t > 0, We will now focus on the existence and uniqueness of periodic solutions.Let G be a Tperiodic source in C T (H 1 (Γ)), and A a m-accretive operator in L(H 1 (Γ)).We will show that there exists a unique periodic solution u satisfying : Let u 0 and v 0 ∈ H 1 (Γ) be two initial states and u, v the associated solutions to problem (18).Let w = u − v. w satisfies the homogeneous problem: We have, for t 0, and since A is m-accretive, We will show that for each conductivity model, the term • In the case of the β-model, it is easy to verify that for any λ 1 , λ 2 with λ 1 = 0, is positive, since β is an increasing function on [0, +∞).
• For the X-model, the ODE has the following solution Therefore the mapping λ → λX(t, λ(t)) is monotone.Using the same argument as for the β-model, we infer that for all λ 1 , λ 2 , λ 1 = 0 For each conductivity model, we have, for t 0, Let Φ(u 0 ) := u(T ).Writing the previous equation in t = T , we have Inria with e − S 0 Cm T < 1.Therefore, Φ is a contraction on L 2 (Γ) and there exists a unique initial state u 0 satisfying Eq. ( 18).We have also proven the asymptotic convergence of any solution to problem (10) to the periodic solution.
Remark 1.Note that the H 1 regularity is required on u and v so as their values are L ∞ , thanks to the injection H s ֒→ L ∞ for s > d 2 .If we had considered the three dimensional case, then more regularity would have been necessary: u 0 and v 0 should be in H 3/2 (Γ) and G should be in and In this case, the convergence speed depends not only on the base conductivity of the membrane, but also on the shape of the cell and the extracellular conductivity, that are present in the definition of A. We will detail in the next section a way to highlight the influence of these parameters for a circular cell.
4 Estimations of the convergence speed in the case of the circular cell.

Coercivity of the operator A
Under specific conditions on the cell shape, it is possible to refine the previous estimation of the convergence rate to periodic solutions, by studying the operator A.
Let O c be a circular disk of radius R 1 and O e a concentric ring around O c of outer radius R 2 > R 1 .Let H 1/2 p (Γ) be the space of functions of H 1/2 (Γ) with a zero mean value.Then the operator (H

M. Leguèbe
with constant terms u 0 = v 0 = 0, the boundary conditions on Γ and ∂Ω lead to Using the same method for Λ c leads to Then e ikθ e −ilθ dθ, according to the definition of • H 1/2 (Γ) given by [10], p. 141, with Therefore A is coercive.
Remark 3. Remark that we cannot extend this property to functions with a non-zero mean value, as a consequence of the non-invertibility of Λ c .For example, constant functions u = u 0 = 0 verify Au, u = 0 < u Remark 4. Also note that it is possible to consider a cell embedded in a uniform electric field by making R 2 grow to infinity and considering only the case |k| = 1.Under these conditions, that are the most frequently found in experiments on cells, the ratio (R Remark 5.It is possible to extend the coercive property of A to the three dimensional case, using spherical harmonic functions instead of usual Fourier series.Following the same computation, Inria the estimation on the coercivity constant differs from the 2D case by a factor 2 : Going back to problem (19), the coercivity of the operator A leads to and for all t 0 In this case, the convergence speed also depends on parameters included in the constant C A : the extracellular conductivity and the size of the cell.These dependencies will be illustrated in the next section.

Numerical validation
In order to solve numerically equation ( 2), the same scheme as in [5] was used.It is based on a method developped by Cisternino and Weynans [2], using second order finite differences for spatial discretization and adding unknowns at the intersection points of the interface Γ on the cartesian grid for a special treatment of fluxes.The details of the method are fully explicited and its convergence proven in a another previous publication dedicated to the numerical aspects of the problem [6].
Simulations were all performed in a two dimensional domain with the parameters presented in table 1.Since the numerical method is based on a cartesian grid, the simulated domain is a square [x min , x max ] × [y min , y max ] that is large enough compared to the radius of the cell, which is circular.
The chosen boundary conditions simulate a uniform electric field of amplitude E. It consists in a Dirichlet condition in the x direction, and a homogeneous Neumann condition in the y direction.In this case, effects due to the presence of the square-shaped boundary on the transmembrane potential difference (∆TMP) can be neglected.
The problem to solve reads (27)

Results of the linear model
Before studying the influence of the electropermeabilization model on the convergence towards a periodic solution, we will focus on the variation of the coercivity constant C A with the extracellular conductivity σ e .In this section, we will consider a constant surface conductivity of the membrane.
The chosen pulse parameters to study the linear effects are not commonly used.Using a clinical pulse protocol with a linear model does not lead to relevant convergence results since Table 1: Parameters of numerical simulations, that are mostly taken from [5].EP stands for electropermeabilization, EPd for electropermeabilized. pulses are separated by a time that is long enough to let the cell retrieve its initial state.The solution is therefore already periodic for any pulse intensity.To prevent this, the pulses frequency was considerably augmented to avoid the complete discharge of the membrane between two pulses.Moreover, the membrane capacitance C m was also changed to 1 (instead of ∼0.01) to extend the duration of the charge and discharge of the membrane.This allows also direct comparison with the estimation of C A in equation ( 25).The periodic solution u to problem (18) is obtained by performing long-time simulations of 50 periods, and verifying that the last obtained period does not vary from the previous one more than a ratio of 10 −10 : if T f designates the final time of the simulation, Under these conditions, the last period is considered as the periodic solution.We then compute the L 2 error e(t) := u(t) − u(t) L 2 (Γ) .The error is fitted with an exponential function where A and B are constants, as illustrated in figure 2. The constant B can be considered as an estimation of the convergence speed C A .
The process is repeated for different extracellular conductivities σ e and radiuses R 1 .The Fourier coefficients of the solution are computed to check whether the influence of the squareshaped simulation box is important or not.All simulations showed a solution with greatest Inria coefficients u k , |k| = 1 being u |3| , 1000 times lower than the main harmonic u |1| .Figure 3 shows that the dependency on σ e and R 1 of the decay rate B S 0 + C A is satisfied.The evaluated constant B has the same order of magnitude as the estimation of C A in (25).).The first fitting constant has the same order of magnitude as R −1 1 ∼ 167 000.The other constant is of the order of (σ −1 e + σ −1 c ) −1 ∼ 0.24.

Results of the static conductivity models
The same method was used to study the decay rate to periodic solutions when the membrane conductivity is a non-linear function of the transmembrane potential difference.Simulations showed that the convergence of the solution to the periodic state is still an exponential function of time: e(t) ∼ e −Ct e(0).
even in the non-linear case.The constant C still includes the base conductivity of the membrane and the coercivity of the operator A. We propose to note C = S 0 + C A + C Sep , where C Sep accounts for the non-linearity of the conductivity.
Figure 4 shows the behaviour of the quantity C when the extracellular conductivity varies.The dependency on σ e is similar to the linear case, the difference being only a constant factor.This can be explained by the increase in the maximum value of the ∆TMP that is obtained for each simulation.Lower conductivities induce a limited rise of the ∆TMP, contrary to high conductivities, for which the ∆TMP is close to the threshold voltage.Then the Lipschitz constant of λ → λS ep vary accordingly to the maximum value of the ∆TMP λ m : since λ m < V rev in all our simulations.Figure 5 shows that the simulation results are in agreement with this estimation, although the proposed constant for the β model is one order of magnitude higher.

Dynamical model
Figure 6 shows a different evolution of the non-linearity for the dynamical model.For lower conductivities, the behaviour is the same as the static cases, since the applied pulse train is not  sufficient to reach the voltage threshold.However, as the extracellular conductivity rises, a drop of the convergence speed occurs, which depends on the dynamics of X.If τ ep is low enough, the dynamical model has a behaviour which is the same as the β-model, whereas slower dynamics do not affect the convergence speed more than the linear constant.Since the choice of the model has a quantitative impact on the convergence speed to a periodic state, it may be important to translate the results of figure 6 in terms of pulse application.Figure 7 gives the number of pulses that are necessary to obtain a periodic response of the cell.This number is computed so as the solution is no more than 1% different from the periodic solution in L 2 (Γ) norm:

Conclusion
In this paper we study the evolution of solutions to several models of electropermeabilization derived from Kavian et al. [5] in the specific case of periodic sources.We show existence and uniqueness of periodic solutions, as well as convergence to these solutions given any initial condition.Using a simulation tool that was created specifically to solve electropermeabilization problems, we validate theoretical estimations of the convergence speed to periodic solutions.
We emphasize that static models of membrane conductivity give an overestimation of the convergence speed, compared to the dynamical model.Since the permeabilization is linked to the conductivity, those models predict a higher permeabilization level of the cell for short pulse treatments.The number of pulses needed to obtain a given conductivity is then lower than predictions of the dynamical model, and the cell could not be actually porated.

Figure 2 :B (s − 1 )
Figure 2: 2(a) : Transmembrane potential difference (solid) and periodic solution (dotted) of the linear problem, with σ e = 0.2 S.m −1 .2(b) : Fit of the L 2 error e(t) between the ∆TMP and the periodic solution, giving the constant C A .

Figure 3 :
Figure 3: Linear conductivity model : evaluated decay rate B of the L 2 error depending on the extracellular conductivity σ e with R 1 = 6 µm (3(a)) and depending on the radius of the cell R 1 with σ e = 0.5 S/m (3(b)).The first fitting constant has the same order of magnitude as R −11 ∼ 167 000.The other constant is of the order of (σ −1 e + σ −1 c ) −1 ∼ 0.24.

Figure 5 :
Figure 5: 5(a): Maximum ∆TMP during simulations for the β (solid) and λ 2 (dashed) static models.5(b): Comparison of simulated and estimated Lipschitz constants C Sep for the β model.Solid lines : difference between non-linear and linear decay rates.Dashed lines : estimations from eq (30) with λ m from figure 5(a).

Figure 6 :
Figure 6: Decay rate to the periodic solution for the dynamical case S m (λ) = S 0 + S 1 X(t), compared to the evaluations already performed for the linear and the β models.Simulations of the dynamical model were made for several characteristic durations: τ ep = 1 µs(+), τ ep = 10 µs(△).