Home | About | Journals | Submit | Contact Us | Français |

**|**Biomed Eng Online**|**v.12; 2013**|**PMC3614452

Formats

Article sections

- Abstract
- Background
- Materials and methods
- Results
- Discussion
- Conclusions
- Competing interests
- Authors’ contributions
- References

Authors

Related links

Biomed Eng Online. 2013; 12: 16.

Published online 2013 February 21. doi: 10.1186/1475-925X-12-16

PMCID: PMC3614452

Selma Corovic,^{1} Igor Lackovic,^{2} Primoz Sustaric,^{3} Tomaz Sustar,^{3} Tomaz Rodic,^{3} and Damijan Miklavcic^{}^{1}

Selma Corovic: is.jl-inu.ef@civoroc.amles; Igor Lackovic: rh.ref@civokcal.rogI; Primoz Sustaric: is.m3c@zomirp; Tomaz Sustar: is.m3c@ratsus.zamot; Tomaz Rodic: is.m3c@cidor.zamot; Damijan Miklavcic: is.jl-inu.ef@cicvalkim.najimad

Received 2012 August 30; Accepted 2012 December 10.

Copyright © 2013 Corovic et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

Electroporation based therapies and treatments (e.g. electrochemotherapy, gene electrotransfer for gene therapy and DNA vaccination, tissue ablation with irreversible electroporation and transdermal drug delivery) require a precise prediction of the therapy or treatment outcome by a personalized treatment planning procedure. Numerical modeling of local electric field distribution within electroporated tissues has become an important tool in treatment planning procedure in both clinical and experimental settings. Recent studies have reported that the uncertainties in electrical properties (i.e. electric conductivity of the treated tissues and the rate of increase in electric conductivity due to electroporation) predefined in numerical models have large effect on electroporation based therapy and treatment effectiveness. The aim of our study was to investigate whether the increase in electric conductivity of tissues needs to be taken into account when modeling tissue response to the electroporation pulses and how it affects the local electric distribution within electroporated tissues.

We built 3D numerical models for single tissue (one type of tissue, e.g. liver) and composite tissue (several types of tissues, e.g. subcutaneous tumor). Our computer simulations were performed by using three different modeling approaches that are based on finite element method: inverse analysis, nonlinear parametric and sequential analysis. We compared linear (i.e. tissue conductivity is constant) model and non-linear (i.e. tissue conductivity is electric field dependent) model. By calculating goodness of fit measure we compared the results of our numerical simulations to the results of *in vivo* measurements.

The results of our study show that the nonlinear models (i.e. tissue conductivity is electric field dependent: σ(E)) fit experimental data better than linear models (i.e. tissue conductivity is constant). This was found for both single tissue and composite tissue. Our results of electric field distribution modeling in linear model of composite tissue (i.e. in the subcutaneous tumor model that do not take into account the relationship σ(E)) showed that a very high electric field (above irreversible threshold value) was concentrated only in the stratum corneum while the target tumor tissue was not successfully treated. Furthermore, the calculated volume of the target tumor tissue exposed to the electric field above reversible threshold in the subcutaneous model was zero assuming constant conductivities of each tissue.

Our results also show that the inverse analysis allows for identification of both baseline tissue conductivity (i.e. conductivity of non-electroporated tissue) and tissue conductivity vs. electric field (σ(E)) of electroporated tissue.

Our results of modeling of electric field distribution in tissues during electroporation show that the changes in electrical conductivity due to electroporation need to be taken into account when an electroporation based treatment is planned or investigated. We concluded that the model of electric field distribution that takes into account the increase in electric conductivity due to electroporation yields more precise prediction of successfully electroporated target tissue volume. The findings of our study can significantly contribute to the current development of individualized patient-specific electroporation based treatment planning.

Electroporation is an increase of cell membrane permeability to molecules that otherwise lack efficient transmembrane transport mechanisms. Due to membrane electroporation cell membrane conductivity is also increased. Electroporation can be achieved if cells are exposed to a sufficiently high electric field induced by application of high voltage pulses (i.e. electroporation pulses) [1-3]. The efficiency of electroporation strongly depends on pulse parameters such as pulse amplitude, pulse duration, number of pulses and repetition frequency [4]. In order to successfully electropermeabilize a cell membrane the induced transmembrane voltage needs to exceed a given value that depends on cell or tissue type (e.g. cell geometry, cell density, cell’s shape and orientation with respect to the electric field) [5-9]. It has been demonstrated that all types of cells (e.g. human, animal, plant and microorganisms) can be effectively electroporated, provided that appropriate electroporation parameters are selected for each of the cell or tissue type [10]. One of the key parameters for effective cell and tissue electroporation is however adequate electric field distribution within the treated sample [11]. The local electric field induces the transmembrane voltage on the cell [12,13], which is in turn correlated with transmembrane flux [3]. A cell or tissue is reversibly electroporated when the local electric field exceeds a reversible threshold value (E_{rev}) [14-16]. If the local electric field however exceeds an irreversible threshold (E_{irrev}) value the membranes of the exposed cells are disrupted to such an extent that cells fail to recover and they die (i.e. irreversible electroporation) [17-19].

Increase in cell membrane permeability allows for both administration of different molecules into electroporated cells (such as therapeutic drugs or molecular probes, which otherwise can not enter the cytosol) and for extraction of different cell components out from the cells [20,21]. Since, electroporation has been demonstrated to be a versatile method for both *in vitro* and *in vivo* settings it has been successfully applied to various domains such as medicine (e.g. human and veterinary oncology) [22-27], food processing [28] and microbial inactivation in water treatment [29].

Medical electroporation based therapies and treatments such as electrochemotherapy in local tumor control [27,30], tissue ablation method by means of irreversible electroporation [31,32] and gene therapy and DNA based vaccination by means of gene electrotransfer [33,34] are paving their way into clinical practice. In particular, electrochemotherapy has been demonstrated as an efficient tumor treatment as 85% objective response was obtained in treating skin metastasis of various cancers. In 2006 electrochemotherapy standard operating procedures (SOP) have been defined for the treatment of easily accessible cutaneous and subcutaneous tumor nodules of different histologies [35]. The SOP for electrochemotherapy were defined based on external parameters such as amplitude of applied voltage and distance between electrodes based on preclinical results and the results of numerical calculations in 3D models (1300 V/cm and 1000 V/cm for plate and needle electrodes respectively) [13,15]. The efficacy of SOP for easily accessible cutaneous tumors and subcutaneous tumor (depth: less than 2 cm) nodules were later confirmed by numerous clinical studies [30,36-38]. The *in vivo* studies that were based also on numerical simulations further demonstrated that the protocols of cutaneous tumors can be refined and that the electrochemotherapy can be successfully applied to deep-seated tumors based on analysis of local electric field distribution within the tumors [39,40]. In our previous study [41] we demonstrated that better electrochemotherapy outcome of cutaneous tumors can be obtained with increase of contact surface between plate electrodes and treated tissue. Similarly, Ivorra and colleagues [42] proposed the use of conductive gels in order to ‘homogenize’ the electric field within the treated tissue of easily accessible cutaneous tumors. However, for deep-seated tumors and tumors in internal organs, such as for example liver, brain or muscle tissue that are surrounded by other tissues having different electric properties an individualized patient-specific treatment planning is required as the heterogeneity of local electric field distribution within all exposed tissues can not be neglected [43,44]. Namely, several preclinical and clinical studies have shown that deep-seated tumors located in internal organs (breast cancer [27,37], bone metastasis [45], brain [46], liver metastases [40], melanoma metastasis in the muscle [39]) are also treatable. The outcome of the treatment however greatly depends on positioning of the electrodes, on applying of adequate pulses, electrode geometry and electrical properties of the treated tissues. These findings were either predicted or confirmed with numerical simulations, which resulted in development of computer-based treatment planning of both reversible electroporation based therapies (such as electrochemotherapy [43,44,47-49] and gene therapy and vaccination [50]) and irreversible electroporation based therapies [47,51-53]. Development of treatment planning of electroporation based therapies is currently focused on analysis of the shape and position of electrodes used and applied voltage to the electrodes [40,44,54,55]. Numerous studies demonstrated that the conductivity is increased due to membrane electroporation [2,56,57] and that in calculating the local electric field distribution within treated tissue one needs to account for these tissue conductivity increases [9,16,47,58,59]. In preclinical and clinical studies authors however often consider the treated tissues as being linear electric conductors (i.e. with constant tissue conductivities).

