Improving Prediction of the Potential Distribution Induced by Cylindrical Electrodes within a Homogeneous Rectangular Grid during Irreversible Electroporation

Featured Application: The Electric-Potential Estimation (EPE) method provides the possibility to considerably shorten the simulation time of treatment planning within a uniformly distributed rectangular grid by using large voxel-sized models without largely affecting the accuracy of the electric-ﬁeld distribution. Thus, real-time clinical IRE treatment planning in realistic heterogeneous target volumes becomes feasible. Abstract: Background: Irreversible electroporation (IRE) is an ablation technique based on the application of short, high-voltage pulses between needle electrodes (diameter: ~1.0 × 10 − 3 m). A Finite Difference-based software simulating IRE treatment generally uses rectangular grids, yielding discretization issues when modeling cylindrical electrodes and potentially affecting the validity of treatment planning simulations. Aim: Develop an Electric-Potential Estimation (EPE) method for accurate prediction of the electric-potential distribution in the vicinity of cylindrical electrodes. Methods: The electric-potential values in the voxels neighboring the cylindrical electrode voxels were corrected based on analytical solutions derived for coaxial/cylindrical electrodes. Simulations at varying grid resolutions were validated using analytical models. Low-resolution heterogeneous simulations at 2.0 × 10 − 3 m excluding/including EPE were compared with high-resolution results at 0.25 × 10 − 3 m. Results: EPE signiﬁcantly reduced maximal errors compared to analytical results for the electric-potential distributions (26.6–71.8% → 0.4%) and for the electrical resistance (30% → 1–6%) at 3.0 × 10 − 3 m voxel-size. EPE signiﬁcantly improved the mean-deviation (43.1–52.8% → 13.0–24.3%) and the calculation-time gain (>15,000 × ) of low-resolution compared to high-resolution heterogeneous simulations. Conclusions: EPE can accurately predict the potential distribution of neighboring cylindrical electrodes, regardless of size, position, and orientation in a rectangular grid. The simulation time of treatment planning can therefore be shortened by using large voxel-sized models without affecting accuracy of the electric-ﬁeld distribution, enabling real-time clinical IRE treatment planning.


