Comparison and calibration of different electroporation models-Application to rabbit livers experiments

Electroporation is a complex phenomenon that occurs when biological tissues are subjected to electric pulses. The clinical interest for the phenomenon has constantly increased for the last decades. Indeed, electroporation makes it possible to either kill directly the cells in the target region (tumor) or to introduce molecules into living cells. However, one of the main limitation of using electroporation in the clinical routine comes from the technical difficulties raised by such therapies, in particular it is difficult to well determine the treated zone. Numerical modeling of the electric field magnitude could provide a powerful strategy to assess the treatment efficacy: thanks to well-designed models, the computation of the distribution of the electric field is achievable, providing a numerical evaluation of the treatment. The main objective of this work is to go further on the patient-adapted numerical modeling of the electric field magnitude by laying the ground of the possible electroporation models – which will be compared qualitatively – and their calibrations. This will be done in the framework of bioelectrical measurements on rabbit livers that come from the literature.


Introduction
Electroporation (EP) is a complex phenomenon that occurs when biological tissues are subjected to short, high intensity electric pulses. The physical rationale of EP at the cell scale consists of the increase of membrane permeability when the cell is subjected to an electric field strong enough. More precisely, when the cell transmembrane voltage (the difference of the electric potential across the cell membrane) reaches around 1V, then defects are created in the membrane, which thus becomes permeant to extracellular molecules [15,16].
The clinical interest for the phenomenon has constantly increased for the last decades. Indeed, EP makes it possible to either kill directly the cells in the target region (tumor) by a nonthermal mechanism named irreversible electroporation (IRE) [5,6,8,11], or to introduce non permeant molecules (ions, cytotoxic drugs like bleomycin, DNA plasmids, etc.) into living cells, which is referred to as reversible electroporation [9,10,17]. The clinical interests in electroporation-based ablation therapies (EPAs) -either IRE alone or electroporation combined with drug delivery, which is referred to as electrochemotherapy (ECT) -have dramatically increased recently, after publications showing that EPAs provide interesting alternatives to standard non surgical ablative techniques -radiation therapy, radiofrequency ablation (RFA), or cryoablation (CrA)-to treat deep-seated non resectable tumors.
However EPAs are currently mainly limited to cutaneous and subcutaneous tumors. One the main reason of the limitation in the clinical routine use of EPAs comes from the technical difficulties raised by such therapies. Unlike standard ablative techniques which mainly deal with one needle, EPAs need at least two and usually three to four needles (even more for complex tumor shapes): the a priori determination of the treated zone is thus trickier than for standard ablative techniques. Moreover unlike thermal based ablation therapies, as RFA or CrA, for which the monitoring of the temperature in the tumor gives to physician an instantaneous assessment of the treatment, for EPAs the treated region is linked to the amplitude of the electric field, which cannot be monitored instantaneously since the pulse durations are very short (several tens of microseconds per pulse).
Electroporation is tightly linked to the electroquasistatic description of the electric field in the tissue. Therefore numerical modeling of the electric field magnitude provides a powerful strategy to assess the treatment efficacy: thanks to well-designed simulations, the numerical distribution of the electric field is achievable, providing a numerical evaluation of the treatment. However the choice of the model and the patient-dependent calibration, and particularly the patient-dependent organ conductivity, are crucial challenges which are still open.
The aims of the project is to compare the solutions of the different electroporation models found in the literature and to investigate their calibrations for electrostatic descriptions of a tissue from bioelectrical measurements on rabbit livers performed by Sel et al. [19]. After a presentation of 3 different static and 1 dynamical models of electroporation performed in Section 1, Section 2 presents the experimental set-up of Sel et al. and shows how to link the models outputs to the data. After the presentation of numerical schemes used to solve the models, Section 3 proposes a comparison of the solutions in the framework of Sel et al.. Section 4 is devoted to the investigation of parameters estimation from data measurements.

Electroporation models
In this section, we present three static models and one dynamical model, which provide the electrical description of tissue. To compare the different models, the following simple configuration is considered (see Figure 1). The domainΩ is a connected, open smooth subset of R d , d ∈ {2, 3}, and consists of two non intersecting subsets: Ω, which corresponds to the tissue under consideration (rabbit liver for instance in the experiment of Sel et al.) and the needles, represented by two parallel cylinders E ± . The needles are set at the isopotential ±g respectively. The outer boundary is denoted by Γ out := ∂Ω \ E ± . Figure 1. The typical geometrical configuration consists of the tissue domain Ω in gray, deprived of the two needles E ± , in white. The outer boundary Γ out is represented in bold.