The aim of our study was to investigate whether the increase in electric conductivity of tissues needs to be taken into account when modeling tissue response to the electroporation pulses and how it affects the local electric field distribution within electroporated tissues. We also attempted to identify the nature of functional dependency of tissue conductivity on electric field. We therefore compared linear (i.e. tissue conductivity is constant) and non-linear (i.e. tissue conductivity is electric field dependent) model to the experimental data *in vivo*[56]. The study was performed for both single tissue (one type of tissue, e.g. liver) and composite tissue (several types of tissues, e.g. subcutaneous tumor). Our computer simulations were performed by using three different modeling approaches: 1. inverse analysis (IA) [60,61], 2. nonlinear parametric analysis (NPA) [62] and 3. sequential analysis (SA) [58]. The algorithm for IA modelling approach has been previously developed by [60,61] to model different physical phenomena in wide range of research and industry including biomechanics, pharmaceutical industry and space [63]. In this paper we introduced and successfully applied the IA approach for the first time in the field of modelling of biophysical processes that occur during tissue electroporation.

The results of all three modeling approaches are then compared to the results of our previous experimental measurements *in vivo*[56]. The results of our present study show that the nonlinear models fit experimental data better than linear models. We also demonstrate that tissue conductivity as a function of local electric field need to be taken into account when an electroporation based therapy or treatment is planned.

In order to investigate the change of electric properties of tissues on the electric field distribution during electroporation we built 3D numerical models of single tissue and composite tissue.

Modeling of electric field distribution in a single tissue during electroporation was carried out in liver tissue. We analyzed experimental data for plate and needle electrodes. Original experiments with plate electrodes were performed in rat liver, while the experiments with needle electrodes were performed in rabbit liver [56].

Modeling of electric field distribution in composite tissue (i.e. tissue composed of different types of tissues) was carried out in cutaneous tumor tissue. *In vivo* experiments were performed in mice with plate electrodes for two different types of cutaneous tumors: B16 and LPB sarcoma within the study by [56].

The high voltage pulses in all *in vivo* experiments were applied to the electrodes according to the electroporation protocol of eight 100 μs long pulses delivered to the tissues with 1 Hz repetition frequency. Within the experiments [56] the real time tissue electroporation control was performed by monitoring and measuring the electric current during the delivery of electric pulses. In order to evaluate the tissue electroporation efficiency during the pulse delivery the ^{51}CrEDTA indicator was used for all tissues. The system for the real time tissue electroporation control and ^{51}CrEDTA uptake measurements were described in our previous work [56]. The electric current measured during the *in vivo* experiments provides the information about the bulk tissue conductivity changes during the pulse delivery.

In order to compare and evaluate the results of our numerical models to the experimental data we calculated the goodness-of-fit measure R^{2} between numerically calculated and *in vivo* measured electric current I [A] at the end of the last delivered pulse (i.e. the 8-th pulse) to the tissue.

Electric field distribution in the models of tissue exposed to electroporation pulses U [V] can be determined by solving the equation (Eq. 1) for scalar electric potential. Namely, if we neglect the capacitive transient and time course of conductivity increase during the pulse, we may assume that the current density in tissue is divergence free and electric potential satisfies:

$$-\phantom{\rule{0.6em}{0ex}}\xb7\left(\sigma \xb7=\phantom{\rule{0.4em}{0ex}}0\right)$$

(1)

where σ and represent tissue conductivity [S/m] and electric potential [V], respectively.

Applied voltage (model input) was modeled as Dirichlet’s boundary condition on the contact surface between electrode and tissue geometry. For the model input values we used the amplitudes of the electroporation pulses applied *in vivo *[56]. In order to mathematically separate the conductive segment from its surroundings we applied Neuman’s boundary condition (J_{n}=0, where J_{n} is the normal electric current density [A/m^{2}]) on all outer boundaries of our models.

As long as the applied voltage U [V] is low enough the tissue conductivity σ [S/m] can be treated as a constant (i.e. σ=const. in Eq. 1), the problem is therefore described by the Laplace’s equation which is a linear partial differential equation that can be easily solved numerically. In this case the tissues can be modeled as a linear conductor with linear current–voltage I(U) relationships due to the constant tissue conductivities since the amplitude of the applied electroporation pulses is too low to produce an E above the reversible electroporation threshold (E<E_{rev}).

However, to account for experimentally observed tissue conductivity increase due to electroporation, Eq. 1 becomes nonlinear partial equation as σ is no longer constant. If the local electric field in the tissues exceeds the E_{rev} value, electric properties of treated tissue change (i.e. tissue conductivity increases due to the electroporation process). In this case the σ in Eq.1 depends on the electric field intensity E: σ=σ(E). Namely, during the application of electroporation pulses, tissue conductivity increases according to the functional dependency of the tissue conductivity on the local electric field distribution σ(E), which in our study describes the dynamics of the electroporation process. This subsequently results in nonlinear electric current voltage relationship I(U) (i.e. numerical model of tissue is non-linear).

Measurement of electric current I [A] during delivery of eight rectangular electroporation pulses of 100 μs duration and 1 Hz repetition showed that after rapid capacitive transient (which is due to tissue capacitance) the current increased towards the end of the pulse [56,52]. The rate of increase of electric current I [A] depends on the applied voltage U [V], which provides information about the changes in electric properties of the treated tissue (i.e. changes in electric conductivity of the tissue σ [S/m]). The rates of increase of I and σ thus allow us to identify the degree of tissue electroporation.

Numerical calculations of Eq.1 and tissue electroporation dynamics σ(E) in our study were performed with three different approaches that were based on finite element method: Inverse analysis (IA), nonlinear parametric analysis (NPA) and sequential analysis (SA). The software for IA was previously developed by Center for Computational Continuum Mechanics (Ljubljana, Slovenia) [61], while for the NPA and the SA we used Comsol Multiphysics software (version 3.5a, COMSOL, Sweden).

The three abovementioned modeling approaches are described here below:

The IA is based on standard Newton–Raphson method, which is widely used for calculating nonlinear problems based on finite element method. The tissue electroporation dynamics i.e. the dependence of the electrical conductivity with respect to electrical field strength σ(E) is described by linear and exponential relationships. Unknown parameters of these relationships are calculated by IA [60]. The IA is performed by using a symbolic system AceGen (Multi-language, Multi-environment Numerical Code Generation) [64] with AceFeM environment (The Mathematica Finite Element Environment) [65] for development of finite element equations and generation of corresponding finite element user subroutines. This enables symbolic formulation of the problem at a high abstract level. The derivation presented here follows a general approach to the automation of primal and sensitivity analysis of transient problems introduced by [66]. Below is a brief description of the formulation.

The response of the steady-state non-linear problem is defined by the global residual vector *R* and the vector of unknown variables of the problem (Eq. 2):

$$R\left(,=,*,{\displaystyle \sum _{e}{R}_{e}\left(\right)=0}\right)$$

(2)

where *R*_{e }represents the contribution of the *e-*th element to the global residual and *e* is the element number. The operator **Σ* is chosen here instead of the ordinary *Σ* symbol in order to denote that an assembly process takes place in which all element contributions have not only to be added up but also the kinematical compatibility between the elements has to be fulfilled [60]. The residual vector *R* is a part of standard notation used in finite element analysis. The finite element problems are represented by a set of equations in the residual form, which are solved iteratively using the Newton–Raphson method. The system given by the Eq. 2 can therefore be solved by the Newton–Raphson method, in which the following iteration is performed according to Eq. 3 and Eq. 4:

$$\frac{\mathit{dR}}{\mathit{d}}$$

(3)

$${\mathrm{\left(i+1\right)}}^{=}$$

(4)

The terms *R*(^{(i)}) and $\frac{\mathit{dR}}{\mathit{d}}$ are referred to as the residual (or load) vector and the tangent operator (or tangential stiffness matrix), respectively. The independent solution vector ^{(i)} is updated after the linear system given by the Eq. 3 has been solved for the unknown increment *Δ*. The upper index ^{(i)} refers to the quantities at the *i-*th iteration step of the outer loop.

The distribution of static electric field is determined by solving the Eq. 1. By using integral form and Galerkin approximation one can obtain

$$\underset{V}{\int}\mathit{\delta .\sigma .\mathrm{\phantom{\rule{0.12em}{0ex}}\mathit{dV}}-\phantom{\rule{0.4em}{0ex}}{\displaystyle \underset{S}{\int}J.n\phantom{\rule{0.12em}{0ex}}\mathit{\delta \phantom{\rule{0.12em}{0ex}}\mathit{dS}}=0}}$$

(5)

where is the electric potential, *V* is volume, *S* is surface and *J* is current density.

Discretization of Eq. 5 by =*N*_{e}._{e} and *δ*=*N*_{e}.*δ*_{e}=*N*_{e} leads to the Eq. 6:

$${R}_{e}\left(,=,\left[{\displaystyle \underset{{V}_{e}}{\int}{N}_{e}.\sigma .{\left({N}_{e},T\right)}^{\phantom{\rule{0.12em}{0ex}}}.{e}_{}}\right]\right)$$

(6)

where is the electric potential, _{e} are the electrical potentials in the element nodes, *δ* is the test function, *N*_{e }are the nodal shape functions and σ is the electrical conductivity.

