 Research
 Open Access
 Published:
Anistropically varying conductivity in irreversible electroporation simulations
Theoretical Biology and Medical Modelling volume 14, Article number: 20 (2017)
Abstract
Background
One recent area of cancer research is irreversible electroporation (IRE). Irreversible electroporation is a minimally invasive procedure where needle electrodes are inserted into the body to ablate tumor cells with electricity. The aim of this paper is to propose a mathematical model that incorporates a tissue’s conductivity increasing more in the direction of the electrical field as this has been shown to occur in experiments.
Method
It was necessary to mathematically derive a valid form of the conductivity tensor such that it is dependent on the electrical field direction and can be easily implemented into numerical software. The derivation of a conductivity tensor that can take arbitrary functions for the conductivity in the directions tangent and normal to the electrical field is the main contribution of this paper. Numerical simulations were performed for isotropicvarying and anisotropicvarying conductivities to evaluate the importance of including the electrical field’s direction in the formulation for conductivity.
Results
By starting from previously published experimental results, this paper derived a general formulation for an anistropicvarying tensor for implementation into irreversible electroporation modeling software. The anistropicvarying tensor formulation allows the conductivity to take into consideration both electrical field direction and magnitude, as opposed to previous published works that only took into account electrical field magnitude.
The anisotropic formulation predicts roughly a five percent decrease in ablation size for the monopolar simulation and approximately a ten percent decrease in ablation size for the bipolar simulations. This is a positive result as previously reported results found the isotropic formulation to overpredict ablation size for both monopolar and bipolar simulations. Furthermore, it was also reported that the isotropic formulation overpredicts the ablation size more for the bipolar case than the monopolar case. Thus, our results are following the experimental trend by having a larger percentage change in volume for the bipolar case than the monopolar case.
Conclusions
The predicted volume of ablated cells decreased, and could be a possible explanation for the slight overprediction seen by isotropicvarying formulations.
Background
Electroporation is an electricallydriven biological process that uses externally applied electric fields to briefly open the nanosized pores in a cell’s membrane. On one hand, electroporation can be used to transport drugs and genetic materials into cells [1]. On the other hand, a strong enough electrical field can irreversibly damage the cell’s membrane and cause cell death. This process is called irreversible electroproation (IRE), and has gained interest as a possible tumor ablation method [2].
A cell’s viability is directly related to the electrical field it is exposed to. Each cell is characterized by a critical value so that when the cell is exposed to an electrical field above this critical value it loses its ability to reseal the pores in its membrane and eventually dies. Any electrical field below this critical value, either does not form pores or the pores close when the electrical field is removed. The critical value depends on various factors such as tissue type, tumor stage, electrode configuration and size, and the direction of the applied electric field. The success of an IREbased therapy for tumors relies on this critical value and thus knowledge about the optimal configuration and size of the electrodes as well as the applied voltage needed to completely ablate the tumor while minimizing the damage to healthy tissue is relevant to clinicians.
Mathematical models of electroporation and corresponding numerical simulations able to predict the electric field distribution in a tissue of interest for various geometries and positions of electrodes could play an important, complementary role in the design of individual treatment protocols [3, 4]. In particular, it is well established in the literature that electroporation can affect the conductivity of tissue [5, 6]. The large increase in conductivity is a result of the membrane’s pores opening during electroporation [7]. Corovic et al. [8] recommended to incorporate in mathematical models a conductivity that is dependent on the electrical field strength to increase accuracy. Experiments performed by Mezeme et al. suggest that electroporation does not always isotropically increase the conductivity of the tissue, but conductivity is increased more in the direction tangent to the electric field [9]. Their experiments determined the conductivity of electroporated chicken livers by using magentic resonance electrical impedance tomography. The result of an anistropic conductivity tensor can be explained by the pores forming in the cell membrane with a bias towards the direction of the electrical field [10]. This is because the voltage drop across the cell membrane is not equal around the entire cell, but is largest where the cell membrane is perpendicular to electrical field [10].
Therefore, in this paper we propose a mathematical model that incorporates a tissue’s conductivity increasing more in the direction of the electrical field. To do so, it was required for us to derive a formulation for the conductivity tensor such that conductivity increases by different amounts in the direction tangent and normal to the electrical field and can be easily implemented in numerical methods such as finite elements. Numerical simulations were performed for isotropic and anisotropic varying conductivities to evaluate the importance of including the electrical field’s direction in the formulation for conductivity.
Method
This section will discuss the details of the mathematical models used for comparing the well established isotropically increasingconductivity formulation and this paper’s anistropically increasingconductivity formulation. The simulations used for comparison will be based off the irreversible electroporation experiments and simulations performed by Castellvi et al. [11]. This work will use their experimentally determined parameters to compare an isotropic formulation to an anisotropic formulation.
Governing equations
In numerical simulations for treatment planning, it is assumed that the applied direct current electric field is at steady state [12]. The electrical potential, U, satisfies the differential form of Gauss’s law at steady state:
where σ is the tissue’s conductivity. Then the electric field E is given by: E=−∇U. If σ is a constant, then Eq. (1) becomes the wellknown Laplace equation. However, the conductivity of tissue is not in general constant but increases during electroporation [6], and therefore, in this paper, we model the conductivity as a function of the electric field, σ=σ(E).
We model the isotropic varying conductivity as a sigmoid Gompertz curve [8, 13, 14]:
where A and B are unitless coefficients, σ _{0} is the the base conductivity before electroporation, σ _{ max } is the maximum conductivity a tissue can achieve after electroporation, and ∥E∥ is the l ^{2}− norm of the electrical field. By definition, the l ^{2}− norm of the electrical field is \(\E \ = \sqrt {E_{x}^{2} + E_{y}^{2} + E_{z}^{2} }\) where E _{ x }, E _{ y }, E _{ z } is the x, y, z component of the electrical field respectively. The values of A and B are found by fitting the curve for conductivity to experimental data.
If the conductivity is isotropic, then σ is a scalar. However, the conductivity is anisotropic when the conductivity increases more in the direction of the electrical field. This results in the conductivity, σ, being represented by a rank2 tensor.
We will compare two different cases. The first being when σ does not take into account the direction of the electrical field (isotropicvarying) and when σ does take into account the direction of the electrical field (anisotropicvarying). Both cases will solve Eq. 1; the only difference will be in how σ is defined.
Isotropicvarying formulation
This is the well established model that is often used for IRE simulation predictions. It will be used for comparison with our anisotropicvarying formulation. The isotropicvarying formulation has the conductivity increase equally in all directions. Hence, it being refereed to as the isotropic formulation in this paper. The advantage of using an isotropic formulation is that the conductivity is represented by a single scalar, and the direction of the electrical field does not need to be considered. It is also less computationally intensive.
In the isotropic formulation, the changes in conductivity from electroporation can be taken into account with a sigmoid Gompertz curve for the conductivity [8, 14, 15]. The same expression for conductivity as was found in [11] which is
will be used for the isotropic formulation.
Anisotropicvarying tensor derivation
In this paper, we propose a conductivity that takes into account the electrical field’s direction and magnitude. We wish to formulate the conductivity tensor such that the conductivity increases more in the direction of the electrical field. Since the electrical field does not increase equally in all directions, we will refer to this formulation as the anisotropic formulation in this paper.
We now wish to derive the matrix representation for the conductivity tensor for the anisotropicvarying case such that it can be implemented into a numerical scheme. This is the main contribution of the paper.
We begin by assuming the conductivity tensor is a real valued symmetric matrix. This is common assumption and is justified by assuming conductivity is independent of current flow being positive or negative. All symmetric matrices are diagnolizable. Therefore, the conductivity tensor is diagonalizable.
The anisotropic case can have a different conductivity in the direction of the electrical field than perpendicular to it. If we make the assumption that the conductivity in the direction of the electrical field is always equal or greater than any other direction, and that the conductivity in a direction perpendicular to the electrical field is always equal or less than any nonnormal direction then the directions tangent and normal to the electrical field are principal directions. These assumptions are justified according to the experimental work found in [9].
Let σ _{ t } be an arbitrary function representing the conductivity in the direction tangent to the electrical field. Also let σ _{ n } and σ _{ b } be arbitrary functions representing conductivity in the direction normal and binormal to the electrical field respectively. By treating σ _{ t }, σ _{ n }, and σ _{ b } as arbitrary functions, we are providing a derivation that is valid for any expression for σ _{ t }, σ _{ n }, and σ _{ b } which may someday be ascertained through experiments.
A tensor is diagonal when the principal directions are chosen as the basis representation. We will therefore consider a coordinate system comprised of a unit vector tangent to the electrical field and two unit vectors perpendicular to the electrical field. A coordinate system consisting of a tangent vector, normal vector, and binormal vector, is sometimes referred to as a FrenetSerret frame. These three vectors form an orthonormal basis in which the basis vectors point in the principal directions of the conductivity. Since the FrenetSerret basis vectors are aligned with the principal directions then the conductivity tensor matrix representation is diagonal in the FrenetSerret coordinate system and can be represented as
where σ _{ f } is the matrix representation of the conductivity tensor in the FrenetSerret coordinate system. It will be assumed that the conductivity is the same for any direction normal to the electrical field. This implies that both the normal and binormal direction conductivities are equal,
Equation 5 simplifies Eq. 4 to
where σ _{ f } represents a linear transformation that takes a vector from the FrenetSerret frame and outputs a vector in the FrenetSerret frame. This can be restated as
where \(\mathfrak {R}_{f}^{3}\) is the space of three dimensional vectors represented by the FrenetSerret basis. However, for calculations, we would like to have the conductivity tensor in the Cartesian coordinate system as it is often impractical to implement the FrenetSerret matrix representation into a numerical scheme such as finite elements. Thus, we wish to derive
where σ _{ c } is the matrix representation of the conductivity tensor in the Cartesian coordinate system, and \(\mathfrak {R}^{3}_{c}\) is the space of three dimensional vectors represented by the standard Cartesian basis. We will use the known form of the matrix representation of the conductivity tensor in the FrenetSerret frame to derive the matrix representation of the conductivity tensor in the Cartesian frame.
Since the FrenetSerret frame and the Cartesian coordinate system are both orthonormal, a change of basis between the two can be represented by a rotation through an angle θ about the zaxis and a rotation ϕ about the yaxis. Note only two rotations instead of three are necessary because the conductivity has been assumed to be equal in all normal directions.
The rotation tensor, R(θ,ϕ), is a linear transformation such that
A function diagram outlining the relationship between σ _{ c }, σ _{ f }, and R can be found in Fig. 1.
We can decompose R into two subsequent rotations. The first rotation will be about the zaxis, R _{ z }(θ), and can be expressed by the matrix,
where the angle θ is obtained through the relationship
The quantities E _{ x } and E _{ y } represent the x and y components of the electrical field respectively. Similarly, the second rotation will be about the yaxis, R _{ y }(ϕ), and is expressed by the matrix,
where ϕ is the angle described by the function
The quantity E _{ z } represents the z component of the electrical field. The rotation tensor can be then be expressed as
Using properties of rotation matrices, R ^{−1} is
We are now able to express the conductivity tensor in matrix form with the domain and codomain being the Cartesian frame. This is accomplished by applying the rotation matrices in the following order:
where
The form of an anisotropic varying conductivity tensor has been derived with arbitrary functions for the conductivity in the tangent and normal directions. Therefore, this formulation is valid for any functions of conductivity for σ _{ t } and σ _{ n }, including those dependent on electrical field strength.
For the anisotropic case simulations, the conductivity used for the tangent direction, σ _{ t }, will have the same form as the conductivity for the isotropic case:
This was chosen because the measurements used to experimentally determine the form of the conductivity is often done using electrodes that are aligned with the electrical field.
As shown in [9], the conductivity tensor becomes more anisotropic as the electrical field strength increases. This ratio will be modelled using a sigmoid function, σ _{ Δ }, where
is an adaptation of the data obtained in the experiments reported in [9]. Due to the limited amount of experimental data, this was the authors’ best guess as to the form for σ _{ Δ }. But since the derivation was done in a general setting, it will be possible to plug in different formulations at a later time as they become available from experimental data. The approximation can be combined with the previous expression for σ _{ t } to obtain the following expression for the conductivity in the perpendicular direction:
It is worth noting, that our form for σ _{ Δ } results in σ _{ t } and σ _{ n } being equal until the onset of electroporation. In other words, the conductivity is isotropic until the onset of electroporation. During electroproation, the conductivity tensor becomes more anisotropic with increasing electrical fields up until a value of σ _{ t } being 30% larger than σ _{ n }.
Boundary conditions
Boundary conditions need to be specified before Eq. 1 can be solved. It is common for boundaries of the electrode to be Dirchlet type and the remaining boundaries to be Neumann type [16]. Specifically, the boundary condition for a charged electrode will be
where V _{0} is the applied voltage of the electrode. For a grounded electrode, V _{0} would be zero and the boundary condition would read
The boundaries not in contact with an electrode, are often modeled as electrically insulating,
where n is the outward pointing normal [16].
Monopolar geometry
The geometry used for the monopolar simulations was the same as used by Castellvi et al. It consisted of two monopolar electrodes spaced 15 mm apart embedded in a 60 mm sphere. A schematic of the simulation geometry is shown in Fig. 2. The outer surface and electrode sleeves were modelled as electrically insulating. That is to say, on those boundaries the condition
where n is the outward pointing normal is enforced. The remaining boundary conditions is to set the active area of one electrode to ground by enforcing the condition
and to enforce the condition that the active area of the other electrode is held at a constant voltage by enforcing the equation
Bipolar geometry
The domain used for the bipolar simulations was a 60 mm sphere. The electrode has a diameter of 1.5 mm and a length of 40.5 mm. The last 7 mm of the electrode are set to the positive voltage. The next 7 mm of the electrode is insulated, and then followed by 7 mm of the electrode set to ground. The remaining portion of the electrode is set to an insulating boundary condition as is the outer sphere. A schematic of the domain is shown in Fig. 3.
Finite element solver
All simulations were performed using the commercial finite element software COMSOL. The meshes used for the simulations consisted of approximately 500,000 cells, and were run using 2.2 GHz Intel Xeon processors.
Results
Monopolar simulations
Both formulations were run with various voltages. A comparison of the the ablation volume can be found in Fig. 4. It can be seen that the anistropoic formulation predicts a slightly smaller ablation zone than the isotropic formulation. For all simulations, a volume was considered irreversibly electroporated if the magnitude of the electrical field was 184 V/cm or greater. This value was determined experimentally in [11] for the cases this work recreated for comparison.
The shape between the anisotropic and the isotropic formulation are similar except for the anisotropic formulation resulting in a slightly less spherical ablation zone than the isotropic formulation. Similar results were seen for other voltages.
Bipolar simulations
For the bipolar geometry, the volume predicted to be ablated is much lower for the anisotropic simulations than for the isotropic simulations. A summary of the results are displayed in Fig. 5. Notice that there is a larger percentage change between the isotropic and the anisotropic cases for the bipolar probe than for the monopolar probes. Ablation zones for the bipolar geometry are displayed in Figs. 5 for an applied voltage of 750 V. This was done to show that the anistotropicvarying formulation produces an ablation shape that is qualitatively similar, but quantitatively slightly smaller. Similar ablation shapes were seen for the other voltages as well.
Discussion and conclusions
This paper derived a general formulation for an anistropicvarying tensor for implementation into irreversible electroporation modeling software. The anistropicvarying tensor formulation allows the conductivity to take into consideration both electrical field direction and magnitude, as opposed to previous published works that only took into account electrical field magnitude. It was derived for arbitrary functions for σ _{ t } and σ _{ n }, and is therefore applicable irregardless of the form chosen for σ _{ t } and σ _{ n }. The formulation was compared to more commonly used isotropically varying formulations, and was found to decrease the predicted ablation zone for both monopolar and bipolar electrode setups.
Further experimental and numerical work is necessary before any definitive conclusions can be drawn on the importance of including an anisotropicvarying formulation in ablation zone predictions. The hope of this work is to encourage further research on how electrical fields can affect the conductivity tensor in different directions because more experimental data is necessary before any definitive conclusions can be drawn.
Abbreviations
 IRE:

Irreversible electroporation
 UQ:

Uncertainty quantification
References
 1
Dev SB, Dhar D, Krassowska W. Electric field of a sixneedle array electrode used in drug and dna delivery in vivo: analytical versus numerical solution. IEEE Trans Biomed Eng. 2003; 50(11):1296–300.
 2
Jourabchi N, Beroukhim K, Tafti BA, Kee ST, Lee EW. Irreversible electroporation (nanoknife) in cancer treatment. Gastrointest Interv. 2014; 3(1):8–18. doi:10.1016/j.gii.2014.02.002.
 3
Miklavčič D, Šemrov 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 Gen Subj. 2000; 1523(1):73–83.
 4
Labarbera N. Uncertainty quantification in irreversible electroporation simulations. Bioengineering. 2017; 4(2).
 5
Pavlin M, Miklavčič D. Theoretical and experimental analysis of conductivity, ion diffusion and molecular transport during cell electroporation—relation between shortlived and longlived pores. Bioelectrochemistry. 2008; 74(1):38–46.
 6
Ivorra A, AlSakere B, Rubinsky B, Mir LM. In vivo electrical conductivity measurements during and after tumor electroporation: conductivity changes reflect the treatment outcome. Phys Med Biol. 2009; 54(19):5949.
 7
Joshi R, Schoenbach K. Electroporation dynamics in biological cells subjected to ultrafast electrical pulses: a numerical simulation study. Phys Rev E. 2000; 62(1):1025.
 8