The standard static model
The most used model to describe tissue electroporation consists of the phenomenological electrostatic problem. The tissue is described as a conductive medium, whose conductivity tensor σ depends on the amplitude of the electric field −∇u. The model reads then The tissue conductivity consists of a 4 parameters sigmoid function. Typically where σ 0 is the conductivity of the non electroporated tissue, σ 1 is the tissue conductivity of the fully porated tissue, E th is the threshold amplitude for electroporation, and k ep is the slope of the nonlinearity. Here, erf is the Gauss error function. The qualitative behaviour of σ is depicted in Figure 2. It should be noted that this choice of σ is largely phenomenological, and as shown in [12], the available experimental data does not seem sufficient to characterize the dependence of the conductivity on the electric field.
1.2. The electric circuit approach of Voyer et al.
In [20], Voyer et al. proposed a biphasic dynamical model based on the description of an individual cell and surrounding matrix as an electric circuit. The ODEs at the cell level are formally generalized to PDEs at the tissue level. It describes the electric potential outside cells φ e and the electric field inside cells J c . The parameters are the extracellular and intracellular electric conductivities, respectively σ e and σ c . The conductivity of the cell membrane σ m depends on time in a way which mimics the effects of poration, i.e. the appearance of holes on the membrane, and permeabilisation, that is the degradation of membrane molecules. Both phenomena increase the conductivity. The resulting system reads where The functions X 1 and X 2 are the respective degrees of poration and degrees of permeabilisation. They satisfẏ where β 1 and β 2 are 2-parameters sigmoid functions of the form where k i is the stiffness of the sigmoid and x trans i the transition threshold. A simplified static version of this model allows to compare it with the standard model. Indeed, the static version of this model writes ∇ · (σ e ∇φ e + J c ) = 0, where σ m is a sigmoid function similar to (2). Written in φ e and E m , the equations become Problem (5) makes it possible to define the equivalent tissue conductivity. Actually, thanks to (5b), the potential φ e satisfies The equivalent tissue conductivity can thus be defined as A simplified version of the problem consists in assuming that σ m depends on ∇φ e instead of E m . Then the simplified problem reads which is similar to the standard model (1) since the function σ eq defined by is a sigmoid function.

The static bidomain model
In electrocardiology, the so-called bidomain model has been proven to be the homogenisation limit of the cell scale electric potential. It consists of 2 electric potentials, u e and v, which satisfy the following PDE system proposed in [7]: where S m = A mSm , with A m (m −1 ) the ratio of membrane area by unit volume andS m (S.m −2 ), is the conductance of the cell membranes. This model enables to link rigorously the membrane conductance S m to the electric potential thanks to an homogenization procedure presented in [3,4]. In this sense, it is more physiological than the previous phenomenological models (see Remark 1). To account for the electroporation phenomenon -that is the increase of the membrane conductance -we assume thatS m follows a sigmoid function of the norm of the outer electric field −∇u ẽ Remark 2 (On the nonlinearity). In cardiac electrophysiology, the nonlinearity usually depends on the homogenised transmembrane voltage u e − u c . However such a nonlinearity is inconsistent with the electroporation phenomenon. Actually, assuming that g − = −g + is constant, the symmetry along the axis (Oy) implies that u e , u c and the transmembrane voltage vanishes along (Oy), which is inconsistent with the fact that electroporation occurs between the needles. For this reason, a nonlinearity depending on ∇u e is chosen, which is relevant with the nonlinearity of the static nonlinear model.