We considered linear (Eq. 7) and exponential (Eq. 8) dependences of tissue conductivity *σ* with respect to electric field strength *E* as follows:

$$\sigma \left(E\right)=\left({\sigma}_{0}+{k}_{\sigma}\phantom{\rule{0.12em}{0ex}}E\right)\phantom{\rule{0.12em}{0ex}}\mathbf{I}$$

(7)

$$\sigma \left(E\right)=\left({\sigma}_{0}+{c}_{1}\phantom{\rule{0.12em}{0ex}}\left({e}^{{c}_{2}\phantom{\rule{0.12em}{0ex}}E}-1\right)\phantom{\rule{0.12em}{0ex}}\right)\mathbf{I}$$

(8)

where electric field strength *E* is defined as:

$$\begin{array}{l}E=-E=\sqrt{E.E}\end{array}$$

(9)

and **I** is the identity matrix. The constant in Eq. 7 and Eq. 8 stands for the baseline tissue conductivity of non-electroporated tissue (i.e. σ_{0}=σ(0)), while the constants *k*_{σ}, *c*_{1} and *c*_{2} stand for the parameters of linear (Eq. 7) and exponential (Eq. 8) functions describing electroporated tissue (i.e.σ(E)).

Unknown coefficients *σ*_{0}, *k*_{σ}, *c*_{1} and *c*_{2} were determined with IA [60] by minimizing the objective function *f* of the following form:

$$f\left({\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right)={\displaystyle \sum _{j}{\left(\frac{{I}_{\mathit{FEM}}-{I}_{\mathit{\text{measured}}}}{U}\right)}_{\left(j\right)}^{2}}$$

(10)

Determination of unknown coefficients *σ*_{0}, *k*_{σ}, *c*_{1} and *c*_{2} allowed us to detect the properties of the treated tissue for both situations σ=const. (i.e. σ_{0}=σ(0)) and σ(E)).

By using the Newton's method in the Mathematica (version 8.0) programming environment the minimization problem is defined as:

$$g\left(x=\left\{{\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right\}\right)=\left\{\begin{array}{l}\frac{f\left({\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right){\sigma}_{0}={\displaystyle \sum _{j}{\left(2\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}-{I}_{\mathit{\text{measured}}}}{{U}^{2}}\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}{\sigma}_{0}}{}\left(j\right)\right)}_{},}\frac{f\left({\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right){k}_{\sigma}={\displaystyle \sum _{j}{\left(2\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}-{I}_{\mathit{\text{measured}}}}{{U}^{2}}\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}{k}_{\sigma}}{}\left(j\right)\right)}_{}}\frac{f\left({\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right){c}_{1}={\displaystyle \sum _{j}{\left(2\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}-{I}_{\mathit{\text{measured}}}}{{U}^{2}}\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}{c}_{1}}{}\left(j\right)\right)}_{}}\frac{f\left({\sigma}_{0},{k}_{\sigma},{c}_{1},{c}_{2}\right){c}_{2}={\displaystyle \sum _{j}{\left(2\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}-{I}_{\mathit{\text{measured}}}}{{U}^{2}}\phantom{\rule{0.12em}{0ex}}\frac{{I}_{\mathit{FEM}}{c}_{2}}{}\left(j\right)\right)}_{}}}{}=0}{}}{}}{}\end{array}\right\}$$

(11)

where *U* is the voltage between the electrodes, *I*_{FEM} is the calculated current and *I*_{measured }the measured current in the *j*-th measurement. The derivatives $\frac{{I}_{\mathit{FEM}}{\sigma}_{0}}{}$, $\frac{{I}_{\mathit{FEM}}{k}_{\sigma}}{}$, $\frac{{I}_{\mathit{FEM}}{c}_{1}}{}$ and $\frac{{I}_{\mathit{FEM}}{c}_{2}}{}$ were calculated within the finite element method analysis.

The NPA has been developed in our previous study [62]. In the present study we used NPA to solve the equation Eq. 1 for both conditions σ=const. (for E<E_{rev}) and σ=σ(E) (for E≥E_{rev}) by using Comsol Multiphysics’ stationary nonlinear solver. Namely, we treat the problem as stationary since we analyzed the tissue behavior by calculating the σ(E) dependency at the end of the applied pulse U [V]. The parameters defined by NPA are the values of the applied voltages U (i.e. model inputs). The electric current I [A] is calculated as a model output providing the value of I [A] at the end of pulse. The Comsol Multiphysics’ nonlinear solver uses an affine invariant form of damped Newton method. The user must supply an initial guess for the dependent variable U (i.e. U_{0}). Starting with the initial guess, the software forms the linearized model using U_{0} as the linearization point i.e. the discrete form of the equations as f(U)=0, where f(U) is the residual vector and U is the solution vector. The iterative Newton method can be applied in the usual way, expanding the residual vector in a Taylor series, neglecting all but the first order term. The software solves the discretized form of the linearized model f ’(U_{0})δU=−f(U_{0}) for the Newton step δU by using the selected linear system solver (the f ’(U_{0}) is the Jacobian matrix). It then computes the new iterate U_{1}=U_{0}+λδU, where λ is the damping factor. The modified Newton correction then estimates the error E_{error} for the new iterate U_{1} by solving f’(U_{0})E_{error}=−f(U_{1}). If the relative error corresponding to E is larger than the relative error in the previous iteration, the code reduces the damping factor λ and recomputes U_{1}. This algorithm repeats the damping-factor reduction until the relative error is less than in the previous iteration or until the damping factor underflows the minimum damping factor. When it has taken a successful step U_{1}, the algorithm proceeds with the next Newton iteration. A value of λ=1 results in Newton’s method, which converges quadratically if the initial guess U_{0} is sufficiently close to the solution. Comsol Multiphysics will reduce the damping factor and re-compute U_{1} if the relative error is larger than the previous value. The purpose of this damping factor is to make the process converge for a broader range of initial values (initial guesses U_{0}).

It is necessary to increase the voltage U from 0 V in small steps and use the solution from a previous step as the initial condition for the next step. This is taken care of automatically by the parametric solver. The parametric solver consists of a loop around the stationary solver. If the nonlinear solver does not converge it tries a smaller parameter step; if it does converge it determines the size of the next parameter step based on the speed of the convergence of the Newton iteration.

The SA in our study was based on the sequential permeabilization model proposed by [58], where changes in tissue conductivity were used as an indicator of tissue permeabilization. Namely, the dynamics of electroporation was modeled as a discrete process with the sequence of static finite element models. Each of the finite element models described the process in one discrete interval (each of the discrete intervals relates to a real, discrete, but undetermined time interval). In each static finite element model in sequence, the tissue conductivity was determined based on the electric field distribution from the previous model in the sequence according to the Eq. 12:

$$\sigma \phantom{\rule{0.25em}{0ex}}\left(k\right)=f\phantom{\rule{0.25em}{0ex}}\left(E\left(k-1\right)\right)$$

(12)

where *k* stands for the number of static finite element models in sequence.

Model input is the applied voltage pulse, and model outputs are the electric field distribution E and total electric current I in each specific sequence k. The increase in electrical current I from I_{0} to I_{k} simulates the tissue response during the delivery of the electroporation pulses U in each discrete interval k (static finite element model in sequence) to the tissue electroporation (i.e. due to the functional dependency of σ on the electric field distribution E). If the electroporation does not occur, σ remains constant, thus I=I_{0}. More detailed description of SA is given in [16,58,59].

The numerical calculations obtained with NPA and SA by using COMSOL Multiphysics 3.5a software package were performed in the 3D Conductive Media DC application mode on a computer running Windows XP (64 bit) with a 3.00 GHz Intel Pentium D processor and 2 GB of RAM.

For the numerical calculations with IA we used AceGen (version AceFEM 3.101) system [61,64,65]. The calculations were performed on a computer running Windows XP (64 bit) with 2.67 GHz Intel Core i7 CPU and 6 GB of RAM.

Results of numerically calculated model outputs (such as total current and local electric field distribution) with all three approaches were controlled for numerical errors by increasing the mesh density until the output results changed by less that 0.5%. The mesh of the resulting finite element models of liver tissue electroporated with plate and needle electrodes were composed of approximately 30 000 and 25 400 elements, respectively. Number of generated finite elements of the tumor mesh model was 30 268.

The 3D geometry of the numerical model and its 2D cross-sectional view in XY plane of a single tissue (i.e. liver tissue) is given in Figure 1. The model of liver tissue electroporated with plate electrodes is given in Figure 1A, while the model of liver tissue electroporated with a pair of needle electrodes is given in Figure 1B.

By using the IA numerical simulations we identified both baseline electric conductivity of tissue σ_{0} (for the condition σ(0)) and linear and exponential σ(E) (Eqs. 7 and 8) that fitted best the experimental data *in vivo*[56]. The determined σ_{0} and σ(E) for plate and needle electrodes were then compared to the data from the literature.