Introduction
Irreversible Electroporation (IRE) is a local ablation modality that is currently investigated for the treatment of unresectable tumors, including liver, lung, pancreatic, and urological cancers [1][2][3][4][5]. Contrary to conventional thermal ablation techniques, in which ablation is achieved by exceeding the thermal-damage temperature threshold, the IRE procedure is based on the application of short, high-voltage pulses across needle or plate electrode pairs implanted in and surrounding the target volume to cause nanopores in the cell membranes [6]. The pore formation induces disturbances in the cell's homeostasis, with as a consequence cell death through accidental and regulated cell death mechanisms, and pore formation also improves solute mobility which translates in increased electrical conductivity in the bulk tissue [7,8].
IRE ablation depends on several parameters, such as the electrical pulse parameters (voltage, shape, duration, number, and frequency), temperature, and tissue properties of the target volume [9][10][11][12]. Therefore, it is important to optimize the applied IRE protocols to mainly obtain ablation through electroporation only or the combined effect of electroporation and mild hyperthermia, without achieving thermal damage [12]. In addition to pre-clinical in vitro and in vivo experiments, computational modeling provides a powerful tool for investigating and predicting the outcome of IRE-protocols. Using computational models, a wide variety of IRE protocol variations can be explored. Specifically, the influence of individual parameters can easily be investigated by varying a single parameter, while maintaining other parameters constant. This way, computational modelling can assist in the optimization of IRE protocols.
Currently, several commercial and non-commercial software packages are available to calculate distributions of physical key quantities that occur during an electroporation treatment, such as electric field and temperature, and the change in behavior of tissue properties as a function of the former quantities [12,13]. Examples of packages capable of implementing field and temperature dependent tissues properties include COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden), QuickField (Terra Analysis, Svendborg, Denmark), Eview (https://eview.upf.edu/; accessed on 4 May 2021) [14], IRENA (https: //team.inria.fr/monc/software/; accessed on 4 May 2021) [15], and OpenEP (https:// github.com/LSC-UBA/OpenEP/; accessed on 28 February 2021) [13]. Depending on the software package, the numerical simulations are generally performed using the Finite Element Method and/or the Finite Difference Method (FDM).
Some software packages allow for local adaptive mesh refinement in the vicinity of small or round structures to improve the modelling accuracy of needle electrodes or cylindrically shaped electrodes [16,17]. However, a uniform rectangular grid is usually used in FDM-based software packages [13,15,18,19]. The small diameter of the needle electrode applied in the clinic may result in reduced accuracy of the calculated electric potential in the vicinity of the needle electrode model in the case of an inappropriate choice of mesh properties in FDM-based software packages. Particularly, this is the case when large grid voxels (from here on the term "voxels" will be used instead of "mesh elements" as a reference for the elements of a uniform rectangular grid) are obtained from Computed Tomography (CT) scans without increasing the grid resolution. Specifically, the pixel size of a CT scan can vary from 0.5 × 10 −3 to 1.0 × 10 −3 m, while the slice thickness can vary from 2.0 × 10 −3 to 6.0 × 10 −3 m [20]. This means that a slight shift in the orientation or position of the electrode or the grid could result in different discretization of the needle electrode (see Figure 1), resulting in an inaccurate electric-potential distribution and, consequently, an inaccurate electric-field distribution. This can be solved by increasing the grid resolution. However, due to the uniformity of the rectangular grid, a decrease in the voxel size will be at the expense of the computer/GPU memories and speed. This would especially be the case when large models are used (i.e., additional memories would be required), and/or when model properties are updated in time due to the change in tissue properties as a function of the electric field and temperature (i.e., additional calculation time would be required). Appl. Sci. 2022, 12, x FOR PEER REVIEW 3 of 20 function of the electric field and temperature (i.e., additional calculation time would be required). In this study we present a method to increase the accuracy of the electric-field distribution within a low-resolution grid by estimating the electric-potential values in the nonmetallic 'estimation set' voxels neighboring the IRE-needle electrode voxels based on an analytical model. Electric-Potential Estimation (EPE) is a method that considers the circular surface of a needle electrode regardless of its size, position, and orientation in a uniformly distributed rectangular grid. In this study, we will determine and validate the effects of EPE on the electric-potential distribution and on the electrical resistance. This will be performed for needle shaped models such as coaxial cable (e.g., for modeling of radiofrequency ablation [21]) and cylindrical electrode pair (e.g., for modeling of electroporation [22,23]). In addition, we will present a case of rotated needle electrodes around a tumor within a low-resolution grid to illustrate the improved distribution of the electric field.

Summary of the Applied Finite Difference Method in Plan2Heat
In this study the software package Plan2Heat was used to simulate the electric-potential distribution and electrical resistance for IRE applications. Plan2Heat is an in-house developed treatment planning package for a variety of hyperthermia applications [24][25][26], also supporting IRE simulations. Since the energy density of the electric field resulting from the IRE pulses is much larger than the energy density of the magnetic field, the electro-quasi-static approximation can be applied [12]. The quasi-static solver of Plan2Heat calculates the electric-potential distribution using the Maxwell's equations by combining Faraday's law of induction, and Ampère's law In this study we present a method to increase the accuracy of the electric-field distribution within a low-resolution grid by estimating the electric-potential values in the non-metallic 'estimation set' voxels neighboring the IRE-needle electrode voxels based on an analytical model. Electric-Potential Estimation (EPE) is a method that considers the circular surface of a needle electrode regardless of its size, position, and orientation in a uniformly distributed rectangular grid. In this study, we will determine and validate the effects of EPE on the electric-potential distribution and on the electrical resistance. This will be performed for needle shaped models such as coaxial cable (e.g., for modeling of radio-frequency ablation [21]) and cylindrical electrode pair (e.g., for modeling of electroporation [22,23]). In addition, we will present a case of rotated needle electrodes around a tumor within a low-resolution grid to illustrate the improved distribution of the electric field. In this study the software package Plan2Heat was used to simulate the electricpotential distribution and electrical resistance for IRE applications. Plan2Heat is an in-house developed treatment planning package for a variety of hyperthermia applications [24][25][26], also supporting IRE simulations. Since the energy density of the electric field resulting from the IRE pulses is much larger than the energy density of the magnetic field, the electro-quasi-static approximation can be applied [12]. The quasi-static solver of Plan2Heat calculates the electric-potential distribution using the Maxwell's equations by combining Faraday's law of induction, and Ampère's law where:

Materials and Methods
is the angular frequency; • µ r is the dimensionless relative permeability; is the permittivity of free space.
Here, the vector quantities are expressed in bold and italic. Assuming that the magnetic field is negligible, we can transform Equation (1) into: allowing the electric-field component to be expressed as: with Φ (V) as in the scalar electric potential. Using the electro-quasi-static assumption, the continuity equation: can be turned into: where J (A·m −2 ) is the electrical current density, and κ = σ + jωε r ε 0 (S·m −1 ) is the complex admittance.
To calculate the electric-potential distribution, Equation (6) was converted to a parabolic partial differential equation by introducing a pseudo-time τ (s), and a spatial-dependent relaxation function F(x, y, z) to smooth large differences in admittance, such that Equation (6) was turned into: 1 F(x, y, z) with (x, y, z) as Cartesian coordinates expressed in (m) [19]. Subsequently, Equation (7) was solved by marching on in pseudo-time, allowing the potential Φ to diffuse in space starting from an initial distribution, until a stationary state (∂Φ/∂τ = 0 V·s −1 ) was reached.

Discretization and Gridding in Plan2Heat
Plan2Heat performs calculations using the FDM within a model that is discretized on a uniform rectangular (tissue) grid. Examples of such grids are shown in Figure 1. Each voxel is assigned appropriate dielectric tissue properties for calculations. Regarding the discretization shown in Figure 1, if an electrode covers the voxel node (black dot; the center of the voxel), then the voxel is modelled as an electrode voxel (grey voxel).
Furthermore, in this study we assumed two types of grids: • Centralized grid ( Figure 1A): the grid is chosen such that the center of a voxel coincides with the center of the electrode. In this example, the electrode is discretized by five voxels; • Shifted grid ( Figure 1B): a grid shift of 1 2 equilateral voxel size in both x-and ydirections was realized by reducing the grid size by one voxel in both x-and ydirections and centralizing the grid again.

Electric-Potential Estimation: Correction Using Estimation Voxels
In this study an Electric-Potential Estimation method is used to increase the accuracy of the simulated electric-field distribution close to the needle surface independent of the grid resolution. This is achieved by estimating and correcting electric-potential values near the IRE-needle by using a so-called 'estimation set' of non-metallic voxels neighboring the IREneedle combined with an analytical model that considers the exact location of the electrode surface within the grid. This means that the position and effect of the circular surface of a needle electrode are correctly represented regardless of its size, position, and orientation in any uniformly distributed rectangular grid by adjusting the electric-potential values in the non-metallic voxels neighboring the electrode voxels (the estimation voxels). Examples of estimation voxels for different discretized electrodes are illustrated as orange voxels in Figure 2. The Electric-Potential Estimation (EPE) values were assessed by assigning spatially-depended values derived from analytical models [27,28]. The incorporation of such models in EPE is explained in the next subsections.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 20 • Shifted grid ( Figure 1B): a grid shift of ½ equilateral voxel size in both x-and y-directions was realized by reducing the grid size by one voxel in both x-and y-directions and centralizing the grid again.

Electric-Potential Estimation: Correction Using Estimation Voxels
In this study an Electric-Potential Estimation method is used to increase the accuracy of the simulated electric-field distribution close to the needle surface independent of the grid resolution. This is achieved by estimating and correcting electric-potential values near the IRE-needle by using a so-called 'estimation set' of non-metallic voxels neighboring the IRE-needle combined with an analytical model that considers the exact location of the electrode surface within the grid. This means that the position and effect of the circular surface of a needle electrode are correctly represented regardless of its size, position, and orientation in any uniformly distributed rectangular grid by adjusting the electric-potential values in the non-metallic voxels neighboring the electrode voxels (the estimation voxels). Examples of estimation voxels for different discretized electrodes are illustrated as orange voxels in Figure 2. The Electric-Potential Estimation (EPE) values were assessed by assigning spatially-depended values derived from analytical models [27,28]. The incorporation of such models in EPE is explained in the next subsections.

Coaxial Cable Model for a Single IRE Electrode in Tissue Cylinder
For simplicity we start with the description of an analytical coaxial cable model with a cylindrical inner electrode radius R1 (m), and a cylindrical outer electrode radius R2 (m) with negligible thickness (see Figure 3). The length of these electrodes were assumed to be much larger than R and their centers were assumed to be located in the origin. The electric-potential distribution within this model can be analytically expressed as:

Coaxial Cable Model for a Single IRE Electrode in Tissue Cylinder
For simplicity we start with the description of an analytical coaxial cable model with a cylindrical inner electrode radius R 1 (m), and a cylindrical outer electrode radius R 2 (m) with negligible thickness (see Figure 3). The length of these electrodes were assumed to be much larger than R 2 and their centers were assumed to be located in the origin. The electric-potential distribution within this model can be analytically expressed as: Φ(x, y) = I 2πLκ ln R 1 where Φ(x, y) is the electric-potential value of a voxel node at position (x, y) with respect to the origin, Φ BEM|1 (V) is the electric potential of the inner electrode, L (m) as the modeled length of the electrodes, and I (A) is the electric current that is calculated by: where Φ BEM|2 (V) is the electric potential of the outer electrode. The electrical resistance was calculated using: with the conduction shape factor S (m): (11) length of the electrodes, and I (A) is the electric current that is calculated by: where ΦBEM|2 (V) is the electric potential of the outer electrode. The electrical resistance was calculated using: with the conduction shape factor S (m): To numerically model the coaxial cable, the length of the electrodes was virtually increased to infinity in the z-direction by isolating the top and the bottom surfaces of the model (κ = 0 S·m −1 ; illustrated by the color aqua blue in Figure 3B), allowing the modeled electrode length L to have the size of a single voxel. More details about the domain properties are provided in Section 2.3. Analytical Validation. Next, the EPE was incorporated by selecting the x,y-coordinates of the estimation voxels for both the inner and the outer electrodes ( Figure 3C). Subsequently, these x,y-coordinates were used as arguments in Equation (8) to determine the EPE values of the estimation voxels, after which the EPE values were assigned to the estimation voxels as a Dirichlet boundary condition before calculating the electric-potential distribution of the model. (C) Discretization of the inner and outer electrodes of the coaxial cables (represented by the circular lines). Electrode voxels are grey and their corresponding estimation voxels are represented by the orange voxels. The red contour represents the surfaces of the voxels that were used to calculate the electrical resistance. (C) Discretization of the inner and outer electrodes of the coaxial cables (represented by the circular lines). Electrode voxels are grey and their corresponding estimation voxels are represented by the orange voxels. The red contour represents the surfaces of the voxels that were used to calculate the electrical resistance.
To numerically model the coaxial cable, the length of the electrodes was virtually increased to infinity in the z-direction by isolating the top and the bottom surfaces of the model (κ = 0 S·m −1 ; illustrated by the color aqua blue in Figure 3B), allowing the modeled electrode length L to have the size of a single voxel. More details about the domain properties are provided in Section 2.3. Analytical Validation. Next, the EPE was incorporated by selecting the x,y-coordinates of the estimation voxels for both the inner and the outer electrodes ( Figure 3C). Subsequently, these x,y-coordinates were used as arguments in Equation (8) to determine the EPE values of the estimation voxels, after which the EPE values were assigned to the estimation voxels as a Dirichlet boundary condition before calculating the electric-potential distribution of the model.

Cylindrical IRE Electrode Pair Model
For the needle electrode model a pair of cylindrical electrodes were assumed with radius R (m) and distance d (m) between the electrodes axis that are parallel to the z-axis. The electrode pair were positioned at (+d/2, 0, 0) and (−d/2, 0, 0) with respect to the origin Appl. Sci. 2022, 12, 1471 7 of 20 (see Figure 4). Furthermore, Φ(x, y) was calculated with respect to the origin assuming that the length of the electrodes are much larger than x 2 + y 2 and d/2. Based on the analytical electric-potential distribution derived by Assis and Mania [28], the electric potential of the estimation voxels were estimated using: with Φ BEM|1 (V) as the electric potential of the left electrode and Φ BEM|2 (V) as the electric potential of the right electrode. The electrical resistance was calculated using the Equation (10) and the conduction shape factor [29]: tential of the estimation voxels were estimated using: with ΦBEM|1 (V) as the electric potential of the left electrode and ΦBEM|2 (V) as the electric potential of the right electrode. The electrical resistance was calculated using the Equation (10) and the conduction shape factor [29]: To numerically model the cylindrical electrode pair, the length of the electrodes was virtually increased to infinity in the z-direction by isolating the top and the bottom surfaces of the model (κ = 0 S·m −1 ; illustrated by the color aqua blue in Figure 4B), allowing the electrode length L to have a size of a single voxel. Again, more details about the domain properties are provided in Section 2.3. Analytical Validation. As in Section 2.2.1. Coaxial Cable Model for a Single IRE Electrode in Tissue Cylinder the EPE was incorporated by selecting the x,y-coordinates of the estimation voxels for both electrodes ( Figure 4C). Subsequently, these x,y-coordinates were used as arguments in Equation (12) to determine the EPE values of the estimation voxels, after which the EPE values were assigned to the estimation voxels as a Dirichlet boundary condition before calculating the electricpotential distribution of the model.  To numerically model the cylindrical electrode pair, the length of the electrodes was virtually increased to infinity in the z-direction by isolating the top and the bottom surfaces of the model (κ = 0 S·m −1 ; illustrated by the color aqua blue in Figure 4B), allowing the electrode length L to have a size of a single voxel. Again, more details about the domain properties are provided in Section 2.3. Analytical Validation. As in Section 2.2.1. Coaxial Cable Model for a Single IRE Electrode in Tissue Cylinder the EPE was incorporated by selecting the x,y-coordinates of the estimation voxels for both electrodes ( Figure 4C). Subsequently, these x,y-coordinates were used as arguments in Equation (12) to determine the EPE values of the estimation voxels, after which the EPE values were assigned to the estimation voxels as a Dirichlet boundary condition before calculating the electric-potential distribution of the model.

Analytical Validation
The analytical solutions for the electric-potential distribution and electrical resistances were compared to the numerical solutions with and without using EPE for both centralized and shifted grids (see Figure 1). The comparison was performed for equilateral voxel dimensions that varied between 0.5 × 10 −3 m and 3.0 × 10 −3 m in the (x, y)-plane, which are comparable with the pixel size range (0.5 × 10 −3 -1.0 × 10 −3 m) and slice thickness size range (2.0 × 10 −3 -6.0 × 10 −3 m) used for CT scans [20]. Choosing 3.0 × 10 −3 m as the maximal voxel dimension yields sufficient voxels between the electrodes to provide an adequate overview of the electric-field distribution. Furthermore, the Dirichlet boundary condition (Φ = 0 V) was applied to the remaining BOS of the model.
As mentioned before, the electric-potential values of the estimation voxels are fixed. To calculate the electrical resistance considering EPE, the total electrical current was calculated using: with n as a unity vector normal to the red contour shown in Figures 3C and 4C. In a discrete manner [19,27], the electric current density between two adjacent voxel nodes i and j can be calculated by: with: where: • i is the node (the center) of an estimation voxel; • j is the voxel node neighboring the voxel node i across the red contour; • I ij (A) is the current flow between voxel nodes i and j; is the current density between voxel nodes i and j; • A ij (m 2 ) is the area of the surface between voxel nodes i and j; is the admittance of the voxel with voxel node i, it is assumed that each voxel has a uniform admittance; • κ j (S·m −1 ) is the admittance of the voxel with voxel node j; • κ ij (S·m −1 ) is the effective admittance between the voxel nodes i and j; • h ij (m) is the distance between the voxel nodes i and j; is the electric-potential value at the voxel node i.
Then, along the red contour (or for all estimation voxels) I total is calculated using: The analytical models described in the previous section were used to validate EPE, since their cylindrical shape is similar to the shape of IRE needle electrodes applied in the clinic. In the next paragraph, the choices for the model parameters are presumed based on parameters used in the clinic: For the coaxial cable model, the inner electrode was assumed to have a diameter of 1 × 10 −3 m, which is comparable with the clinical IRE needle electrodes [30]. The outer electrode diameter was assumed to be 40 × 10 −3 m to attain a distance of 20 × 10 −3 m between the inner and outer electrode. The voltage across the electrodes was presumed to be 1500 V. Similarly, the electrode diameters of the cylindrical electrodes were assumed to be 1 × 10 −3 m with an interdistance of 20 × 10 −3 m. The voltage across the electrodes was presumed to be 3000 V. For both models, healthy pancreatic tissue was assumed with κ = 0.511 S·m −1 [31], resulting in analytical electrical resistances of 2297.9 Ω and 4594.9 Ω for the coaxial cable model and the cylindrical electrodes model, respectively. Please note that in both models, L was assumed to have the length of the voxel dimension in the z-direction, which was 0.5 × 10 −3 m.

Application in Heterogeneous Tumor Models
The next step is to assess the performance of EPE in realistic tumor models with heterogeneous tissue properties and different electrode orientations, as would occur in Appl. Sci. 2022, 12, 1471 9 of 20 clinical use. An example of EPE application to a 3D model was presented to demonstrate the performance of EPE in treatment planning. A spherical pancreatic tumor illustrated in Figure 5 was modeled with a pair of cylindrical electrodes at its extreme ends. The diameter of the tumor and the electrodes lengths were assumed to be 20 × 10 −3 m, and the interdistance between the electrodes was assumed to be 19 × 10 −3 m with the presumed voltage of 2850 V across the electrodes (voltage-to-distance ratio of 1500 V·cm −1 ) [30]. The electric-field threshold of the pancreatic cancer was assumed to be 500 V·cm −1 , based on pancreatic cell line experiments [32]. We also changed the electrode orientation: the electrode-tumor model was rotated around the y-axis by 45 • (see Figure 5A) and placed in the center of a cube representing healthy pancreatic tissue, with dimensions that vary between 90 × 10 −3 and 94 × 10 −3 m depending on voxel size. Furthermore, a centralized grid was applied to the cubic model with an equilateral voxel dimension of 2 × 10 −3 m (diameter of a clinical needle electrode is~1 × 10 −3 m) to simulate a worst-case scenario regarding electrode diameter, position, and orientation. Next, two heterogeneous model cases were presented: • of 2297.9 Ω and 4594.9 Ω for the coaxial cable model and the cylindrical electrodes model, respectively. Please note that in both models, L was assumed to have the length of the voxel dimension in the z-direction, which was 0.5 × 10 −3 m.

Application in Heterogeneous Tumor Models
The next step is to assess the performance of EPE in realistic tumor models with heterogeneous tissue properties and different electrode orientations, as would occur in clinical use. An example of EPE application to a 3D model was presented to demonstrate the performance of EPE in treatment planning. A spherical pancreatic tumor illustrated in Figure 5 was modeled with a pair of cylindrical electrodes at its extreme ends. The diameter of the tumor and the electrodes lengths were assumed to be 20 × 10 −3 m, and the interdistance between the electrodes was assumed to be 19 × 10 −3 m with the presumed voltage of 2850 V across the electrodes (voltage-to-distance ratio of 1500 V·cm −1 ) [30]. The electric-field threshold of the pancreatic cancer was assumed to be 500 V·cm −1 , based on pancreatic cell line experiments [32]. We also changed the electrode orientation: the electrodetumor model was rotated around the y-axis by 45° (see Figure 5A) and placed in the center of a cube representing healthy pancreatic tissue, with dimensions that vary between 90 × 10 −3 and 94 × 10 −3 m depending on voxel size. Furthermore, a centralized grid was applied to the cubic model with an equilateral voxel dimension of 2 × 10 −3 m (diameter of a clinical needle electrode is ~1 × 10 −3 m) to simulate a worst-case scenario regarding electrode diameter, position, and orientation. Next, two heterogeneous model cases were presented: • Case 1: the electrical conductivities of the human cancerous (sphere) and healthy (cube) tissues were assumed to be 0.20 and 0.35 S·m −1 , respectively [32]; • Case 2: the assumed conductivity values of these tissue types were interchanged. Along the length of the electrodes the electric-potential values of EPE were assumed to be constant. Furthermore, a Dirichlet boundary condition (Φ = 0 V) was applied to all BOS of the model. The effect of EPE was then evaluated by determining the electric-field distributions after the exclusion and inclusion of EPE, and compare the results to the simulation results of the same model along the lines O1 and O2 (see Figure 5) with (1) highresolution (an equilateral voxel dimension of 0.25 × 10 −3 m), and (2) with a model in which the electrode axes were positioned on the x-axis and was parallel to the z-axis Figure 5B. Along the length of the electrodes the electric-potential values of EPE were assumed to be constant. Furthermore, a Dirichlet boundary condition (Φ = 0 V) was applied to all BOS of the model. The effect of EPE was then evaluated by determining the electric-field distributions after the exclusion and inclusion of EPE, and compare the results to the simulation results of the same model along the lines O 1 and O 2 (see Figure 5) with (1) highresolution (an equilateral voxel dimension of 0.25 × 10 −3 m), and (2) with a model in which the electrode axes were positioned on the x-axis and was parallel to the z-axis Figure 5B. Outcomes were expressed as a mean of the deviations of the electric-field distributions of the low-resolution models ex/including EPEs relative to the electric-field values at the corresponding position within the high-resolution model within a radius of 15 × 10 −3 m from the center. The radius 15 × 10 −3 m was selected to include the whole tumor and a part of the surrounding healthy tissue, which mimics the clinical ablation zone that aims to include the tumor-free margin [33]. The electric field was calculated using the method described in Appendix A.

Validation of the Single Electrode Coaxial Cable Model
In Figure 6 the analytical and numerical solutions of electric potential within a central grid are presented for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m. A comparison between both solutions was illustrated as a relative error with respect to the analytical solution in Figure 7. The datapoints of the numerical solutions lined up with the analytical solution after the inclusion of Electric-Potential Estimation (EPE) for all voxel sizes. Specifically, the inclusion of EPE significantly reduced the maximal errors from 87.8%, 88.0%, and 47.4% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.1%, 1.3%, and 0.4%, respectively. Note that these maximal errors without EPE occur in the far field, but the maximal errors without EPE near the inner electrode are also very high, 7.5%, 3.3%, and 26.6% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, respectively, see Figure 7. These errors near the electrode also reduced to values smaller than 0.4% after inclusion of EPE.
Furthermore, in Figure 8 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. Inclusion of EPE significantly reduced the maximal errors from 7.2%, 2.4%, and 30.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.2%, 1.6%, and 2.8%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 31.3% for the voxel sizes 0.5 × 10 −3 m and 3.0 × 10 −3 m, to 0.1% and 6.2%, respectively, except for the voxel size 1.0 × 10 −3 m. While the error remains low, it slightly increased from 2.4% to 3.5%, possibly due to the fact that the voxel size has the same width as the diameter of the cylindrical electrode. Outcomes were expressed as a mean of the deviations of the electric-field distributions of the low-resolution models ex/including EPEs relative to the electric-field values at the corresponding position within the high-resolution model within a radius of 15 × 10 −3 m from the center. The radius 15 × 10 −3 m was selected to include the whole tumor and a part of the surrounding healthy tissue, which mimics the clinical ablation zone that aims to include the tumor-free margin [33]. The electric field was calculated using the method described in Appendix A.

Validation of the Single Electrode Coaxial Cable Model
In Figure 6 the analytical and numerical solutions of electric potential within a central grid are presented for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m. A comparison between both solutions was illustrated as a relative error with respect to the analytical solution in Figure 7. The datapoints of the numerical solutions lined up with the analytical solution after the inclusion of Electric-Potential Estimation (EPE) for all voxel sizes. Specifically, the inclusion of EPE significantly reduced the maximal errors from 87.8%, 88.0%, and 47.4% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.1%, 1.3%, and 0.4%, respectively. Note that these maximal errors without EPE occur in the far field, but the maximal errors without EPE near the inner electrode are also very high, 7.5%, 3.3%, and 26.6% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, respectively, see Figure 7. These errors near the electrode also reduced to values smaller than 0.4% after inclusion of EPE.   Furthermore, in Figure 8 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. Inclusion of EPE significantly reduced the maximal errors from 7.2%, 2.4%, and 30.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.2%, 1.6%, and 2.8%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 31.3% for the voxel sizes 0.5 × 10 −3 m and 3.0 × 10 −3 m, to 0.1% and 6.2%, respectively, except for the voxel size 1.0 × 10 −3 m. While the error remains low, it slightly increased from 2.4% to 3.5%, possibly due to the fact that the voxel size has the same width as the diameter of the cylindrical electrode.

Validation of the Cylindrical Electrode Pair
In Figure    Furthermore, in Figure 8 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. Inclusion of EPE significantly reduced the maximal errors from 7.2%, 2.4%, and 30.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.2%, 1.6%, and 2.8%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 31.3% for the voxel sizes 0.5 × 10 −3 m and 3.0 × 10 −3 m, to 0.1% and 6.2%, respectively, except for the voxel size 1.0 × 10 −3 m. While the error remains low, it slightly increased from 2.4% to 3.5%, possibly due to the fact that the voxel size has the same width as the diameter of the cylindrical electrode.

Validation of the Cylindrical Electrode Pair
In Figure 9 the analytical and numerical solutions of electric potential within a central grid are presented for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m. A comparison between both solutions was illustrated as a relative error with respect to the analytical solution in Figure 10. According to the results, the datapoints of the numerical solutions again lined up with the analytical solution after the inclusion of EPE for all voxel sizes.

Validation of the Cylindrical Electrode Pair
In Figure 9 the analytical and numerical solutions of electric potential within a central grid are presented for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m. A comparison between both solutions was illustrated as a relative error with respect to the analytical solution in Figure 10. According to the results, the datapoints of the numerical solutions again lined up with the analytical solution after the inclusion of EPE for all voxel sizes. Specifically, the inclusion of EPE significantly reduced the maximal errors from 9.3%, 4.0%, and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.5%, 1.7%, and 0%, respectively, in the target volume between the cylindrical electrodes. Without EPE, the errors near the electrode were 7.9%, 2.7% and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, respectively, which also significantly reduced to values equal to or smaller than 0.4% after inclusion of EPE.
12, x FOR PEER REVIEW 12 of 20 Specifically, the inclusion of EPE significantly reduced the maximal errors from 9.3%, 4.0%, and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.5%, 1.7%, and 0%, respectively, in the target volume between the cylindrical electrodes. Without EPE, the errors near the electrode were 7.9%, 2.7% and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, respectively, which also significantly reduced to values equal to or smaller than 0.4% after inclusion of EPE.  Furthermore, in Figure 11 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. EPE significantly reduced the maximal errors in the region in between the electrodes from 7.1%, 3.2%, and 29.6% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, to 0.5%, 0.8% and 1.0%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 25.5% for the voxel sizes 0.5 Specifically, the inclusion of EPE significantly reduced the maximal errors from 9.3%, 4.0%, and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m, and 3.0 × 10 −3 m, to 0.5%, 1.7%, and 0%, respectively, in the target volume between the cylindrical electrodes. Without EPE, the errors near the electrode were 7.9%, 2.7% and 71.8% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, respectively, which also significantly reduced to values equal to or smaller than 0.4% after inclusion of EPE.  Furthermore, in Figure 11 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. EPE significantly reduced the maximal errors in the region in between the electrodes from 7.1%, 3.2%, and 29.6% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, to 0.5%, 0.8% and 1.0%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 25.5% for the voxel sizes 0.5 Furthermore, in Figure 11 the relative error of the numerically calculated resistance values with respect to the analytical value were shown as a function of voxel sizes and grid types. EPE significantly reduced the maximal errors in the region in between the electrodes from 7.1%, 3.2%, and 29.6% for the voxel sizes 0.5 × 10 −3 m, 1.0 × 10 −3 m and 3.0 × 10 −3 m, to 0.5%, 0.8% and 1.0%, respectively, for the central grid. For the shifted grid, EPE also significantly reduced the maximal errors from 2.3% and 25.5% for the voxel sizes 0.5 × 10 −3 m and 3.0 × 10 −3 m, to 0.6% and 1.3%, respectively. For the voxel size 1.0 × 10 −3 m the error slightly increased from 1.8% to 2.5%.

Application in Heterogeneous Tumor Models
In the next section the results are provided for the example cases described in Section 2.4. Application in Heterogeneous Tumor Models. Additional simulation data regarding the electrical resistance values and the calculation time are presented in Table S1. 3.3.1. Case 1 Figure 12A illustrates an example of the electric-field distribution of the high-resolution model provided in Figure 5B, showing that the tumor is fully covered by an electric field larger than the electric-field threshold (~500 V·cm −1 ). Note that this example was intended to highlight the impact of a well-defined form of heterogeneity, rather than to demonstrate how a therapeutically optimal field distribution can be achieved. The electric-field distributions of the high-and low-resolution models along O1 and O2 lines are presented in Figure 12B, with the associated relative mean deviation shown in Table 1. These results demonstrate a significant reduction in the relative mean deviations after the inclusion of EPE by 30.7% along O1 and 25.6% along O2. In addition, the inclusion of EPE significantly improved the electrical resistance by reducing the deviation relative to the high-resolution model from 28.4% to 18.6%. Furthermore, the simulation of the high-resolution model required a calculation time of 22,156 s (369 min), while the simulation of the low-resolution models only required 1.4 s (a gain factor of ~16,000).

Application in Heterogeneous Tumor Models
In the next section the results are provided for the example cases described in Section 2.4. Application in Heterogeneous Tumor Models. Additional simulation data regarding the electrical resistance values and the calculation time are presented in Table S1. 3.3.1. Case 1 Figure 12A illustrates an example of the electric-field distribution of the high-resolution model provided in Figure 5B, showing that the tumor is fully covered by an electric field larger than the electric-field threshold (~500 V·cm −1 ). Note that this example was intended to highlight the impact of a well-defined form of heterogeneity, rather than to demonstrate how a therapeutically optimal field distribution can be achieved. The electric-field distributions of the high-and low-resolution models along O 1 and O 2 lines are presented in Figure 12B, with the associated relative mean deviation shown in Table 1. These results demonstrate a significant reduction in the relative mean deviations after the inclusion of EPE by 30.7% along O 1 and 25.6% along O 2 . In addition, the inclusion of EPE significantly improved the electrical resistance by reducing the deviation relative to the high-resolution model from 28.4% to 18.6%. Furthermore, the simulation of the highresolution model required a calculation time of 22,156 s (369 min), while the simulation of the low-resolution models only required 1.4 s (a gain factor of~16,000). Appl

Case 2
The electric-field distributions of the high-and low-resolution models along O1 and O2 lines were presented in Figure 13, with the associated relative mean deviation shown in Table 1. Similar to the results of Case 1, the relative mean deviations demonstrate a significant reduction after the inclusion of EPE by 30.1% along O1 and 30.0% along O2. In addition, the inclusion of EPE significantly improved the electrical resistance by reducing the deviation relative to the high-resolution model from 29.9% to 16.5%. Furthermore, the simulation of the high-resolution model also required a long calculation time of 19,983 s (333 min), while the simulation of the low-resolution models only required 1.3 s (a gain factor of ~15,000).

Discussion
An Electric-Potential Estimation (EPE) method for accurate prediction of the potential distribution in the vicinity of cylindrical electrodes was developed. This study demonstrated that the use of EPE can significantly improve the accuracy of both the electricpotential distribution and the electrical resistance for any grid resolution and electrodes configuration (diameter, orientation, and position). The improvement is particularly large at low-resolution, which can be attributed to the fact that EPE creates a grid-independent rendition of the electrode in the grid. Specifically, the EPE considers the circular surface of the applied electrodes by reshaping the electric-potential distributions in the direct vicinity and between the electrodes, resulting in a more accurate prediction of the electricfield distribution, particularly, within low-resolution models with a voxel size larger than the electrode diameter. The improvement of the accuracy means the low-resolution model

Case 2
The electric-field distributions of the high-and low-resolution models along O 1 and O 2 lines were presented in Figure 13, with the associated relative mean deviation shown in Table 1. Similar to the results of Case 1, the relative mean deviations demonstrate a significant reduction after the inclusion of EPE by 30.1% along O 1 and 30.0% along O 2 . In addition, the inclusion of EPE significantly improved the electrical resistance by reducing the deviation relative to the high-resolution model from 29.9% to 16.5%. Furthermore, the simulation of the high-resolution model also required a long calculation time of 19,983 s (333 min), while the simulation of the low-resolution models only required 1.3 s (a gain factor of~15,000).

Discussion
An Electric-Potential Estimation (EPE) method for accurate prediction of the potential distribution in the vicinity of cylindrical electrodes was developed. This study demonstrated that the use of EPE can significantly improve the accuracy of both the electricpotential distribution and the electrical resistance for any grid resolution and electrodes configuration (diameter, orientation, and position). The improvement is particularly large at low-resolution, which can be attributed to the fact that EPE creates a grid-independent rendition of the electrode in the grid. Specifically, the EPE considers the circular surface of the applied electrodes by reshaping the electric-potential distributions in the direct vicinity and between the electrodes, resulting in a more accurate prediction of the electric-field distribution, particularly, within low-resolution models with a voxel size larger than the electrode diameter. The improvement of the accuracy means the low-resolution model becomes reliable, which is associated with a significant gain in the calculation time of, e.g., 15,000× as was presented in Section 3.3. Application in Heterogeneous Tumor Models, which makes it attractive to be implemented into the numerical treatment planning of electrical energy-based ablation modalities, such as IRE.
Next to IRE, the EPE can also be utilized in numerical calculations of other electrical energy-based ablation modalities such as reversible electroporation, and thermal ablation, where no local adaptive mesh refinement is utilized. For example, Zorbas and Samaras investigated the effects of realistic geometry and blood vessels on the treatment outcome of radiofrequency ablation using the Finite Difference Method by modeling a single cylindrical electrode with a diameter of 1.5 × 10 −3 m within an uniform grid (voxel size of 0.5 × 10 −3 m in the plane normal to the electrode) [21,34]. As illustrated in Figure 1, a slight shift in the electrode position can alter the electrode discretization resulting in inaccurate electric-field distribution. In such cases, the EPE can be applied to improve the accuracy of the electric-potential distribution in the vicinity of the electrode. In this example, the coaxial cable model is recommended to validate the modeling of a single electrode, though one may also use different analytical models. For example, in a study performed by Gallinato et al. [15], in which numerical modeling of electroporation was provided in the liver, the authors designated a fixed electric-potential amplitude of Φ 1 /5 (V) to the estimation voxels to avoid a large increase in electric conductivities in the vicinity of the needle electrodes due to the overestimation of the electric field. For a low-resolution grid, a fixed electric-potential can distort the electric-field distribution in between the needle pairs, and therefore, one must be cautious with the selection of the appropriate electric-potential values near the electrode.
In this study the presented homogeneous and heterogeneous examples, e.g., in Figures 4 and 5, are less complex than the real anatomies, since the real anatomies display far more spatial heterogeneity in tissue properties (different tissue types). Still, reconstruction of these real anatomies is based on 1-2 mm resolution CT or MR scans, thus we can expect that our low-resolution simulations will yield reliable results in realistic clinical heterogeneous tumor models. In addition to the simplification of the tissues, the electrical conductivities in this study were assumed to be constant, as we focused on evaluating the mechanism of EPE. The changing characteristics of the electrical conductivity during an electroporation treatment are often taken into consideration in the numerical calculations by spatially updating the electrical conductivity, mainly, as a function of electric field solely, or electric field and temperature [8,12,22,[35][36][37]. Thus, in future studies, further investigations should be conducted into the incorporation of the heterogeneous dynamic spatial changes in the electrical conductivity during electroporation in EPE, and investigate the effect on the electric field distributions.
Note that in Section 2.4. Application in Heterogeneous Tumor Models, the electricpotential values of EPE were assumed to be constant along the longitudinal axis of the electrodes. In addition, these values were only varied along transversal planes of the electrodes. However, the electric-potential distributions would vary along the longitudinal axis of the electrode when a relatively short electrode is utilized with a length in the order of the distance between the electrodes, [23,38]. In future studies, this can be addressed, for example, by multiplying the EPE values with a function that estimates the electricpotential distribution along the longitudinal axis of the electrodes. Further research should be performed to determine a reliable estimation.
IRE is a relatively novel ablation technique, and most of the simulation studies published use idealized conditions which assume homogeneous tissue properties and perfectly parallel alignment of the electrodes, which are not representative of realistic clinical conditions [12]. Some validation studies have been performed to achieve robust treatment planning suitable for introduction in the clinical practice [39,40]. Considering various quantities and parameters (e.g., electrode configuration and electrical conductivity dependence on electric field and temperature) during such validation processes, the EPE can significantly reduce the calculation time, allowing for the inclusion of additional parameters for further improvement of the accuracy. In addition, due to the reduction in calculation time, the inclusion of EPE makes it possible to perform personalized IRE treatment planning as it is possible to evaluate several different treatment plan options using the most recently registered patient images shortly prior to IRE treatment delivery, or even after the electrode implantation is performed. This increases the likelihood of an optimal treatment delivery and more accurate prediction of the treated region and measured quantities, such as electrical current and resistance, which can be compared with the actual values during treatment delivery. Fast EPE-based low-resolution planning can also be clinically applied for a more reliable assessment of the size of the possible thermal ablation region in the high field region close to the electrode [12], for assessing the IRE effect close to specific structures in the treatment region (e.g., bile ducts, blood vessels, etc.) [41], and to assess the most favorable needle orientation in the presence of stents or other small metallic structures that may present a risk on local thermal damage [5,42]. Further research should be performed to determine the reliability of the assessment of fast low-resolution EPE in these clinical applications.

Conclusions
This study presented the Electric-Potential Estimation (EPE) method, an approach that can accurately predict the potential distribution in the vicinity of cylindrical electrodes, regardless of electrode size, position, and orientation in a uniformly distributed rectangular grid. High grid size-independent accuracy is achieved by more correctly estimating the electric-potential values of the voxels neighboring on the electrode voxels based on the actual physical distance to the electrode. This EPE method provides a possibility to considerably shorten the simulation time of treatment planning by using large voxel-sized models without largely affecting the accuracy of the electric-field distribution. Thus, real-time clinical IRE treatment planning in realistic heterogeneous target volumes becomes feasible.

Conflicts of Interest:
Authors declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: K. P. van Lienden and T. M. van Gulik are paid consultants for Angio-Dynamics. Angio-Dynamics had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
Assume an electric-potential distribution within an (x, y)-plane as shown in Figure A1. Within the grid, every voxel is indicated by the letters (m, n) that represent the location of the voxel inside (x, y)-plane; m and n represent the location of the voxel along the x-and y-axis, respectively. For simplicity, the main focus is to calculate the electric-field value within the voxel (m, n). Next, considering the adjacent voxels to voxel (m, n), we separately calculate the electric-field values along the x-and y-axis by: and: where E x (V·m −1 ) and E y (V·m −1 ) are the electric-field components along the x-and y-axis, and h x (m) and h y (m) are the distance between the voxel nodes along the x-and y-axis, respectively. Then, the combined electric-field value E x,y (V·m −1 ) is calculated by: and: where Ex (V·m −1 ) and Ey (V·m −1 ) are the electric-field components along the x-and y-axis, and hx (m) and hy (m) are the distance between the voxel nodes along the x-and y-axis, respectively. Then, the combined electric-field value Ex,y (V·m −1 ) is calculated by: Figure A1. A uniform rectangular grid with voxel (m, n) at its center. Figure A1. A uniform rectangular grid with voxel (m, n) at its center.