Data measurements
The experimental set-up of Sel et al. consists of an ex-vivo cubic piece of rabbit liver in which two needles are inserted (see Figure 3). Square pulses of different amplitudes are applied on the needles and the corresponding intensities are recorded. The data set consists of the measurements of the electric intensity that flows through one needle, say E + for instance. tion in tissue [17]- [19]. Electrode sue permeabilization are parallel late electrodes as well as needle ys [19]- [21]. have shown that the extent of cell depends on electric field intensity, n with surrounding cells [14], [15]. eabilized cells such as cells volume dium, membrane conductivity, cell smembrane potential (TMP) affect . The relationship between these ductivity of cells in suspension has 16]. Besides effective conductivity abilization on a tissue level is the As the field results from a voltage odes, the electrode configuration tion in tissue [17]- [19]. Electrode sue permeabilization are parallel late electrodes as well as needle ys [19]- [21]. have shown that the extent of cell depends on electric field intensity, n with surrounding cells [14], [15]. eabilized cells such as cells volume dium, membrane conductivity, cell smembrane potential (TMP) affect . The relationship between these ductivity of cells in suspension has 16]. Besides effective conductivity abilization on a tissue level is the As the field results from a voltage odes, the electrode configuration tion in tissue [17]- [19]. Electrode sue permeabilization are parallel late electrodes as well as needle ys [19]- [21]. ssue permeabilization it is required ssue to E intensities between the se in advance a suitable electrode rameters for the effective tissue , electric field distribution in tissue e treatment, which can be achieved id tests [18], [24] with models of owever, modeling of electric field anding due to heterogeneous tissue plex geometry. Analytical models simple geometries. Usually they nsional problems and tissue with rties [25]. Therefore, in most cases ques are still more acceptable as ling three-dimensional geometries ies. For that purpose mostly finite nite difference method are applied. ve been successfully applied and computed and measured electric In this paper, we present the first model which describes tissue permeabilization by taking into account tissue conductivity change. The model consists of a sequence of static models (steps), which describe E distribution in discrete time intervals during permeabilization. In this way model presents dynamics of electropermeabilization since in each step the tissue conductivity is changed according to distribution of electric field intensities from the previous step. For that purpose tissue conductivity in the model is expressed as a function of electric field intensity. Sigmoid dependency between specific conductivity and electric field intensity is used. Estimation of the sigmoid function parameters is based on current measurements. Namely current indicates the extent of permeabilization [26], [27] through the change in the tissue conductivity [28], [29], [46]. The proposed model of tissue electropermeabilization is then validated on experimentally obtained total current measurements and areas of reversibly and irreversibly permeabilized rabbit liver tissue.