Since the NPA and SA required the baseline tissue conductivity σ_{0}, the conductivity increase factor due to electroporation, the reversible and ireversible threshold values (E_{rev} and E_{irrev}) to be predefined, in our models we selected the data from [58,62,67,68]. The functional dependencies σ(E) predefined in the models were sigmoid-like function implemented as smoothed Heaviside's step function [62] with continuous second derivative and sigmoid function σ(E) that was implemented by Sel and colleagues [58].

The geometry of composite tissue i.e. subcutaneous tumor is shown in Figure 2A. The distance between two plate electrodes used in experiments was d=5.2 mm. The tumor geometry was created so as to match as accurate as possible the electrode placement and the geometry of the tissue treated within the *in vivo* experiments. The electrodes in Figure 2A only sketch the site of electrodes’ position in the experiments and do not take part of the numerical model. The body of the electrodes was not modeled in our numerical model in order to reduce the computation time of numerical simulations. The electrodes were modeled rather as boundary conditions defined on the electrode-tissue contact surfaces (Figure 2B).

The numerical model of the tumor was built as a heterogeneous 3D model composed of four types of tissues with different electric properties (i.e. electric conductivities). The skin in our subcutaneous tumor model was modeled as two-layer tissue: the first layer comprised the electric properties of stratum corneum and epidermis (σ_{SE}), while the second layer comprised dermis and fat tissue (σ_{DF}). The underlying muscle tissue was modeled as anisotropic tissue with σ_{M} (Y and Z direction -Figure 2) and σ_{M||} (X direction - Figure 2) in perpendicular and parallel direction of E with respect to the muscle fibers. The baseline conductivities of the modeled tissues (i.e. the initial σ of non-electroporated tissues for E<E_{rev}), the conductivity increase factor due to tissue electroporation and the threshold values E_{rev} and E_{irrev} for each of the tissues were selected based on data we published previously [9,16,59,69]and based on the comparison of I [A] calculated in our models to the *in vivo* measurements. Our numerical and experimental analysis was performed for two different types of cutaneous tumors: B16 melanoma and LPB sarcoma.

The relationship of electric conductivity with respect to the local electric field σ(E) for all modeled tissues was considered to be smoothed Heaviside step function taking into account the electric conductivities and the threshold values data listed in Table 1. For the numerical calculations we used NPA. We examined the electric current and local electric field distribution within the model of subcutaneous tumor for three different electric properties of target tumor tissue σ_{T1}, and σ_{T2} and σ_{T3}.

To compare the outputs of our models to the experimentally measured data we calculated coefficient of determination R^{2}, which represents a statistical measure of how well the model data approximates the real data points (i.e. the measure of the goodness-of-fit between the modeled and measured data). The R^{2}=1 means that the model data perfectly fit the real data points. If the R^{2} has a value closer to the 0 or a negative value the model data do not fit the real data points.

The coefficient of determination is given by the Eq.13.

$${R}^{2}=1-\frac{S{S}_{\mathit{err}}}{S{S}_{\mathit{tot}}}$$

(13)

where *SS*_{err} and *SS*_{tot} are residual sum of squares and the total sum of squares, respectively.

The "variability" of the data set is measured through the sums of squares *SS*_{err} and *SS*_{tot}.

The residual sum of squares *SS*_{err} is defined by the Eq. 14:

$$S{S}_{\mathit{err}}={{\displaystyle \sum _{i}\left({y}_{i}-{f}_{i}\right)}}^{2}$$

(14)

where *y*_{i} data stand for the measured values of electric current at the end of the pulse, each of which has an associated modeled value *f*_{i}.

The total sum of squares SS_{tot}, which is proportional to the sample variance, is given by the Eq. 15:

$$S{S}_{\mathit{tot}}={{\displaystyle \sum _{i}\left({y}_{i}-\stackrel{\_}{y}\right)}}^{2}$$

(15)

where the $\stackrel{\_}{y}$ is the mean of the measured data *y*_{i}, which is given by the Eq. 16:

$$\stackrel{\_}{y}=\frac{1}{n}{\displaystyle \sum _{i}^{n}{y}_{i}}$$

(16)

where *n* is the number of observations.

Modeling of electric field distribution during electroporation with plate electrodes was performed with IA and with NPA for the model geometry given in Figure 1A. The geometry details are identical as in *in vivo* experiments [56] and in our previous work where the NPA was first implemented [62].

By using IA we first identified linear and exponential relationship σ(E) (Eqs. 7 and 8) and determined corresponding baseline tissue conductivity σ(0) (Figure 3A).

The calculated I(U) and G(U) for both identified relationships σ(E) (linear and exponential) by means of IA are shown in Figure 4. In Figure 4 the I(U) and G(U) characteristics with σ(E) were also compared to the corresponding I(U) and G(U) obtained in the model with constant electric conductivity (i.e. with baseline tissue conductivities: σ=const.=0.093 S/m and σ=const.=0.088 S/m).

We then performed numerical modeling with NPA with smoothed Heaviside relationships of σ(E) for two predefined baseline electrical conductivities: 0.126 S/m [67] and 0.091 S/m [68]. Within the smoothed Heaviside relationships of σ(E) the conductivity increase factor was chosen to increase by factor of two and three [62] which is in agreement with previously experimentally observed σ increase at the end of pulse [2] For the reversible and irreversible threshold values we have chosen 200 V/cm and 1000 V/cm, respectively [62]. The examined smoothed Heaviside relationships σ(E) with corresponding baseline conductivities and the conductivity increase factors are shown in Figure 3B.

The obtained I(U) and G(U) characteristics in the models with relationships σ(E) and constant σ (from Figure 3B) are shown in Figure 5A and B, respectively.

Calculated goodness of fit R^{2} measure for the results obtained with IA and NPA are given in Tables 2 and and3,3, respectively. Almost identical fit was obtained for exponential and linear σ_{L}(E) relationship with IA (Table 2). The results of NPA performed in the model with smoothed Heaviside relationship σ_{L}(E) and baseline conductivity 0.126 S/m better fitted the experimental data when the conductivity increase factor was chosen to increase by a factor of two - σ_{L1}(E), while the model with baseline conductivity 0.093 S/m better fitted the experimental data when the conductivity increase factor was chosen to increase by a factor of three - σ_{L2}(E), Table 3.

In order to identify the baseline liver conductivity and functional dependency of the tissue conductivity on local electric field σ(E) we first performed numerical modeling by using inverse analysis IA on a set of experimental data obtained with needle electrodes with diameter 0.7 mm (Figure 6B). The identified relationship σ(E) was determined to be linear with baseline liver conductivity σ_{0}=0.046 S/m (Figure 7B). In order to verify the identified parameters of IA we further implemented the same identified relationship on sets of experimental data obtained with needle electrodes with diameters: 0.3 mm, and 1.1 mm, given in Figure 6A and C, respectively and calculated the goodness of fit measure R^{2} for all three models (Table 4). The numerically calculated I(U) relationship taking into account the linear relationship σ(E) was compared to the I(U) with σ=constant corresponding to the identified baseline tissue conductivity (σ=const. = 0.046 S/m) (Figure 6A, B and C).

In order to compare the results of NPA and SA to the experimental data we also calculated I(U) relationships for predefined relationships σ(E): smoothed Heaviside with NPA and sigmoid with NPA and SA (Figure 7B). We predefined the baseline liver conductivity, reversible and irreversible threshold values based on previously reported data by Sel and colleagues [58,54]: σ_{L0}=0.067 S/m, E_{rev}=460 V/cm and E_{irrev}=700 V/cm. The NPA and SA were also performed for tissue model with needle electrodes of 1.1 mm diameter given in Figure 1B. Our model corresponds to the one used in [58] in which the SA-analysis was first developed and implemented. The I(U) relationships obtained with σ(E) were compared to the models with constant σ (i.e. σ=0.067 S/m and σ=0.046 S/m), Figure 7.

Comparison of U(I) characteristics obtained in the model (needle electrodes, diameter 1.1 mm) with relationships σ(E) and σ=const. to the experimental data by using IA, NPA and SA is shown in Figure 7.

Calculated goodness of fit measure R^{2} for models with needle electrodes obtained with all three different analyses is given in Table 4.

By using SA analysis we calculated the I(U) characteristics for two different discretizations of modeled pulses: 5 steps and 15 steps. In Figure 7 we showed only the I(U) for 5 discretization steps since the goodness-of-fit for I(U) obtained with 15 steps did not significantly differ from the goodness-of-fit obtained with 5 steps: R^{2}=0.8958 (15 steps) vs. R^{2}=0.8938 (5 steps), Table 4.

