NUMERICAL MODELLING CHALLENGES FOR CLINICAL ELECTROPORATION ABLATION TECHNIQUE OF LIVER TUMORS

. Electroporation ablation is a promising non surgical and minimally invasive technique of tumor ablation, for which no monitoring is currently available. In this paper, we present the recent advances and challenges on the numerical modeling of clinical electroporation ablation of liver tumors. In particular, we show that a nonlinear static electrical model of tissue combined with clinical imaging can give crucial information of the electric ﬁeld distribution in the clinical conﬁguration. We conclude the paper by presenting some important questions that have to be addressed for an eﬀective impact of computational modeling in clinical practice of electroporation ablation


Introduction
Electroporation consists in applying high voltage short pulses to cells (typically several hundreds of volts per centimeter during one hundred microseconds) in order to create defects in the cell membrane.Since the last decades, the phenomenon has been widely studied both in vivo and in vitro experiments.Interestingly it has been shown that this method selectively affects the cell membrane, which becomes permeable if the high voltage is high enough.It thus makes it possible to either introduce non permeant molecules (ions, cytotoxic drugs like bleomycin, DNA plasmids) into living cells, which is referred to as reversible electroporation [3,46] or to kill directly the cells in the targeted region (tumor) by a nonthermal mechanism (this is named irreversible electroporation [11,12,15].

Electroporation-based therapies: an alternative to non surgical ablation therapies
From the clinical point of view, electroporation-based therapies (EBTs) provide interesting alternatives to standard ablative techniques.Actually radiotherapy, microwave and radiofrequency ablation, or cryotherapy for instance, are widely used in treatment of patients with deep-seated solid tumors, when surgery cannot be used.Such methods rely upon the indiscriminate use of thermal energy to induce necrosis of the affected cells (cancerous or not).In particular, due to the thermal diffusion in tissues, these techniques may result in dramatic damages to the nearby structures including blood vessels, vital organs close to the target, bile ducts, and nerves.Compared to these standard ablative therapies, EBTs -mainly irreversible electroporation (IRE) and electrochemotherapy (ECT) -have several advantages: -As electroporation is not based on thermal mechanisms -since pulses duration is too short for Joule heating-the tissue scaffold is preserved: this makes such therapies suitable for tumors in proximity of vital organs [4,38,39], and it improves the healing.-It has been proved that electroporation preserves the integrity of the main vessels, while it disrupts neovascularization [46], leading to tumor hypoxia.-Due to a vascular lock effect [23], electroporation greatly reduces bleeding.-As proteins are not denatured, antigens may stimulate an immune response against remaining cancer cells, including metastasis located far from the target tissue.Note that this question is still unclear with controversial publications: in [1], no evidence of immune response has been found (unlike for electrochemotherapy as reported in [6] while Thomson et al. [49] claim that some evidence of an immune response has been shown in animals. In spite of these advantages, that have been well-established in the experiments, the clinical use of EBT is still mainly limited to cutaneous and subcutaneous tumors.It is worth noting that the hospital CHU J. Verdier has become a leader in the use of IRE for liver tumors, in particular hepatocellular carcinoma (HCC) to the liver: we thus benefit from this knowledge to focus our proposal on HCC.A first reason of the limitation in the use of EBTs comes from the technical difficulties raised by such therapies: unlike standard ablative techniques which mainly deal with one needle, EBTs 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.In addition, these therapies suffer from the lack of appropriate post treatment protocols to evaluate the tumor response.For standard ablative techniques, the MRI control is usually performed several weeks after the treatment, while recent publications suggest [8,38] that for IRE, early post procedure MRI at 1 day could improve the understanding and the evaluation of the treatment.Another reason for the limitation of IRE in clinic resides in the remaining lack of knowledge of the effect of IRE on cells.More precisely, there is currently no standard criterium to evaluate the procedure on early MRI.Indeed, the image exhibits a white enhancement in the vicinity of the needle location on the T2 weighted MRI modality, which raises several crucial questions (see Fig 1).For instance, is it a reversibly permeabilised or irreversibly permeabilised region?How to assess the success of the procedure on this early image?It is possible to predict this white enhancement with to the needle location?

Objectives and outline of the paper
The aim of this paper is to present our numerical strategies that can provide highly valuable answers to the above clinical questions.In the next section, we present briefly the clinical workflow which is routinely performed at the University Hospital J. Verdier at Bondy.Then we present the most used numerical model of tissue electroporation we investigate numerical the sensitivity to needle location and to conductivity value.Section 5 proposes a strategy which has been presented in [20].In conclusion, we draw some clinically relevant perspectives to increase the impact of the numerical modeling on the clinical practice of EBTs.

Clinical electroporation of liver tumors
The clinical workflow as routinely performed at the University Hospital J. Verdier at Bondy has been detailed and validated by Sutter et al. [47,48].It is somehow the minimal clinical procedure, which can be performed in most of interventional radiology operative room equipped with new generation of angiographic suit including 3D cone beam CT acquisition capacity [48].The clinical workflow is split into 3 steps as described in [20].The treatment region exhibits a white enhancement for which radiology interpretation is still controversial [45].

Clin. Step 1
The diagnostic stage provides a preoperative CT scan or MRI where the hepatic capsule, the tumor and the main liver structures are determined several days before the treatment (Fig. 1 (Left)).

Clin. Step 2
The day of the procedure, IRE ablation is performed under general anesthesia, as described by Martin et al. [36].As previously described, the needles are percutaneously inserted around the tumor by the interventional radiologist with free-hand technic under combination of real-time ultrasound (US) and 3D Virtual Target Fluoroscopic Display such that the electric field covers the target region.The needle positioning is performed as parallel as possible and the needles location is recoverd by cone beam CT (CBCT) imaging (see Fig. 2).Then the pulse parameters are set in the pulse generator (number and duration of pulses, voltage per needle pairs, etc.).

Clin.
Step 3 Three days after the procedure, an MRI is performed to assess the treatment efficacy (Fig. 1 (Right)).
The clinical question can be summarized as follows.Is it possible, thanks to the pretreatment imaging and the data acquired during the procedure to interpret the posttreatment imaging?In the following, we investigate a model of tissue electroporation that can be used to address the clinical question.

Standard tissue electroporation model
At the tissue scale, the electroporation modeling is in progress, since the phenomenon is still unclear.The most current models are based on nonlinear static descriptions of the electric field.

The standard tissue electroporation model: a nonlinear static model
The mostly used model to describe tissue electroporation consists of the classic electrostatic problem [35,37,44].The tissue is described as a conductive medium, whose conductivity σ depends on the amplitude of the where E ± denotes respectively the activated anode and cathode needles used for the pulse delivery, Γ out is the boundary of the domain Ω deprived of the activated needles, and g ± are the Dirichlet data at the needle E ± .The classical description of the tissue conductivity consists of a of 4-parameter sigmoid function.For instance 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.
Remark 3.1.The choice of the sigmoid function does not seem essential according to Ivorra et al. [28], but it is probably because of the lack of experimental data.In the code we took a hyperbolic tangent but the function erf can be preferred.
Standard arguments lead to the well-posedness of the nonlinear static model.
Theorem 3.2.Assume that Ω is a Lipschitz domain.We assume that the boundary of Ω has three connected components: the outer boundary ∂Ω, and 2 inner disjoint boundaries E ± .Let g ± ∈ H 1/2 (E ± ) be the applied potentials on each inner boundary.Then there exists a unique potential V in H 1 (Ω) solution to problem (3.1)-(3.2).
Sketch of the proof.Thanks to the monotonicity of σ one can easily show that the above nonlinear PDE is well-posed in H 1 (Ω) if g ± ∈ H 1/2 (E ± ).Actually, denote by U a lift of the Dirichlet condition g ± .For instance let U be the solution to the linear elliptic problem: U obviously belongs to H 1 (Ω).Let V(Ω) be defined by Define by J the functional by − The hypotheses on g ± clearly imply that for any w ∈ V(Ω), |J(w)| < +∞.− Moreover J is differentiable in V(Ω), since for any (w, h) ∈ (V(Ω)) 2 , one has − The strict convexity of J comes from the strict mononicity of σ: where in the above inequality the equality holds only for w = w.− Finally, since one has lim Then the solution V to (3.1) is given by V = W + U .

Numerical method dedicated to clinical applications
The evaluation of the treatment efficacy is crucial information, which is still not monitored in clinical use of IRE.Usually, the treatment planning numerical strategies aim to provide to physicians the optimal positioning of needles to perform the electroporation ablation.Our approach changes this paradigm by letting the physician place percutaneously the needles as best as possible, and then by using this effective clinical positioning of the needles as a starting point of the simulation.We benefitted from a tight collaboration with Olivier Seror, interventional radiologist at Jean Verdier Hospital, to develop the finite difference method-based software IRENA for the IRE Numerical Assessment in clinical routine [22].

IRENA: a finite difference method-based dedicated to clinical IRE
IRENA provides a numerical tool to assess the distribution of the electric field and check whether the tumor is included in the estimated treatment area or not.For this purpose, numerical simulations are performed from the real clinical data: the medical images which give the position of the organ and the tumor, the real position of the needles, and the current graphs of test pulses for tissues conductivity calibration.These data are the input parameters (see the numerical workflow detailed in Gallinato et al. [20]).
In a first step of the procedure, the software makes it possible to setup the generator (Nanoknife R ) by computing the input voltages or nominal electric field to be applied.This essential step avoids over currents and The circle corresponds to the needle/tissue interface, the grid is the numerical grid.makes easier the procedure.In a second step, the electric field is computed, giving the immediate assessment of the procedure.Finally, a retrospective validation can be done afterwards, by quantitatively comparing the numerical results with the observations on MRI a few days after the procedure.Difference in volumes and shape indicators (Hausdorff distance) are provided.As the description of electrical current propagation in biological tissues remains an open issue, IRENA also provides a numerical framework for confronting the existing different mathematical models, as discussed in Section 6.4.For now, the standard nonlinear model is implemented in IRENA but other models will be added in a near future.The discretization is performed thanks to the finite difference method on Cartesian grid.This is the natural choice for computing from the voxel structure of medical imaging.This also gives a significant advantage for future code parallelization.

Discretization close to the needle location
Far from the interfaces, the Laplace operator is discretized with the standard 5-points stencil.As the needles are too thin to be discretized by the grid, the computation of the operator at the point close to a needle relies on a modified ghost fluid method.As shown in Figure 3b, computing the operator at one point u i,j of the first layer around a needle (Fig. 3a) may involve some neighboring points inside the needle (u i,j+1 ), some which are themselves in the first layer u i−1,j and u i+1,j and farther points (u i,j−1 ) which do not require any specific treatment.Neighboring points inside the needle are considered as standard ghost point and are handled as in [17] with the exact boundary condition, whereas neighboring points in the first layer are handled as specific ghost points (Fig. 3c), even if they are in the tissue.Let consider the projection u N on the needle of the specific ghost point.Then, in the Neumann case, the derivative in the normal direction at the specific ghost point is a first order accurate approximate of the known normal derivative at the projected point.In the Dirichlet case, a first order expansion about the specific ghost point towards its projection gives a second order accurate approximate of the known value at the projected point.Hence, the Neumann and Dirichlet boundary conditions on the needle surface are accounted for thanks to Neumann and Robin conditions at the specific ghost points, so that the overall accuracy of the solution is preserved in [21].
The above ghost fluid methods (GFM) require to know the distances to the immersed interfaces (liver capsule, domain bounds, needle surfaces), These interfaces are implicitly described by the level 0 of the signed distance function, obtained from a mask and the Fast Marching algorithm.As the geometrical configuration is static, the redistanciation step is required only in case of change in needle position (pullback, pushforward or new needle insertion) during the procedure.As the electric field derives from the electric potential, the main PDE variable,

Intensity computation
An efficient and accurate computation of the electrical intensity which flows through one needle -say E + -is required to calibrate the tissues conductivity from initial test pulses and their recorded intensities (IRE current graphs) [20].It is defined as However the direct computation of the above integral leads to numerical errors due to the fact that the intensity is along the needle which is placed in the region where the nonlinearity has the most influence on the conductivity.Therefore a small error on the electric field leads to a large error on the conductivity and thus on I E + .To prevent such drawbacks, we propose to compute I E + along an potential closed isosurface C u0 , where u 0 is the value of the chosen potential.This potential is such that C u0 is around and far enough from the needle E + , the Green's formula shows that it leads to the same intensity: The surface C u0 is approximate by the face of the mesh cells, and the fluxes are computed along these faces (see Fig. 4).The harmonic mean of the conductivity and the normal electric field are directly computed at the center of the faces with the known conductivity and the potential values are computed at the center of the adjacent meshes.

Numerical investigation of model sensitivity
It is crucial to investigate the model sensitivity to parameters.Two main types of uncertainties are crucial to analyze.In appendix, we also illustrate the reason why 3D simulations are crucial to provide relevant results.On the one hand, uncertainties can come from the conductivities obtained by the calibration.The second source of uncertainties is due to the errors of the geometric data.In the numerical sensitivity tests, the Hausdorff distance as defined in [43] between the isolines of the electric field magnitude is computed, and the volumes are compared.We first study the influence of liver and tumor conductivities on the electric field distribution.The tissue conductivity is patient-dependent, with more or less variability depending on the tissue and the pathology.Thus, the conductivity value ranges found in literature are wide as stated above.For both tests on liver and tumor conductivity, σ X 0 and σ X 1 stand for the minimal (tissue at rest) and the maximum (fully porated tissue) conductivity values, respectively.As a second step, we investigate the influence of a small translation (3 mm) or inclination (5 • ) of a needle on the electric field distribution.

Model sensitivity to liver and tumor conductivities
Regarding the sensitivity of the electric field distribution to the liver conductivity, equations (3.1)-(3.2) suggest that only the ratio between both linear and nonlinear parts of expression (3.2) matters.Therefore σ L 0 is set to the median value found in literature (σ L 0 = 1 mS cm −1 [19,24,25,31,41]), and σ L 1 ranges from 2.25 to 4.5 mS cm −1 [10,26,31].Figure 5a shows the variation of the electric field distribution for a given position of the needles with respect to the liver conductivity (without any tumor).As mentioned above the strict parallel placement of the needles is not always achievable because of anatomical constraints faced by the operator using percutaneous approaches.In this section, the choice of the needle positions is the typical positioning performed by the physicians in the case of liver tumors located near vital regions.The Hausdorff distance between the two 3D areas (delineated by the dotted and filled lines in the planar views of Fig. 5a) is 2.1 mm, and the difference in volumes is 994 mm 3 .
Both results, in transverse and longitudinal views, highlight the relatively weak influence of changes in conductivity (in the range of values found in literature for liver and HCC) on the electric field distribution.This is corroborated by the corresponding Hausdorff distances.However, the differences in volumes highlight two points: the influence slightly impacts the whole isosurface resulting in a significant difference in volumes, and changes in liver conductivity have a stronger impact than changes in tumor conductivity, as previously highlighted in [30].
Note that an underestimation of the tumor conductivity could result in an overestimation of the treated region (see area below the IRE threshold in Fig. 5b), which may lead to a treatment failure.However, only one pulse per pair of needles has been simulated.The accumulation of pulses in clinical procedure leads to a better coverage of the tumor.In addition, the radiologist may also perform a pullback to enlarge the treated region.

Model sensitivity to the needle locations
As known, the electric field is very sensitive to the needle positions.Figure 6 shows that a small translation (3 mm) or inclination (5 • ) of a needle has a significant impact on the electric field distribution.A translation of 3 mm leads to a Hausdorff distance of 3.4 mm between the 650 V.cm −1 -isolines, and the difference in volumes is 209 mm 3 .An inclination of 5 • of one needle leads to a Hausdorff distance of 4.4 mm, and the difference in volumes is 247 mm 3 .For each test, the Hausdorff distance is larger than for the tests with extreme values of the conductivity.Figure 6a shows that a small error in the needle location may lead to large errors in the prediction of the region affected by the electric field.This aspect shows that the discrepancy between the treatment planning and the clinical procedure may lead to treatment failure.In liver, the actual layout can then vary substantially (sometimes with much larger variations than an inclination of 5 • or a translation of 3 mm).The second step of the numerical workflow proposes a way to overcome this drawbacks by reconstructing the geometrical framework of the actual clinical procedure.
Regarding the tumor location, it has been shown previously that changes in conductivities have only a small influence.However, due to its high conductivity, the tumor behaves as a potential well, which attracts the electric field lines.As a result, the change in the malignant tissue position leads to a local change in the electric field distribution, as shown in Figure 6b.More generally, any structure substantially more conductive than the liver tissue (surgical clips, scar tissue) and located in the vicinity of the target lesion has to be accounted for by the numerical simulations.
Among the structures that may influence the electric field, the liver capsule (Glisson's capsule) behaves like a barrier to the current.It is composed of a fine dense connective tissue layer, which makes it a thin insulating membrane.If the ablation zone is located in the upper part of the right lobe, the lung is adjacent to the capsule.The air contained in the lung then reinforces the insulating property of this part of the capsule, which is modeled by imposing a Neumann condition.However, the influence of the capsule depends strongly on the proximity of the needles and vanishes when the distance increases (Fig. 6c).

How to provide clinically relevant simulations?
The above numerical sensitivity analysis shows that for deep-seated tumors, the numerical treatment planning solution, for the needles are exactly parallel and at the same depth cannot be used for the treatment efficiency assessment.Indeed, a small inclination of the needle due to anatomic constraint modifies substantially the isosurfaces of the electric field, and thus the electroporation ablation zone.In [20], we propose a numerical workflow to provide a numerical evaluation of the clinical electroporation ablation as performed by the interventional radiologist during the procedure.

Numerics for clinical electroporation ablation: radiologists at the center of the procedure
The novelty of this workflow is to mimic the steps of the clinical procedure, in order to the all the information provided by the clinical data.

Num. Step 1 Segmentation of the Regions of Interest (ROIs) from the preoperative image (CT-scan or MRI)
thanks to a semi-automatic segmentation before the day of the procedure.

Num.
Step 2-a Needles segmentation from the CBCT of the day of the procedure (see Fig. 2), and ROIs registration on the CBCT to get the patient-dependent geometrical framework.Num.
Step 3 Registration of the numerical simulation on the early postoperative MRI.

The issue of image registration
For the successful completion of the proposed numerical workflow, several concepts, such as delineation propagation, rely on establishing a spatial coherence between pre-/post-operative images and images acquired at different time instants over the course of the therapy.To meet this need, image-based motion estimation and compensation can relie on fast, automatic, accurate and precise registration algorithms.In the current study we used the evolution algorithm.The latter aims at minimizing a variationnal functional composed by a regularization term (assuming that displacements in neighboring voxels is moderate) added to a data fidelity term (i.e. a measure of the similarity between input images).The data fidelity term here quantifies the angle between the image gradient orientations.Parallel and anti-parallel gradients are favored, in order to match similar local contrast patterns.

Estimation of the conductivities from recorded intensities and feedback to physicians
Literature provides approximations of liver and tumor conductivities.However, as the conductivity parameters are patient-dependent, it is once again possible to use clinical data to adjust them.At the beginning of the procedure, test-pulses are recommended by the manufacturer, to check whether or not the device has been correctly setup.As the device provides current graphs with the intensities measured during the pulses, these tests can be advantageously used to estimate the conductivity parameters from the intensity measurements.
The values used for simulations are obtained thanks to a preliminary simulation, by a rough fitting of the simulated intensities on chronograms, as shown in Figure 7 by a trial error method -starting from conductivity values available in the literature and changing them handly step by step to obtain the fitting-, to fit the numerical intensities with the intensities recorded by Nanoknife R .The tolerance criterium of the fitting procedure is set such that during the first pulse of each pair combination, the relative error between the maximal value of the recorded intensity and the simulated intensity is less than 10%.It is worse noting that the calibration is performed once, with the initial needle location, and not modified afterwards.They are consistent with literature and are reported in Table 1.
The calibration enables to provide to physician the effective distribution of the electric field, as performed during the procedure (see Fig. 8).This new insights to the interventional radiologist could improve drastically Table 1.Estimate of conductivity parameters computed from test-pulse current graphs, from [20].
Parameter Value (mS cm −1 ) Reference the electroporation ablation, it provides to the physicians a visualization of the treated region, enabling possible treatment optimization.

Retrospective validation
The early postoperative MRI shows an impact of the treatment 3 days after the procedure.This area can be segmented and then, the images can be registered on the CBCT, for a retrospective study.Then, the preoperative ROIs and simulations can be superimposed on the registered MRI (Fig. 9) for comparison between observations and simulations.In Figure 9, the 2D area delineated by the isosurface 500 V cm −1 seems to match roughly with the contour of the observed area.In the 3D representation of Figure 10, the shape of the simulated area looks correct and is included in the observed area.This is consistent with the assumption that the observation contains the ablation area, and a margin in which irreversible and reversible phenomena coexist (margin between reversible and irreversible thresholds).Finally, it can be retained that the tumor is inside the area where the maximum electric field is beyond 500 V cm −1 , and inside the observed area, which gives a rather favorable tendency to the success of the procedure chosen for the proof of concept.
It should be noted that the mismatch between simulated and observed areas would indicate a model error or the omission of a crucial influence.As a result, even if what is observed on the MRI is not fully understood yet, this retrospective validation is essential to check the consistency of the simulation, and to give a second assessment of the procedure.

Discussion and still open questions
Each step of the numerical workflow raises crucial numerical questions which are described below and that need further investigations.

The importance of the image registration algorithm
Regarding the image registration, it is important to use non rigid registration algorithms.Indeed, it is worth noting that between two medical exams, the location of the liver may change.This is particularly true when needles are inserted: the insertion imposes a deformation to the liver, and therefore a rigid registration lead to non relevant results.In [20] we propose to use the EVolution algorithm proposed by Denis de Senneville et al. [14].Roughly speaking, the velocity field between the two gradient of the images is regularized by a diffusion term.The minimization is perfomed by solving the reaction-diffusion equation on the velocity field.This approach has been validated on clinical data, for soft tissue registration: bladder, liver, etc. However the use of constant in space diffusion coefficients is very restrictive.In particular, in a rigid body such a bone has to be registered, the algorithm fails to provide relevant results, because the bone would be deformed as the soft tissues.Therefore it is necessary to develop new strategies for clinical image registration.

The issue of the choice of boundary conditions
When several needles (more than 2) are used, the choice of the boundary conditions on the inactivated needles depends on the procedure.The interventional radiologist has the possibility to isolate the inactivated needles with sheathing.Then the standard Neumann boundary condition has to be imposed on the inactivated needles.However, due to a lack of time, the inactivated needles are sometimes left not insulated.The device imposes a high impedance so the net flux through the needle is zero, however Neumann condition is not appropriate to account for this situation.The floating potential boundary condition should be used: the potential needle is set at a constant value chosen such that the total net flux through the needle is zero [2].Then the isosurface of the eletric field are modified compared with the homogeneous Neumann conditions.Even though this situation often happens in clinical application, there is still a lack of appropriate numerical study of the influence of such boundary conditions.

Automatic model calibration procedure
The strategy of the model calibration lies in the fact that the model parameters impact the value of the electrical current that flows through the needle.Since the recording intensities are the only electrical data of tissue, it is natural to calibrate the model by fitting the numerical intensity on the recorded intensity.Here we present the results of [20] obtained thanks to a trial-error calibration.However it is important to develop an appropriate calibration strategy to improve and to automate this step.Several approaches will be tested in a near future, from simple Monte-Carlo approach, to more advanced calibration such as variational or sequential approaches.In particular, the sequential strategy of Chapelle, Moireau et al. [9,33,34] could be the most appropriate way to perform the joint state-parameter estimation.It which consists in solving a dynamical system obtained by incorporating in the original model a correction taking into account the discrepancy between the current simulation and the data.

Dynamical model of electroporation
Even though tissue electroporation has been discovered for several decades, there is currently a lack of model that can describe all the features of the phenomenon.Indeed, the model presented in equation (3.1) is standardly used, but is is a static model, which cannot described the dynamical behavior of the increase of tissue conductivity.Remark 6.1.We would like to make clear a misunderstanding which appeared in some previous studies.In some papers (see for instance [7,35,37]) this model is referred to as dynamic.We emphasize that the model is static since time is not involved in the equations.Due to the nonlinearity, the numerical solution is obtained iteratively, but the iterations have nothing to do with any time evolution.[31] a dynamical model of tissue electroporation that they use for clinical application [40].They adopted a numerical view point by introducing a discretized time step in the equations.However it seems difficult to write their model in a continuous way (see equation (11) of [31]) and their solution should hardly depend on their time step ∆t, which raises questions regarding the numerical stability.Voyer et al.

Langus et al. proposed in
proposed in [50] a dynamical two-phase model based on the electric circuit approach, which is written in a continuous way.The model can be written as follows: ∇ • (σ e ∇φ e + J c ) = 0, (6.1a) where E m = ∇φ e − σ −1 c J c , and The functions X 1 and X 2 are the respective degrees of poration and degrees of permeabilisation s described by Leguèbe et al. [32].The well-posedness of this model is still under investigation, but interestingly the model has been validated on experimental data [50].Fitting the dynamical model on the clinical data is a challenging numerical problem which will be addressed in forthcoming studies.

Conclusion
In this paper, we have shown how numerical modeling combined with clinical imaging can assist the clinical procedure of IRE ablation in liver tumor.
The 3D sensitivity analysis shows that the needle location is a crucial datum, which has to be precisely extracted from the clinical imaging to obtain realistic distribution of the elecric field.The proof of concept on a specific clinical case, which is detailed in [20] has been summarized here.It gives an interesting basis for confronting models to clinical cases in order to improve understanding and modeling.The electric field distribution in realistic geometric configuration (and for a simplistic model) shows that it possible to get an estimation of treatment success or failure from the computation.The IRENA software offers a sufficiently accurate computational framework for electric field computation from clinical data.In particular, this gives the possibility to the physician to visualize the estimated area of edema (see Sect. 1.1) during the procedure itself.Obviously, this simulation result has to be validated afterwards thanks to the postoperative imaging.But this is a first (and only) valuable information obtained in real-time, which can help the operator in his decision-making.Future implementations in IRENA will be related to new models, specific far field boundary conditions on the computational domain bounds, automatic calibration from test pulse current graphs, and a GUI for using the software in clinical routine.Forthcoming investigations will deal with more advanced modeling of tissue IRE, for which an energy-based approach could be adopted.As a conclusion, numerical modeling of IRE has a highly valuable potential impact in the clinical practice of IRE ablation.However to achieve this goal, several issues, from the high performance computing to the parameter estimation and from the mathematical modeling to the image registration have to be addressed and solved.

Statement of ethics board approval
This retrospective study is in accordance with ethical principals of the Declaration of Helsinki and has been approved by the local committee on human research of the University Hospital J. Verdier.

Appendix A. On the importance of 3D-simulations
Most of numerical studies are devoted to better understanding and modeling the electric field distribution during electroporation and its impact at cell or tissue level [5,13,16,27,29,30,36,42].For simplicity, these studies generally consider simple geometries, in two dimensions, with two parallel needles.
The disadvantage of the two-dimensional framework lies in the strong assumptions imposed on the geometric configuration: this supposes to consider electrodes parallel to the study plane, with infinite active lengths.though the generator manufacturer (NanoKnife) recommends parallel needles, this is not always possible in clinical practice, especially for treatment of liver deep-seated tumors.The presence of ribs, vena cava, hepatic arteries or important bile ducts may require a more complex needle placement.Then, if the needles are not parallel as in Figure A.1d, the 2D simulation can be very far from the results in three dimensions, and even may incorrectly suggest a total treatment of the tumor and therefore provide a wrong assessment of the procedure.
As clinical configurations are sometimes far from the ideal parallel configuration (see for instance Fig. A.2), a 3D framework is therefore mandatory.Note also that nonparallel needles do not compromise the clinical procedure.The distribution of the electric field is certainly less intuitive, but it is computationally predictable.This is one of the reasons why 3D numerical simulation is essential for clinical procedures.

Figure 1 .
Figure 1.Imaging of the liver.(Left): CTScan performed 1 month before the IRE procedure.The image exhibits a scar tissue left from a previous RF ablation.The patient biomarkers indicated a suspicion fo relapse located by the radiologist by the black arrow.(Right): Imaging of the liver of the same patient at Day3 after the IRE procedure in T2 weighed MRI modality.The treatment region exhibits a white enhancement for which radiology interpretation is still controversial[45].

Figure 2 .
Figure 2. CBCT Imaging of the liver of the same patient as Figure 1.The resolution is quite low but one can see easily the needles and the liver boundary.

Figure 3 .
Figure 3. Scheme of the modified ghost fluid method to address the issue of the needle location.The circle corresponds to the needle/tissue interface, the grid is the numerical grid.

Figure 4 .
Figure 4. Scheme of the intensity computation.The value of the potential u 0 is such that the isosurface C u0 is far enough from the needle.In practice, u 0 is 1/5 of the tension applied to the needle.

Figure 6 .
Figure 6.Numerical study of the sensitivity of the IRE threshold isolines (650 V cm −1 ) with respect to the geometry.(a) Influence of the needle position.The dotted (filled) lines stand for the positions 1 (2 resp.).(b) Influence of the tumor position.The dotted line corresponds to the tumor in the center of the needles.(c) Infuence of the hepatic capsule.The dotted line corresponds to the case without liver capsule.

Figure 7 .
Figure 7. Measured (left) and simulated (right) intensities of the first test-pulse of each pair of probes, from [20].

Figure 8 .
Figure 8. Isovalues of the maximum electric field magnitude.(a), (b) and (c) Simulated 3D isosurfaces 500 V cm −1 at each step of the procedure.

Figure 9 .
Figure 9. Superimposition of simulation results on the registered postoperative MRI.

Figure 10 .
Figure 10.3D comparison between the simulation (isosurface 500V cm −1 ) and the registered ROI of the treatment area (delineation in light color).

Figure A. 1
shows 2D and 3D simulations of the maximum electric field in a homogeneous tissue, after one pulse delivered by each pair of electrodes.For each 3D simulation, the area where the maximum electric field is above the IRE threshold is compared to this computed in the corresponding 2D simulation.If the 2D simulation is performed on the axial plane passing through the middle of the needle active tip, then the 2D simulation gives a result very close to what is observed on the same plane in the equivalent 3D simulation (Fig. A.1a).However, the results are different if the 2D simulation is performed in another plane (Fig. A.1b).Moreover, if the needles are parallel but inclined on the 2D simulation plane, then the results differ significantly as shown in Figure A.1c.These configurations impose the hypothesis of parallel needles, inserted at the same depth.Even

Figure A. 1 .
Figure A.1.Difference between 2D and 3D computation of the electric field magnitude after 1 pulse for each needle pair (V = 3000 V, σ 0 = 1 mS cm −1 , σ 1 = 4 mS cm −1 ).The dotted and filled lines stand for the IRE threshold (650 V cm −1 ) in 2D and 3D simulations, respectively.The median plane passes through the middle of the active part of the needles.

Figure A. 2 .
Figure A.2. Fluoroscopic image acquisition of the needles position.One can see that the needles are neither in exact parallel configuration nor at the exact same depth.