A. Experiments
In vivo experiments were performed at the Institute Gustave-Roussy, France on rabbit liver tissue in accordance with European Commission Directives and French legislation concerning animal welfare. Three rabbits were used in the experiments. Animals were kept anaesthetised for the whole duration of the experiments. A subxypoid incision was made and the liver was gently exteriorised and exposed to electrical treatment. For that purpose two parallel needle electrodes as shown in Fig. 1 were inserted perpendicularly to tissue surface approximately 7 mm in depth. In experiments three different needle diameters were used: , and . The inner distance between the needles was always 8 mm as in [22]. Eight rectangular monophasic pulses of 100 duration and 1 Hz repetition frequency were applied. Pulses were delivered by pulse generator x 2 E ± I V: electric potential (rV electric field) I : conductivity (  Taking advantage of the symmetries of Figure 4, the computation domain is restricted to a quadrant with homogeneous Neumann conditions on the right, top and bottom borders, homogeneous Dirichlet condition on the left boundary and the non homogeneous Dirichlet condition on the needle.
The expression of the numerical electric intensity depends on the choice of the models. Denoting by E + the length of the electrode in the perpendicular plane to the 2D simulation plane, one has where V is the solution to (1), is the solution to (5), or to (3) where φ e is the solution to (7), is the solution to (8).
The above definition of the electrical intensities of the models involve the gradient of the solution along the needle. This is numerically unstable since the electric field is the most intense nearby the electrodes, and thus numerical instabilities may appear. To smoothen this behavior, it is possible to introduce a specific function w such that the surface integral is replaced by a volume integral. More precisely, let w be defined by Then, thanks to appropriate integration by parts one has where φ e is the solution to (7), or to (3) (10c) 3. Numerical discretization and qualitative comparison between the models

Discretization and numerical approximation
In this section, we present the discretizations that have been used to solve the different models. All these strategies have been implemented using the Python finite element library FEniCS [1].

Static models: illustration on the standard static model
We want to solve (1) which we recall here: which is a non-linear problem since σ depends on ∇u. In Breton et al. [2] a modified fixed point has been presented and the numerical convergence of this scheme has been shown. In this work, we compare two other strategies. The first one consists in using iteration schemes coming from the time-discretization of the following non-linear evolution equation ∂ s u − ∇ · (σ( ∇u )∇u) = 0 in Ω , with the same boundary conditions as above, and where s acts as the pseudo time variable. We follow the same strategy for the static version of the biphasic model and for the bidomain model. Depending on the eigenvalue of the linearised problem, one can expect that u behaves asymptotically as the static solution. We explicit the strategy (which can be easily extended to the other static models) for the standard model. First, we consider the decomposition u = v + u lin , where u lin solves the linear static problem associated to (1), that is the Laplace equation on Ω with boundary conditions (1b). Thus, we obtain the following equation on v with homogeneous Dirichlet boundary conditions, since u lin carries the boundary conditions of (1). By discretizing the s derivative by finite differences with step ∆s, multiplying this equation by a test function ϕ ∈ C ∞ c (Ω) and integrating over Ω, we get The term on Γ out vanishes because of the homogeneous Neumann boundary conditions, and we can restrict to test functions with zero trace on E ± by integrating the Dirichlet data in the function space, see below. Hence, the variational problem to solve reads A(v n+1 , ϕ) = L(ϕ) where This iterative method is initialized by taking v 0 = u lin . We then proceed by approximating v, u lin using P 2 finite elements, i.e.
where P 2 is the space of polynomials in two variables, of degree at most 2 and T h is a triangulation of Ω. We also introduce This problem is well-posed according to standard finite element analysis arguments.
The second approach consists in using Newton's method directly on the finite element system. Consider the problem consisting in finding u h ∈ V h g ± such that where A is given by which is the discretized weak formulation of (1). Assuming that some u h 0 satisfying the boundary conditions is known, one can try to solve this system of non linear equations using the Newton's method. Such a u h 0 is given by solving Laplace's equation on Ω with boundary conditions (1b). In FEniCS, the computation of the corresponding Jacobian can be performed automatically. For the models we consider, this approach is faster and more stable compared to the evolution-based iterative method presented above, although it requires a good guess for initialization. In practice, such a good guess can be found by starting from the solution to the linear problem and increasing gradually the value of k ep until the desired value is reached. The Newton's approach is also simpler, in the sense that is does not require a choice of a step size (∆s in the evolution approach).
A qualitative comparison of the two methods can be seen in Figure 5. σ( ∇u ) is projected on piecewise constant discontinuous elements. On the right of the domain, the profile given by the Newton's method is noticeably smoother.

The dynamical biphase model
Concerning the dynamical biphase model, we use an iterative time discretization method. We denote the time step by ∆t, so that t n = n∆t. At the n th step, we compute the membrane electric field, and its amplitude where d c is the mean of the cell diameters. From there, we solve the differential equations for the degrees of poration and permeabilization X 1 and X 2 using an explicit Euler scheme Then, we compute the membrane conductivity Similarly to the static models, we multiply (3) by the test functions ϕ 1 ∈ C ∞ (Ω) and ϕ 2 ∈ C ∞ 0 (Ω) 2 respectively and integrate by parts to obtain We discretize again on P 2 finite elements and integrate the Dirichlet boundary conditions in the finite element space. We then have to find φ n+1,h The functionals A and L are given by:

Qualitative comparison between the models
In this section we show the characteristic solution of all three models presented in Section 1. We consider an electrode diameter of 0.7mm, a voltage of 800V. Results for the standard, the bidomain and the dynamical biphase models are illustrated in Figures 6, 7 and 8 respectively. In all cases, the contour lines of u and u e look similar.

Estimation strategy: illustrations on the standard static model
As we have seen in the previous section, the solutions of the different models are realistic and can be qualitatively compared. They allow to represent the electroporation phenomena. The aim of this section is to investigate numerical calibration strategies of these models. We will focus for this work on the standard static model.
First, the estimation strategy is introduced. Then, it is illustrated with synthetic data and finally, it is applied on real data [19] presented in Section 2.   Table  1], with a 800V voltage step. It is easy to see that σ m is mainly driven by the value of X 1 , i.e. the poration process. At the lower left is the evolution of the current intensity flowing through the electrodes, the spike qualitatively matches experiments, see [20].