The dashed green curve and the solid green curve in the Figure 7A represent the I(U) relationships obtained with SA and NPA, respectively. In both analyses the same σ(E) was used (i.e. sigmoid function represented with green solid curve in Figure 7B). The difference between the two I(U) relationships occurs due to the different discretisation approaches of σ(E). Namely, the discretisation of σ(E) used in SA is usually coarser (i.e. discretisation of σ(E) is obtained by employing 5 step functions) compared to finer automatic discretization of σ(E) used in NPA.

The comparison of numerically calculated total current I [A] in tumor models (with σ(E) and σ=const. given in Table 1) to the total current measured *in vivo* for two different types of tumors B16 melanoma and LPB sarcoma is given in Figure 8A and B, respectively. The examined σ_{T}(E) are shown in Figure 9B.

The calculated goodness of fit R^{2} for two different types of tumors B16 melanoma and LPB sarcoma is given in Table 5.

The calculated volume fraction (i.e. volume percentage V [%]) of the target tumor tissue exposed to the local the electric field distribution above E_{rev} and E_{irrev} for three different σ(E) relationships of target tumor tissue (i.e. σ_{T1}(E), σ_{T2}(E) and σ_{T3}(E)) and for σ independent of E (i.e. σ_{T}=const.), is presented in Figure 9. The entire volume of the target tumor tissue (V_{T_E>Erev}=100%) was exposed to the local electric field above E_{rev} at different applied voltages for different target tumor tissue conductivities σ_{T}(E): U=340 V (for σ_{T1}(E)), U=400 V (for σ_{T2}(E)) and U=440 V (for σ_{T3}(E)). The exposure of the target tumor tissue to the *E* above irreversible threshold value V_{T_E > Eirrev} was also obtained for different U for different σ_{T}(E): V_{T_E >Eirrev}=100% for σ_{T1}(E) was obtained at U=720 V Figure 9A, while the V_{T_E >Eirrev} =100% for σ_{T2}(E) and σ_{T3}(E) was obtained for U>800 V (in our study the U was applied from 20 V to 800 V).

The volume of the target tumor tissue V_{T} in the subcutaneous model with constant conductivities of all tissue σ_{SE}, σ_{DF}, σ_{M} and σ_{M||} and σ_{T}=const. was zero (V_{T_E>Erev}=0) as shown in Figure 9A.

We visualized the calculated local electric field distribution within the model of subcutaneous tumor tissue for two different applied voltages U=176 V and U=276 V (Figure 10) in XY cross-section plane (see Figure 2). At the applied voltage U=176 V the local electric field in the target tumor tissue just started to enter the tumor, while the U=276 V was chosen to show the difference in distribution of E for different electric properties of the target tumor tissue σ_{T}(E).

The electric field distribution is calculated for four different conductivities of target tumor tissue σ_{T}: Figure 10A σ_{T}=const., Figure 10B σ_{T1}(E), Figure 10C σ_{T2}(E) and Figure 10D σ_{T3}(E). The electric field distribution in Figure 10 was displayed in the range from the reversible E_{rev} and irreversible threshold value E_{irrev} of the target tumor tissue 400 V/cm≤E≤800 V/cm.

In a recent study of robustness of electrochemotherapy treatment planning based on analysis of local electric field distribution E the authors show that the uncertainties in predefined dielectric properties of the treated tissues and in the rate of increase in electric conductivity due to the electroporation have large effect on treatment effectiveness, which indicated that more experimental and numerical research is needed [43]. Numerous other numerical and experimental studies also showed that the electric conductivity changes during tissue electroporation. However, in preclinical studies and clinical studies the electroporated tissues are often considered as conductive materials with constant electric conductivities.

In our study we therefore investigated whether the functional dependency of electric properties and of local electric field distribution within treated tissues need to be taken into account when modeling tissue response to the electroporation pulses.

We compared linear (i.e. tissue conductivity is constant σ=const.) and non-linear model (tissue conductivity depends on local electric field σ(E)) to the experimental data *in vivo *[56]. We built our numerical models so as to fit as accurately as possible the tissue geometry, the electrode geometry and the contact surface between the electrodes and tissue treated in *in vivo* experiments. In the first part of our study we performed the comparison for single tissue (i.e. one type of tissue) for two types of electrodes: plate and needle electrodes. For single tissue modeling we analyzed the experimental data performed in rat liver treated with plate electrodes and in rabbit liver treated with needle electrodes (Figure 1).

The main findings of the single tissue analysis were then applied to the model of composite tissue (i.e. the tissue composed of different types of tissues with different electric conductivities). The study of composite tissue was performed for subcutaneous tumor that consisted of stratum corneum, dermis, epidermis, connective and fatty tissue, tumor tissue and muscle tissue and was treated with plate electrodes (Figure 2).

Numerical simulations in our study were performed by using three different modeling approaches: IA [60,61], NPA [62] and SA [58].

In our study the IA was first applied for numerical modeling of electroporation in a single tissue (i.e. rat liver tissue) treated with plate electrodes. The relationship between local electric field distribution and the conductivity of liver tissue σ(E) was described by linear (Eq. 7) and exponential functions (Eq. 8). The parameters of these two functions (parameters σ_{0} and k_{σ} of linear function and parameters σ_{0}, c_{1} and c_{2} of exponential function) were than calculated with IA-analysis algorithm by fitting the data of numerically calculated to the experimentally measured values of total electric current flowing through the tissue (Figure 3). The parameters σ_{0} in both functions σ(E) (linear and exponential) represented the baseline tissue conductivities for the condition σ_{0}=σ(0)) (i.e. the initial conductivity of non-electroporated tissue). We calculated these parameters to be σ_{0}=0.088145 S/m (linear function Eq. 7) and σ_{0}=0.092768 S/m (exponential function Eq. 8). The numerically calculated values of σ_{0} of both functions correspond well to the measured conductivities of rat liver tissue published by Gabriel and colleagues [68]. These results show that the IA allows for identification of the initial tissue conductivity based on measured values of total electric current flowing through the tissue exposed to the electroporation pulses. In order to compare the numerical models of liver tissue that consider exponential and linear σ(E) relationship to the numerical model with constant conductivity (σ(E)=const.) we used the calculated data σ_{0} for development of the models with constant conductivities: σ(E)=σ_{0}=const. The calculation of goodness of fit measure (R^{2}) of the numerical models that take into account the relationship σ(E) fit better experimental data compared to the models with constant conductivity: linear σ(E) R^{2}=0.7993 and exponential σ(E) R^{2}=0.7995 vs. σ(E)=σ_{0}=0.088145 S/m=const. R^{2}=0.1295 and σ(E)=σ_{0}=0.092768 S/m=const. R^{2}=0.1297 (Table 2)). These results show that the IA also allows for identification of functional dependency of tissue conductivity σ(E) of the electroporated tissue. The results of our study thus demonstrate that by using IA both the baseline electric conductivity for non-electroporated tissue (for the condition σ(0)=const. for E=0 or E<E_{rev}) and the σ(E) relationship of the tissue exposed to U that results in tissue electroporation (E>E_{rev}) can be identified.

The calculation of tissue conductance G [S] show that the conductance of the models with σ(E) increased with increase of the applied voltage following the measured data of G Figure 3D. We can observe that the models G=const. do not describe the electroporation process of the tissue which is in agreement with our previous findings [2,56], as well as recent multiscale modeling [57].

In order to compare numerical simulations obtained by nonlinear parametric analysis (NPA) to the numerical simulations obtained with inverse analysis (IA) we used the same geometry of liver tissue and electrodes (Figure 1A) and electric properties of tissue (σ_{L0}=const. and σ(E)) as used in our previous work [62]. Namely, for the initial parameters of tissue conductivity we used two different measured values from the literature σ_{L0}=0.091 S/m (pig liver) [68] and σ_{L0}=0.124 S/m (rat liver) [67]. The relationship between the local electric field and tissue conductivity σ(E) was implemented to the tissue model as smoothed Heaviside’s function. The numerical calculations were performed for both constant conductivity σ=const. and σ(E) with conductivity increase factor of two and three (Figure 3B), which is in agreement with experimentally measured conductivity increase at the end of the pulse [2,58]. Our results show that the numerical models with σ(E) better fit experimental data than models with σ=const. (Table 3). Numerical model with initial conductivity of rat liver better fit experimental data for the conductivity increase factor of two (R^{2}=0.7997), the numerical model with initial conductivity of pig liver better fit experimental data for the conductivity increase factor of three (R^{2}=0.7838). When comparing the numerical results obtained with NPA to the results obtained with IA we can conclude that similar results and good agreement between numerical and experimental data can be obtained with both modeling approaches when σ(E) relationship is taken into account (regardless of the shape of the σ(E)). Both approaches show that the models with σ(E) fit the measured data considerably better than models with constant conductivity. However, more precise measurements are needed for determination of the shape of the σ(E) function that might be tissue type specific.