Corovic S, Lackovic I, Sustaric P, Sustar T, Rodic T, Miklavcic D. Modeling of electric field distribution in tissues during electroporation. Biomed Eng Online. 2013; 12(1):1.
 9
Mezeme ME, Kranjc M, Bajd F, Serša I, Brosseau C, Miklavčič D. Assessing how electroporation affects the effective conductivity tensor of biological tissues. Appl Phys Lett. 2012; 101(21):213702.
 10
Weaver JC, Chizmadzhev YA. Theory of electroporation: a review. Bioelectrochem Bioenerg. 1996; 41(2):135–60.
 11
Castellví Q, Banús J, Ivorra A. 3d assessment of irreversible electroporation treatments in vegetal models. In: 1st World Congress on Electroporation and Pulsed Electric Fields in Biology, Medicine and Food & Environmental Technologies. Singapore: Springer: 2016. p. 294–7.
 12
Golberg A, Rubinsky B. A statistical model for multidimensional irreversible electroporation cell death in tissue. Biomed Eng Online. 2010; 9(1):1.
 13
Laird AK. Dynamics of tumour growth. Br J Cancer. 1964; 18(3):490.
 14
Neal RE, et al. Simulation of in vivo irreversible electroporation renal ablations. In: 6th European Conference of the International Federation for Medical and Biological Engineering. Cham: Springer: 2015. p. 813–6.
 15
Banús Cobo J. 3d assessment of irreversible electroporation treatments in vegetal models. Barcelona: Thesis University Peompeu Fabra; 2015.
 16
Edd JF, Davalos RV. Mathematical modeling of irreversible electroporation for treatment planning. Technol Cancer Res Treat. 2008; 6(4):275–86.
Acknowledgements
The authors would like to acknowledge The Pennsylvania State University for providing computational resources.
Funding
The authors received no funding during this research.
Availability of data and materials
Data sets of the simulation data are available by contacting the corresponding author.
Author information
Affiliations
Contributions
NL derived the anistropic formulation and carried out the simulations. CD provided invaluable advice and help. Both authors read and approved the final manuscript.
Corresponding author
Correspondence to Nicholas Labarbera.
Ethics declarations
Authors’ information
NL and CD are both researchers at The Pennsylvania State University.
Ethics approval and consent to participate
Not applicable.
Consent for publication
The author’s give their consent to publish this work.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Irreversible electroporation
 Anisotropic conductivity
 Tumor ablation
 Treatment planning