Parameters estimation algorithm
We denote by θ ∈ Θ the vector which concatenates all the parameters which have to be estimated. We denote by y ∈ Y the solution and by z ∈ Z the observations. We denote by · Θ,(P ) −1 and by · Z,R the norms defined by · 2 Θ,(P ) −1 = · , (P ) −1 · and · 2 Z,R = · , R · .
The static problem can be rewritten as A(y, θ) = 0, where A corresponds to the model operator. We denote by C the observation operator which relies y to z z = C(y).
and we denote by O the following operator O : θ → y, such that A(y, θ) = 0.
Our objective is to inverse the operator Ψ = C • O. In order to do that, the criterium that we want to minimize is given by where θ corresponds to an a priori value of θ. This criterium corresponds to a likelihood functional where all disturbances are assumed to be Gaussian θ ∼ N (θ , P 1/2 ) and z ∼ N (Ψ(θ), R −1/2 ).
If we want to use N z measurements denoted by (z k ) 1≤k≤Nz associated to (R k ) 1≤k≤Nz , we have to consider the following criterium which can be rewritten as It exists many methods to solve this kind of problem as for example as a non-exhaustive list: gradient descent based methods ; stochastic strategies as for example Monte-Carlo based strategies or expectation-maximization algorithm. In this work, we decided to use a stochastic algorithm based on an unscented transform [14]. Let (ω i ) 1≤i≤p and (e (i) ) 1≤i≤p be some weight parameters and directions in Θ such that This is called a set of sigma points and many different sets of sigma points have been proposed in the literature.
In this work, we will consider the simplex sigma points, which allow to consider only dim Θ +1 points [13]. The algorithm reads as follows: • Initialisation:θ 0 = θ and P 0 = P .
• Generation of p particles of meanθ k and covariance P k : • Computation of the observations generated by the solution of the system for each particle: • Computation of the mean and of the covariance of Z i : • Computation of the covariance between θ (i) and J(θ (i) ) (sensitivity of the observations to the parameters): • Computation of the new value of the parameters and of their covarianceŝ

Estimations on synthetic data
As said previously, we focus on the standard static model in this section. The parameters are then the parameters of the sigmoid function: E th , σ 0 , σ 1 and k ep . We want to validate on synthetic data the estimation strategy. To do that, we consider a framework close from the real data. Following [19], we build synthetic data with intensity measurements for 5 different voltages and 3 electrode sizes (the same of the article). The corresponding observation space has then dimension N z = 15. The synthetic data are obtained using E th = 5.75 × 10 4 V /m, σ 0 = 0.065 S/m, σ 1 = 0.1483 S/m and k ep = 10 −3 . The parameters that we want to estimate are all positive which will be constrained by a log-transformation during the estimation procedure. We assume that Concerning the uncertainties quantification of the observations, we assume that R −1/2 = 0.025 Id Nz,Nz which corresponds to a noise standard deviation of approximately 20%: R −1/2 ≈ 20% 1

Nz
Nz k=1 |z k |Id Nz,Nz ≈ 0.025. First of all, we start by an individual estimation of the parameters and the synthetic data are not noisy. Table 1 presents the results (with 2 digits of precision). The algorithm stops when the convergence is reached.
As a conclusion, it is possible to estimate three parameters: E th , σ 0 and σ 1 . The sensitivity of σ 0 is weaker and more iterations are needed (compared to E th and σ 1 ). Then, k ep is not identifiable: there is no sensitivity associated to this parameter. As k ep is not identifiable, we fix it at its true value and we make a coupled estimation of the three other parameters, see the third column of Table 2. We consider also a case where only one size of electrodes is considered. The results are very encouraging and show that it is possible to estimate the three parameters together. The results are slightly better when three electrodes are considered. Finally, we add a gaussian noise of standard deviation σ n = 0.005, 0.01 (resp. 5% and 10%) on the synthetic data to see the sensitivity of the estimation to measurement noise. In the case of σ n = 0.01, we also try to estimate using  Table 3. Coupled estimation of E th , σ 0 and σ 1 (2 digits of precision). The true and estimated values are given after the log-transformation. The third (resp. fourth) column corresponds to σ n = 0.005 (resp. σ n = 0.01). The last column is for σ n = 0.01 and only one electrode is considered.
only one size of electrodes. The results are given in Table 3. As expected, the more the noise on the data is important, the more it is difficult to well estimate the parameters but the results are still encouraging. More precisely, with noisy data, it is crucial to consider measurements performed with different size of electrodes.