In the second part of our study we built numerical models of single tissue (rabbit liver) treated with needle electrodes of three different diameters 0.3 mm, 0.7 mm and 1.1, that were used also in *in vivo* experiments. By using IA we first identified the σ(E) for the tissue model with needle electrodes with diameter of 0.7 mm, as in this case we analyzed the maximum number of measurements Figure 6. We found that the best fit of numerical data with experimental data was obtained with linear σ(E) function (Figure 6B) with R^{2}=0.875605 (Table 3). The identified value of the initial tissue conductivity was σ_{0}=0.04593 S/m which correspond to the measured conductivity of rabbit liver (summarized in [58]). The identified linear function in the model with 0.7 mm was then applied to the models with 0.3 mm and 1.1 mm and good agreement between calculated data were obtained R^{2}=0.92873 (for 0.3 mm) and R^{2}=0.883214 (for 1.1 mm) (Table 4).

For the model with electrode diameter 1.1 mm we also performed NPA with smoothed Heaviside and sigmoid σ (E) function and sequential analysis with sigmoid function σ(E) with parameters determined by [58]. We also performed numerical analysis with σ(E)=const. The comparison of the model with σ(E)=const. to the numerical models with σ(E) showed that models with σ(E) fit experimental data considerably better then models with σ(E)=const. (Table 4). When comparing numerical results of different modeling approaches (Figure 7) we further observed that the best fit between the numerical and experimental data was obtained with NPA with smoothed Heaviside’s function R^{2}=0.9192 (Table 4).

Based on the comparison of all three modeling approaches validated with experimental measurements we determined the important advantages for electroporation process modeling of each of the approaches. Namely, the IA can be used as a modeling approach for the situations when the conductivity of treated tissue is difficult to be measured before or during electroporation. The NPA allows for parametric analysis of electric current at the end of the applied voltage pulses while the SA allows also an insight into the electroporation process (i.e. the course of σ(E) changes) during the applied pulse. Although, the NPA and SA modeling approaches require the baseline conductivity and the σ(E) need to be predefined, the important advantage is that these two modeling approaches are relatively simple to be used in commercial software which is accessible and widely used.

In the third part of our study we built the numerical model of subcutaneous tumor in order to examine the influence of tissue heterogeneity on the electric field distribution during electroporation of a composite tissue (i.e. tumor model composed of different types of tissues). We compared the numerical calculation of total electric current in the tumor model with constant conductivity of all constituent tissues and the numerical calculation of total electric current in the tumor model that took into account σ(E) of all involved tissues to the experimental measurements *in vivo* carried out in two different types of tumors: B16 melanoma (Figure 8A) and LPB sarcoma (Figure 8B). For all tissues in the model we implemented smoothed Heaviside’s function σ(E).

Based on the comparison of the numerical and experimental data *in vivo* we demonstrated that the model that took into account σ(E) of all involved tissues fits better *in vivo* data than model with constant conductivities (Table 5). Our results further show that the subcutaneous tumor model with higher conductivity of target tumor tissue σ_{T3}(E) fits better the *in vivo* data of B16 melanoma while the model with σ_{T1}(E) better fits the *in vivo* data of LPB sarcoma subcutaneous tumor (Table 5), suggesting that B16 melanoma has higher conductivity than LBP sarcoma. The visualization of local E in the subcutaneous tumor model demonstrates that the more conductive the target tumor tissue the higher the intensity of E in the surrounding tissue and the lower E in the tumor tissue (Figure 10A and B), which is in agreement with our previous findings [70].

We also calculated the tumor volume exposed to the local electric field above E_{rev} and E_{irrev} (Figure 9A) and demonstrated that more conductive target tumor tissue were more difficult to be successfully electroporated, since they require higher voltage U to be applied to the electrodes in order to expose the entire volume of the target tumor above the E_{rev} or higher, which also results in higher intensity of E in the skin and muscle and subsequently may induce more damage to the surrounding healthy tissue (Figure 10). We also demonstrated that the electroporation of more conductive tissues resulted in higher values of total electric current I [A] flowing through tissue (Figure 8). The exposure of the target tumor tissue to the E above irreversible threshold value also depends on the electric properties of the target tumor tissue σ_{T}(E). It needs to be emphasized that the subcutaneous model with constant conductivities of all tissues (i.e. in the subcutaneous tumor model that do not take into account the relationship σ(E)) showed that a very high electric field (E>E_{irrev}) was concentrated only in the stratum corneum while the target tumor tissue was not treated (Figure 10A). Furthermore, the volume of the target tumor tissue V_{T} exposed to the E>E_{rev} in the subcutaneous model with constant conductivities of all tissue was zero (Figure 9A).

In our study we investigated whether the increase in electric conductivity of tissues needs to be taken into account when modeling tissue response to the electroporation pulses and how it affects the local electric distribution within electroporated tissues. We examined electric field distribution during electroporation in both linear models (i.e. tissue conductivity is constant) and non-linear models (i.e. tissue conductivity is electric field dependent). The results show that nonlinear models fit experimental data better than linear models. We demonstrated that the tissue conductivity as a function of electric field (i.e. σ(E)) needs to be considered when modeling tissue behavior during electroporation. Therefore, the σ(E) relationship needs to be taken into account when an electroporation based treatment is planned or investigated. The findings of our study can be of great importance for precise treatment planning of electroporation based therapies and treatments (e.g. for electrochemotherapy, gene electrotransfer for gene therapy and DNA vaccination, tissue ablation with irreversible electroporation and transdermal drug delivery). Particularly, for an electroporation based therapy of deep-seated cutaneous tumors and tumors in internal organs, such as for example bone, brain or muscle tissue that are surrounded by other tissues having different electric properties an individualized patient-specific treatment planning that takes into account heterogeneity of electric conductivity and local electric field distribution within all exposed tissues is required.

The authors declare that they have no competing interests.

SC and DM designed the study. PS, TS and TR performed modeling with IA. IL and SC performed modeling with NPA. SC performed modeling with SA. All authors were involved in the analysis of numerical and experimental data. All authors were involved in the preparing of the manuscript. All authors read and approved the final manuscript.