Estimations on the data set of Sel et al.
We now turn to the problem of parameter estimation for real data.

Dataset
The dataset used is the one provided by [19], the corresponding experimental setup is recalled in Figure 3. It provides intensity measurements for 5 different voltages and 3 electrode sizes, which are collected in Table 4. As for the synthetic data, the corresponding observation space has then dimension N z = 15.

Estimation and validation
As was made clear in Section 4.2, the sensivity on k ep in the case of a 4-parameter sigmoid is very low. Trying to estimate this parameter would most likely result in an unreliable value. To avoid having to choose an arbitrary value for k ep , we consider a 3-parameter sigmoid instead, of the form The remaining three parameters, E th , σ 0 and σ 1 have to be estimated jointly, since we have no a priori value for them. We assume that the parameters are centered around the values given in Section 4.2, i.e.: log E th ∼ N (log(5.75 × 10 4 ), 2), log σ 0 ∼ N (log(0.065), 0.2), log σ 1 ∼ N (log(0.1483), 0.5) .
The noise on the measurements is modeled by R −1/2 = 0.2, which corresponds to a noise of 25%. When all electrodiameters are considered together, the convergence of the estimation procedure is shown in Figure 9. As is the synthetic data case, the convergence of σ 0 is much slower than the other two parameters, which requires a large number of iterations.
For each individual electrode size, the convergence of the estimation procedure is shown in Figure 10. The estimated values differ significantly from one electrode diameter to the next, and σ 0 converges again slower than E th and σ 1 .
As we can see, the estimated parameters differ significantly between electrode diameters, especially for 0.3mm. We solve the direct problem with parameters estimated for each diameter individually ; and for the parameters estimated when considering all diameters simultaneously, see Figure 11. We can see the resulting distribution of σ( ∇u ) for the three electrode diameters and for voltages of 200V and 1000V . The distribution is very different, especially for 1000V . Table 6 presents the maximum relative error on the intensity, and shows that the best results are obtained when estimating the parameters only on the electrode of diameter 1.1mm. This is not  Figure 10. Estimation of E th , σ 0 and σ 1 , after log-transformation. All voltages are considered, but electrode diameters are considered separately. The thin lines denote the standard deviation of the uncertainty on the estimation.
surprising regarding to the relatively low values for the current for a diameter of 0.3mm, the quasi-equivalent values for diameters of 0.7 and 1.1mm and also the relatively low values for low voltages (see Table 4).
Maybe, better results could be achieved by understanding better the impact on the measurement noise using for example synthetic data. One strategy could be for example to consider different values of the standard deviations of the observations in order to give less importance to the most noisy data.    Table 6. Maximum relative error on the intensity across voltages, for each set of estimated parameters (one for each electrode diameter plus one for the estimated on all three simultaneously).

Conclusion and perspectives
In this paper, we have compared different models of electroporation: the standard static model, the dynamical biphase model and the bidomain model in the same framework built using bioelectrical measurements on rabbit livers proposed by Sel et al. [19]. The solutions of the different models are realistic and can be qualitatively compared. Then we investigate a numerical calibration of the parameters only on the standard model. The strategy is based on a stochastic algorithm using an unscented transform for estimating parameters of static systems. The results on synthetic data are convincing. Finally, the estimation strategy is applied on the rabbit measurements and the attempt is promising even if the best fit is not achieved using all geometries simultaneously. As explained previously, this could be due to the measurements noise. This work allows to lay the ground of the calibration of electroporation models but many improvements have to be studied. The work on synthetic data has to be continued in order to better understand the estimation difficulties due to the measurements noise. Furthermore, the parameters estimation has to be applied on the bidomain model. Finally, concerning the dynamical biphase model, a dynamical version of the estimation strategy has to be written. The main difficulty lies in the fact that we need an efficient strategy which does not involve prohibitive computational times. To address this issue, Reduced-order Kalman filter method, introduced by Moireau and Chapelle [18], will be investigated in a near future.