This work was supported by the Slovenian Research Agency and conducted within the scope of LEA EBAM (European Associated Laboratory on the Electroporation in Biology and Medicine). This manuscript is a result of the networking efforts of the COST Action TD1104 (http://www.electroporation.net). The authors thank Dr. David Cukjati for the results of *in vivo* experiments performed within his postdoctoral project.

- Neumann E, Schaefer-Ridder M, Wang Y, Hofschneider PH. Gene transfer into mouse lyoma cells by electroporation in high electric fields. EMBO J. 1982;1:841–845. [PubMed]
- Pavlin M, Kanduser M, Rebersek M, Pucihar G, Hart FX, Magjarevic R, Miklavcic D. Effect of cell electroporation on the conductivity of a cell suspension. Biophys J. 2005;88:4378–4390. doi: 10.1529/biophysj.104.048975. [PubMed] [Cross Ref]
- Kotnik T, Pucihar G, Miklavcic D. Induced transmembrane voltage and its correlation with electroporation-mediated molecular transport. J Membrane Biol. 2010;236:3–13. doi: 10.1007/s00232-010-9279-9. [PubMed] [Cross Ref]
- Pucihar G, Krmelj J, Rebersek M, Napotnik TB, Miklavcic D. Equivalent Pulse Parameters for Electroporation. IEEE T Biomed Eng. 2011;58:3279–3288. [PubMed]
- Teissié J, Rols MP. An experimental evaluation of the critical potential difference inducing cell membrane electropermeabilization. Biophys J. 1993;65(1):409–413. doi: 10.1016/S0006-3495(93)81052-X. [PubMed] [Cross Ref]
- Cemazar M, Jarm T, Miklavcic D, Macek Lebar A, Ihan A, Kopitar NA, Sersa G. Effect of electric-field intensity on electropermeabilization and electrosensitivity of various tumor-cell lines in vitro. Electro Magnetobiol. 1998;17:263–272.
- Valic B, Golzio M, Pavlin M, Schatz A, Faurie C, Gabriel B, Teissie J, Rols MP, Miklavcic D. Effect of electric field induced transmembrane potential on spheroidal cells: theory and experiments. Eur Biophys J. 2003;32:519–528. doi: 10.1007/s00249-003-0296-9. [PubMed] [Cross Ref]
- Pucihar G, Kotnik T, Valic B, Miklavcic D. Numerical determination of transmembrane voltage induced on irregularly shaped cells. Annals Biomed Eng. 2006;34:642–652. doi: 10.1007/s10439-005-9076-2. [PubMed] [Cross Ref]
- Corovic S, Zupanic A, Kranjc S, Al Sakere B, Leroy-Willig A, Mir LM, Miklavcic D. The influence of skeletal muscle anisotropy on electroporation: in vivo study and numerical modeling. Med Biol Eng Comput. 2010;48:637–648. doi: 10.1007/s11517-010-0614-1. [PMC free article] [PubMed] [Cross Ref]
- Prud’homme GJ, Glinka Y, Khan AS, Draghia-Akli R. Electroporation enhanced nonviral gene transfer for the prevention or treatment of immunological, endocrine and neoplastic diseases. Curr Gene Ther. 2006;6:243–273. doi: 10.2174/156652306776359504. [PubMed] [Cross Ref]
- Miklavcic D, Beravs K, Semrov D, Cemazar M, Demsar F, Sersa G. The importance of electric field distribution for effective in vivo electroporation of tissues. Biophys J. 1998;74:2152–2158. doi: 10.1016/S0006-3495(98)77924-X. [PubMed] [Cross Ref]
- Kotnik T, Bobanovic F, Miklavcic D. Sensitivity of transmembrane voltage induced by applied electric fields – a theoretical analysis. Bioelectrochem Bioenerg. 1997;43:285–291. doi: 10.1016/S0302-4598(97)00023-8. [Cross Ref]
- Miklavcic D, Corovic S, Pucihar G, Pavselj N. Importance of tumour coverage by sufficiently high local electric field for effective electrochemotherapy. Eur J Cancer. 2006;4(Suppl):45–51.
- Rols MP, Teissié J. Electropermeabilization of mammalian cells, Quantitative analysis of the phenomenon. Biophys J. 1990;58(5):1089–1098. doi: 10.1016/S0006-3495(90)82451-6. [PubMed] [Cross Ref]
- Miklavcic D, Semrov D, Mekid H, Mir LM. A validated model of in vivo electric field distribution in tissues for electrochemotherapy and for DNA electrotransfer for gene therapy. Biochim Biophys Acta. 2000;1523:73–83. doi: 10.1016/S0304-4165(00)00101-X. [PubMed] [Cross Ref]
- Corovic S, Mir LM, Miklavcic D. In vivo muscle electroporation threshold determination: realistic numerical models and in vivo experiments. J Membrane Biol. 2012;245:509–520. doi: 10.1007/s00232-012-9432-8. [PubMed] [Cross Ref]
- Davalos RV, Mir LM, Rubinsky B. Tissue ablation with irreversible electroporation. Ann Biomed Eng. 2005;33:223. doi: 10.1007/s10439-005-8981-8. [PubMed] [Cross Ref]
- Rubinsky J, Onik G, Mikus P, Rubinsky B. Optimal parameters for the destruction of prostate cancer using irreversible electroporation. J Urol. 2008;180:2668–2674. doi: 10.1016/j.juro.2008.08.003. [PubMed] [Cross Ref]
- Al-Sakere B, Andre F, Bernat C, Connault E, Opolon P, Davalos RV, Rubinsky B, Mir LM. Tumor ablation with irreversible electroporation. PlosOne. 2007;2:1135. [PMC free article] [PubMed]
- Vorobiev E, Lebovka N. Electrotechnologies for Extraction from Food Plants and Biomaterials. New York: Springer Science; 2008.
- Sack M, Sigler J, Frenzel S, Chr E, Arnold J, Michelberger T, Frey W, Attmann F, Stukenbrock L, Mueller G. Research on Industrial-Scale Electroporation Devices Fostering the Extraction of Substances from Biological Tissue. Food Eng Rev. 2010;2:147–156. doi: 10.1007/s12393-010-9017-1. [Cross Ref]
- Cemazar M, Tamzali Y, Sersa G, Tozon N, Mir LM, Miklavcic D, Lowe R, Teissie J. Electrochemotherapy in veterinary oncology. J Vet Intern Med. 2008;22(4):826–831. doi: 10.1111/j.1939-1676.2008.0117.x. [PubMed] [Cross Ref]
- Testori A, Tosti G, Martinoli C, Spadola G, Cataldo F, Verrecchia F, Baldini F, Mosconi M, Soteldo J, Tedeschi I, Passoni C, Pari C, Di Pietro A, Ferrucci PF. Electrochemotherapy for cutaneous and subcutaneous tumor lesions: a novel therapeutic approach. Dermatol Ther. 2010;23(6):651–661. doi: 10.1111/j.1529-8019.2010.01370.x. [PubMed] [Cross Ref]
- Cemazar M, Jarm T, Sersa G. Cancer electrogene therapy with interleukin-12. Curr Gene Ther. 2010;10(4):300–311. doi: 10.2174/156652310791823425. [PubMed] [Cross Ref]
- Heller LC, Heller R. Electroporation gene therapy preclinical and clinical trials for melanoma. Curr Gene Ther. 2010;10(4):312–317. doi: 10.2174/156652310791823489. [PubMed] [Cross Ref]
- Neal RE II, Rossmeisl JH Jr, Garcia PA, Lantz O, Henao-Guerrero N, Davalos RV. Successful treatment of a large soft-tissue sarcoma with irreversible electroporation. J Clin Oncol. 2011;29(13):372–377. doi: 10.1200/JCO.2010.33.0902. [PubMed] [Cross Ref]
- Sersa G, Cufer T, Paulin Kosir S, Cemazar M, Snoj M. Electrochemotherapy of chest wall breast cancer recurrence. Cancer Treat Rev. 2012;38(5):379–386. doi: 10.1016/j.ctrv.2011.07.006. [PubMed] [Cross Ref]
- Morales-de La Peña M, Elez-Martínez P, Martín-Belloso O. Food preservation by pulsed electric fields: an engineering perspective. Food Eng Rev. 2011;3:94–107. doi: 10.1007/s12393-011-9035-7. [Cross Ref]
- Rieder A, Schwartz T, Schön-Hölz K, Marten SM, Süß J, Gusbeth C, Kohnen W, Swoboda W, Obst U, Frey W. Molecular monitoring of inactivation efficiencies of bacteria during pulsed electric field (PEF) treatment of clinical wastewater. J Appl Microbiol. 2008;105:2035–2045. doi: 10.1111/j.1365-2672.2008.03972.x. [PubMed] [Cross Ref]
- Campana LG, Valpione S, Mocellin S, Sundararajan R, Granziera E, Sartore L, Chiarion-Sileni V, Rossi CR. Electrochemotherapy for disseminated superficial metastases from malignant melanoma. BJS. 2012;99:821–830. doi: 10.1002/bjs.8749. [PubMed] [Cross Ref]
- Pech M, Janitzky A, Wendler JJ, Strang C, Blaschke S, Dudeck O, Ricke J, Liehr UB. Irreversible electroporation of renal cell carcinoma: a first-in-man phase I clinical study. Cardiovasc Intervent Radiol. 2011;34(1):132–138. doi: 10.1007/s00270-010-9964-1. [PubMed] [Cross Ref]
- Onik G, Rubinsky B. In: Irreversible Electroporation. Rubinsky B, editor. Heidelberg, Germany: Springer Berlin; 2010. Irreversible electroporation: First patient experience—Focal therapy of prostate cancer; pp. 235–247.
- Daud AI, DeConti RC, Andrews S, Urbas P, Riker AI, Sondak VK, Munster PN, Sullivan DM, Ugen KE, Messina JL, Heller R. Phase I trial of interleukin-12 plasmid electroporation in patients with metastatic melanoma. J Clin Oncol. 2008;26:896–903. [PMC free article] [PubMed]
- Ahmad S, Casey G, Sweeney P, Tangney M, O'Sullivan GC. Optimized electroporation mediated DNA vaccination for treatment of prostate cancer. Genet Vaccines Ther. 2010;8:1. doi: 10.1186/1479-0556-8-1. [PMC free article] [PubMed] [Cross Ref]
- Mir LM. Bases and rationale of the electrochemotherapy. EJC. 2006;4(Suppl):38–44.
- Campana LG, Mocellin S, Basso M, Puccetti O, De Salvo GL, Chiarion-Sileni V, Vecchiato A, Corti L, Rossi CR, Nitti D. Bleomycin-based electrochemotherapy: clinical outcome from a single institution's experience with 52 patients. Ann Surg Oncol. 2009;16(1):191–199. doi: 10.1245/s10434-008-0204-8. [PubMed] [Cross Ref]
- Matthiessen LW, Chalmers RL, Sainsbury DC, Veeramani S, Kessell G, Humphreys AC, Bond JE, Muir T, Gehl J. Management of cutaneous metastases using electrochemotherapy. Acta Oncol. 2011;50(5):621–629. doi: 10.3109/0284186X.2011.573626. [PMC free article] [PubMed] [Cross Ref]
- Gargiulo M, Papa A, Capasso P, Moio M, Cubicciotti E, Parascandolo S. Electrochemotherapy for non-melanoma head and neck cancers: clinical outcomes in 25 patients. Ann Surg. 2012;255(6):1158–1164. doi: 10.1097/SLA.0b013e31824f68b2. [PubMed] [Cross Ref]
- Miklavcic D, Snoj M, Zupanic A, Kos B, Cemazar M, Kropivnik M, Bracko M, Pecnik T, Gadzijev E, Sersa G. Towards treatment planning and treatment of deep-seated solid tumors by electrochemotherapy. Biomed Eng Online. 2010;9:10. doi: 10.1186/1475-925X-9-10. [PMC free article] [PubMed] [Cross Ref]
- Edhemovic I, Gadzijev EM, Brecelj E, Miklavcic D, Kos B, Zupanic A, Mali B, Jarm T, Pavliha D, Marcan M, Gasljevic G, Gorjup V, Music M, Pecnik Vavpotic T, Cemazar M, Snoj M, Sersa G. Electrochemotherapy: A new technological approach in treatment of metastases in the liver. Technol Cancer Res Treat. 2011;10:475–485. [PubMed]
- Corovic S, Al-Sakere B, Haddad V, Miklavcic D, Mir LM. Importance of contact surface between electrodes and treated tissue in electrochemotherapy. Technol Cancer Res Treat. 2008;7:393–400. [PubMed]
- Ivorra A, Al-Sakere B, Rubinsky B, Mir LM. Use of conductive gels for electric field homogenization increases the antitumor efficacy of electroporation therapies. Phys Med Biol. 2008;53:6605–6618. doi: 10.1088/0031-9155/53/22/020. [PubMed] [Cross Ref]
- Kos B, Zupanic A, Kotnik T, Snoj M, Sersa G, Miklavcic D. Robustness of treatment planning for electrochemotherapy of deep-seated tumors. J Membrane Biol. 2010;236:147–153. doi: 10.1007/s00232-010-9274-1. [PubMed] [Cross Ref]
- Pavliha D, Kos B, Zupanic A, Marcan M, Sersa G, Miklavcic D. Patient-specific treatment planning of electrochemotherapy: Procedure design and possible pitfalls. Bioelectrochemistry. 2012;87:265–273. [PubMed]
- Fini M, Tschon M, Alberghini M, Bianchi G, Mercuri M, Campanacci L, Cavani F, Ronchetti M, De Terlizzi F, Cadossi R. In: Clinical aspects of electroporation. Kee ST, Gehl J, Lee EW, editor. New York: Springer; 2011. Cell electroporation in bone tissue; pp. 115–127.
- Agerholm-Larsen B, Iversen HK, Ibsen P, Moller JM, Mahmood F, Jensen KS, Gehl J. Preclinical validation of electrochemotherapy as an effective treatment for brain tumors. Cancer Res. 2011;71(11):3753–3762. doi: 10.1158/0008-5472.CAN-11-0451. [PubMed] [Cross Ref]
- Garcia PA, Pancotto T, Rossmeisl JH, Henao-Guerrero N, Gustafson NR, Daniel GB, Robertson JL, Ellis TL, Davalos RV. Non-Thermal Irreversible Electroporation (N-TIRE) and Adjuvant Fractionated Radiotherapeutic Multimodal Therapy for Intracranial Malignant Glioma in a Canine Patient. Technol Cancer Res Treat. 2011;10:73–83. [PubMed]
- Corovic S, Zupanic A, Miklavcic D. Numerical modeling and optimization of electric field distribution in subcutaneous tumor treated with electrochemotherapy using needle electrodes. IEEE T Plasma Sci. 2008;36:1665–1672.
- Zupanic A, Corovic S, Miklavcic D. Optimization of electrode position and electric pulse amplitude in electrochemotherapy. Radiol Oncol. 2008;42:93–101.
- Zupanic A, Corovic S, Miklavcic D, Pavlin M. Numerical optimization of gene electrotransfer into muscle tissue. Biomed Eng Online. 2010;9:66. doi: 10.1186/1475-925X-9-66. [PMC free article] [PubMed] [Cross Ref]
- Zupanic A, Miklavcic D. In: Irreversible Electroporation. Rubinsky B, editor. Heidelberg: Springer; 2010. Optimization and numerical modeling in irreversible electroporation treatment planning; pp. 203–222.
- Neal RE II, Garcia PA, Robertson JL, Davalos RV. Experimental characterization and numerical modeling of tissue electrical conductivity during pulsed electric fields for irreversible electroporation treatment planning. IEEE T Biomed Eng. 2012;59(4):1077–1085. [PubMed]
- Garcia PA, Rossmeisl JH, Neal RE, Ellis TL, Olson JD, Henao-Guerrero N, Robertson J, Davalos RV. Intracranial Nonthermal Irreversible Electroporation: In Vivo Analysis. J Membr Biol. 2010;236:127–136. doi: 10.1007/s00232-010-9284-z. [PubMed] [Cross Ref]
- Mahmood F, Gehl J. Optimizing clinical performance and geometrical robustness of a new electrode device for intracranial tumor electroporation. Bioelectrochemistry. 2011;81:10–16. doi: 10.1016/j.bioelechem.2010.12.002. [PubMed] [Cross Ref]
- Zupanic A, Kos B, Miklavcic D. Treatment planning of electroporation-based medical interventions: electrochemotherapy, gene electrotransfer and irreversible electroporation. Phys Med Biol. 2012;57:5425–5440. doi: 10.1088/0031-9155/57/17/5425. [PubMed] [Cross Ref]
- Cukjati D, Batiuskaite D, André F, Miklavčič D, Mir LM. Real time electroporation control for accurate and safe in vivo non-viral gene therapy. Bioelectrochemistry. 2007;70:501–507. doi: 10.1016/j.bioelechem.2006.11.001. [PubMed] [Cross Ref]
- Essone Mezeme M, Pucihar G, Pavlin M, Brosseau C, Miklavcic D. A numerical analysis of multicellular environment for modeling tissue electroporation. Appl Phys Lett. 2012;100:143701. doi: 10.1063/1.3700727. [Cross Ref]
- Sel D, Cukjati D, Batiuskaite D, Slivnik T, Mir LM, Miklavcic D. Sequential finite element model of tissue electropermeabilization. IEEE Trans Biomed Eng. 2005;52:816–27. doi: 10.1109/TBME.2005.845212. [PubMed] [Cross Ref]
- Pavselj N, Bregar Z, Cukjati D, Batiuskaite D, Mir LM, Miklavcic D. The course of tissue permeabilization studied on a mathematical model of a subcutaneous tumor in small animals. IEEE Trans Biomed Eng. 2005;52:1373. doi: 10.1109/TBME.2005.851524. [PubMed] [Cross Ref]
- Gresovnik I. A General Purpose Computational Shell for Solving Inverse and Optimisation Problems - Applications to Metal Forming Processes. PhD Thesis: University of Wales Swansea, UK; 2000.
- Center for Computational Continuum Mechanics (Ljubljana, Slovenia) [ http://www.c3m.si]
- Lackovic I, Magjarevic R, Miklavcic D. In: Proceedings of IFMBE (XII Mediterranean Conference on Medical and Biological Engineering and Computing) Magjarevic R, Bamidis Panagiotis D, Pallikarakis N, editor. Berlin Heidelberg: Springer; 2010. Incorporating Electroporation-related conductivity changes into models for the calculation of the electric field distribution in tissue; pp. 695–698.
- Giavazzi S, Ganatea MF, Trkov M, Sustarsic P, Rodic T. Inverse determination of viscoelastic properties of human fingertip skin. Materials and Geoenvironment. 2010;57:1–16.
- Korelc J. AceGEN (Multi-language, Multi-environment Numerical Code Generation) – Users manual. 2009. [ http://www.fgg.uni-lj.si/symech/]
- Korelc J. AceFEM ((The Mathematica Finite Element Environment) – Users manual. 2009. [ http://www.fgg.uni-lj.si/symech/]
- Korelc J. Automation of primal and sensitivity analysis of transient coupled problems. Comput Mech. 2009;44:631–649. doi: 10.1007/s00466-009-0395-2. [Cross Ref]
- Haemmerich D, Staelin ST, Stai JZ, Tungjitkusolmun S, Mahvi DM, Webster JG. In vivo electrical conductivity of hepatic tumours. Physiol Meas. 2003;24:251–260. doi: 10.1088/0967-3334/24/2/302. [PubMed] [Cross Ref]
- Gabriel C, Peyman A, Grant EH. Electrical conductivity of tissue at frequencies below 1 MHz. Phys Med Biol. 2009;54:4863–4878. doi: 10.1088/0031-9155/54/16/002. [PubMed] [Cross Ref]
- Pavselj N, Miklavcic D. A numerical model of permeabilized skin with local transport regions. IEEE T Biomed Eng. 2008;55:1927–1930. [PubMed]
- Corovic S, Pavlin M, Miklavcic D. Analytical and numerical quantification and comparison of the local electric field in the tissue for different electrode configurations. Biomed Eng Online. 2007;6:37. doi: 10.1186/1475-925X-6-37. [PMC free article] [PubMed] [Cross Ref]

Articles from BioMedical Engineering OnLine are provided here courtesy of **BioMed Central**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |